投影技术之线性代数投影：LAP

之前写了投影技术的理论基础，即转动算符、D矩阵、投影算符等理论知识：https://www.cnblogs.com/luyi07/p/14586631.html
又写了两种投影技术的思路与比较：https://www.cnblogs.com/luyi07/p/14823838.html
在这个随笔中，我整理其中一种思路，即线性代数投影（Linear Algebraic Projection: LAP) 的操作细节。以下是参考文献：
[1] Calvin W. Johnson and Keven D. O'Mara, "Projection of angular momentum via linear algebra", Physical Review C 96, 064304 (2017).
[2] Calvin W. Johnson and Changfeng Jiao "Convergence and efficiency of angular momentum projection", Journal of Physics G: Nuclear and Particle Physics 46, 015101 (2019).

1. LAP 投影

LAP 方法的核心是下面的式子，其中\(N^J_{KM} = \langle \Psi | \hat{P}^J_{KM} | \Psi' \rangle\)。

\[\langle \Psi | \hat{R}(\Omega_i) | \Psi' \rangle = \sum_{JKM} D^J_{KM}(\Omega_i) N^J_{KM}. \tag{1} \]

需要选取合适的 \(\Omega_i\)，使得 \(D^J_{KM}(\Omega_i)\) 是一个可逆的矩阵，最好还将式子稍作变形，变成实数域的操作，方便编程。

2. 欧拉角 \(\alpha, \gamma\) 的相关逆运算

2.1 单位圆上复数求和式

如果 \(M = 0, 1, \cdots, N-1\)，则

\[\frac{1}{N} \sum^{N-1}_{k=0} exp( i \frac{2\pi M }{N}k ) = \delta_{M,0}. \]

记 \(\theta = \frac{2\pi M }{N}\)，则其取值在 \([0, 2\pi)\)。上式左边为复平面单位圆上一系列等间距的点。很容易证明上式，即 \(M\neq0\) 时求和为 0，否则为1。

\[S = \frac{1}{N} \sum^{N-1}_{k=0} exp( i k \theta ),\\ e^{i\theta} S = \frac{1}{N} \sum^{N}_{k=1} exp( i k \theta ), \\ ( e^{i\theta} -1 ) S = \frac{1}{N} (e^{iN\theta} - 1 ) = 0. \]

所以，\(M = 1,2,\cdots, N-1\) 时，\(S = 0\)；\(M=0\) 时，则直接计算 \(S = \frac{1}{N} N = 1\)。

顺便多说一句，\(\{ 1, e^{i\theta}, \cdots, e^{i(N-1)\theta} \}\)是 \(C_N\) 群的表示。一共有\(N\)个类，所以正好构成\(N\)个不可约表示，这个求和式就是这些不可约表示的正交归一关系。

2.2 \(\alpha, \gamma\) 相关的矩阵求逆

以 \(\gamma\) 为例，如果要投影的最大角动量为 \(J_{max}\)，则取 \(\gamma\) 的离散值为

\[\gamma_k = k \frac{ 2\pi }{ 2J_{max} +1 }, k = 0, 1, \cdots, 2 J_{max}; \]

然后定义矩阵

\[Z_{K k} = \frac{1}{2J_{max}+1} e^{-i K \gamma_k }. \]

因为 \(D\) 矩阵为

\[D^{J}_{KM}(\alpha, \beta, \gamma) = e^{iK\alpha} d^J_{KM}(\beta) e^{iM\gamma}, \]

所以对方程 (1) 进行操作，得到

\[\sum_{ik} Z_{Ki} Z_{M k} \langle \Psi | \hat{R}( \alpha_i, \beta_j, \gamma_k ) | \Psi' \rangle = \sum_J d^J_{KM}(\beta_j) N^J_{KM}. \tag{2} \]

到了这里，\(j\) 的维数一般不等于 \(J\) 的维数，做逆变换不是特别方便，所以构造了下面的办法。

3. 关于 \(\beta\) 的矩阵求逆

其实可以看出来，因为\(D^J_{KM}\)可以写成分别与 \(\alpha, \beta, \gamma\) 相关的三部分的乘积，所以 \(D^J_{KM}(\Omega_i)\) 可以看做是三个小矩阵的直积。因此先对 \(\alpha, \gamma\) 相关的小矩阵做取逆运算，再做与 \(\beta\) 相关的取逆运算，是可行的。
取 \(\beta\) 的离散点

\[\beta_j = ( j - \frac{1}{2} )\frac{\pi}{N}, j=1,2,\cdots,N， \]

若为偶核子系统，则取 \(N = J_{max} + 1\)，否则取 \(N = J_{max} + \frac{1}{2}\)。构造

\[\Delta^{J' J}_{K M} = \sum_j d^{J'}_{KM} (\beta_j) d^J_{KM}(\beta_j), \]

这里 \(J, J' \geq |M|, |K|\)。
然后对方程 (2) 进行操作，得到

\[\sum_J \Delta^{J' J}_{KM} N^J_{KM} = \sum_{ ijk} d^{J'}_{KM}(\beta_j) Z_{Ki} Z_{Mk} \langle \Psi | \hat{R}(\alpha_i, \beta_j, \gamma_k) | \Psi \rangle. \]

这里有个好处，就是 \(\Delta^{J' J}_{KM}, N^J_{KM}\)还有等号右边整体，这三者都是实数，编程的时候做实数方程组求解即可。如此解出 \(N^J_{KM}\)，类似地解出 \(H^J_{KM}\)，然后求解 Hill-Wheeler 方程。