单极对力下 BCS 方程求解
作者：互联网
BCS方程可用于描述单满壳原子核的基态，这个笔记总结BCS方程的构成和特征，并展示BCS方程数值解的算法，讨论算法中的收敛性问题。
最近实验学家测量了\({}^{136}\rm Xe\)的价质子在不同轨道上的占据数，所以我借此做了实例演示，在唯象的单极对力哈密顿量下解BCS方程，与最新的实验数据相比较。
1. BCS波函数
核子-核子相互作用非常复杂，所以多体波函数也很复杂。
然而，由于核子-核子相互作用基本体现为短程吸引，在单满壳原子核中，对力是剩余相互作用的一个很好的近似。
在对力的作用下，原子核基态上同类核子配成\(M=0\)的对。
所以，单满壳原子核的基态波函数是可以大幅简化的，假设它的形式为
这里\(u_j, v_j\)是与每个单\(j\)轨道相关的变量。
这个波函数中，在每一对轨道\((j,m)\) 和 \((j,-m)\)上，要么没有粒子，要么占据2个核子.
我们定义
可以证明， BCS 波函数可以写作
\[\psi = \prod_{j} u_j^{\Omega_j} \exp{\frac{v_j}{u_j} S_+(j)}, \]其中 \(\Omega_j=\frac{2j+1}{2}\). \(\psi\) 的总角动量为 0.
然而，如果我们展开\(\psi\)，会得到一系列表达式的和，有些表达式有0个粒子(\(\sim|0\rangle\)), 有些有2个粒子(\(\sim| a^\dagger_{jm} a^\dagger_{j-m} \rangle\)), 或者4个粒子......
所以，BCS方程的粒子数是不守恒的，这是广为人知的。
因此，我们需要限制\(v_j\),使得\(\psi\)的粒子数期望值等于真实粒子数，即
\(|\psi \rangle\)的归一化也提供了限制条件，即
\[\langle 0 | \{ u_j + (-1)^{j-m} v_j a_{j -m} a_{jm} \} \{ u_j + (-1)^{j-m} v_j a^\dagger_{j m} a^\dagger_{j -m} \} | 0 \rangle = (u^2_j + v^2_j), \]得到
\[\langle \psi | \psi \rangle = \prod_j (u^2_j + v^2_j)^{\Omega_j} = 1. \]所以我们令
\[u^2_j + v^2_j \equiv 1 \]这样归一化条件会自然满足。实际上\(u_j,v_j\)都可设为正数，进一步简化问题。
2. \(u_j,v_j\)的物理含义
在\(\psi\)中，有多少个粒子占据\((j',m')\)？我们计算相应的粒子数算符\(a^\dagger_{j'm'} a_{j'm'}\)的期望值，得到
\[\langle \psi | a^\dagger_{j'm'} a_{j'm'} | \psi \rangle = v^2_{j'}, \]请注意它不依赖\(m'\)。所以，\(v^2_j\)是\((j,m)\)被占据的概率，其中\(m\)可以是\(-j\)到\(j\)之间的任意数。
这是合理的，因为我们假设的BCS波函数对于不同的\(m\)值是对称的。这也体现在 BCS 波函数角动量是 0 这一点。
由于\(u^2_j + v^2_j =1\)，所以\(u^2_j\)的意义是\((j,m)\)没有被占据的概率。
我们全部的任务，就是在特定的哈密顿量下，解出\(v_j\)，这样就得到了原子核的费米面。
3. 单极对力下的BCS方程
单极对力哈密顿为
\[H = \sum_{jm} \epsilon_j a^\dagger_{jm} a_{jm} - G \sum_{j j'} S_+(j) S_-(j'), \]右边第一项是单粒子能，第二项是单极对力，G是对力的强度。
我们的目标是对\(v_j\)做变分，在粒子数条件限制之下，使得\(\langle \psi | H | \psi \rangle\)取得最小值。
我们使用拉格朗日乘子法，即求
\[L = \langle \psi | H | \psi \rangle - \lambda ( \langle \psi | \sum_j N_j | \psi \rangle ) \]的最小值，即
\[\frac{ \partial L } { \partial v_j } = 0, ~~~ \frac{ \partial L } { \partial \lambda } = 0. \]代入BCS波函数，即得BCS方程
\[2(\epsilon_j - \lambda) u_j v_j - \Delta ( 1 - 2v^2_j) - 2 G v^3_j u_j = 0, \\ \sum_j (2j+1) v^2_j = N, \]这里\(\Delta = G \sum_j \Omega_j u_j v_j\).
4. 数值求解
4.1 迭代逼近
在求解BCS方程之前，不妨简单讨论迭代逼近方法的收敛问题。
假设我们得到了一个关系式
并需要求得精确解，然而这个关系式可能是高阶的，非线性的。
我们可以提出一个初始试探解\(x_0\)，并定义\(x_1 = f(x_0), x_2 = f(x_1), \cdots\)，希望这个数列收敛在\(x_{\rm true} = f( x_{\rm true} )\)。
它会收敛吗？
我发现，如果\(f(x)\)单调且有界，情况就简单多了。
如果\(f(x)\)是单调的，\(f(f(x))\)就一定是单调递增的，所以数列\(x_0, x_2, \cdots x_{2k}, \cdots\)与数列\(x_1, x_3 \cdots x_{2k+1}, \cdots\)必为单调有界，所以它们必然各自存在极限\(x_{\rm even}\)与\(x_{\rm odd}\)。
如果\(x_{\rm even} = x_{\rm odd}\)我们即得到了\(x=f(x)\)的解，否则\(x_0,x_1,\cdots\)将在两个极限之间跳跃。
例如，定义在\((0,1)\)上的\(f(x) = \sqrt{1-x^2}\)，\(x_n\)将会在\(x_0\)和\(\sqrt{1-x^2_0}\)之间跳跃，迭代求解将会进入死循环。
但是如果\(f(x)\)是单调递增且有界的，\(x_0, x_1, x_2 \cdots\)一定会收敛。
4.2 BCS方程的数值解
在 BCS 方程中，\(\epsilon_j,G,\Omega_j\)都是设定的参数，所以\(v_j\)是\(\lambda\)的函数。
对某个特定的\(\lambda\)值，我们可以求得\(v_j\)，将它们代入 BCS 方程, 求得
改变\(\lambda\)的值，当\(g(\lambda)=0\)时，即得正确的\(\lambda\)值，这个值有时被叫做``化学势”。
对于给定的\(\lambda\)值，求得\(\Delta\)的算法是这样的：
- 1.提出一个初始试探值\(\Delta\)（比如6.0）；
- 2.用这个值计算\(v_j(\lambda,\Delta)\)；
- 3.将\(v_j\)值代入\(\Delta=G\sum_j u_j v_j\)重新计算\(\Delta\)；
- 4.用这个新的\(\Delta\)值再次计算\(v_j(\lambda,\Delta)\)，如此循环，直到\(\Delta\)收敛。
4.2.1 可以证明 \(\Delta\) 一定会收敛
\(v^2_j\)越靠近\(\frac{1}{2}\)，\(u_j v_j\)越大。
我们定义\(v_j = \cos \alpha, \alpha \in (0,\frac{\pi}{2})\)，把 BCS 方程写作
由于\(\Delta\)为正，\((\epsilon_j -\lambda - G\cos^2\alpha)\)与\(\tan 2\alpha\)必然同号。
\[\frac{ d\Delta}{d\alpha} = G \sin 2\alpha \tan 2\alpha + (\epsilon_j - \lambda - G \cos^2\alpha)\sec^2\alpha, \]如果\(v^2_j > \frac{1}{2}\)，\(\frac{d\Delta}{d\alpha}<0\)，所以如果\(\Delta\)增大，\(v^2_j\)减小（距离\(\frac{1}{2}\)更近)，所以\(u_j v_j\)增大，\(f_1(\Delta)\)增大。
如果\(v^2_j < \frac{1}{2}\)，\(\frac{d\Delta}{d\alpha}>0\)，所以如果\(\Delta\)增大，\(v^2_j\)减小（距离\(\frac{1}{2}\)更近），所以\(u_j v_j\)增大，\(f_1(\Delta)=G\sum_j \Omega_j u_j v_j\)增大。
所以\(f_1(\Delta)\)是单调递增的。
显然\(f_1(\Delta)\)有边界\(\left[ 0, G\sum_j \frac{1}{2} \Omega_j \right]\)，所以\(\Delta = f_1(\Delta)\)在迭代求解中一定会收敛。
4.2.2 可以证明 \(v_j\) 一定会收敛
根据 BCS 方程
\[v^2_j = f_2(v^2_j) = \frac{1}{2} \left[ 1 - \frac{ (\epsilon_j -\lambda - G v^2_j)}{\sqrt{(\epsilon_j - \lambda - G v^2_j)^2 + \Delta^2} } \right]. \]如果\(\lambda, \Delta\)已经设定，即可用\(v^2_j = f_2(v^2_j)\)迭代求数值解。
由于
所以\(f_2(v^2_j)\)是单调递增的，而\(v^2_j\)显然有界\(\left[0,1\right]\)。
所以\(v^2_j = f_2(v^2_j)\)必然收敛。
现在我们已经确定，对于任意给定\(\lambda\)值，都可以解出\(v^2_j\)和\(\Delta\)确定的数值解。
只需满足BCS方程的粒子数，即\(\langle\psi|\sum_j N_j | \psi \rangle=N\)。
所以我们定义\(g(\lambda) = \langle \psi | \sum_j N_j | \psi \rangle - N\)，并且寻找其数值解，使得\(g(\lambda) = 0\).
割线法可用于完成这个任务。
这样BCS方程即可解出。
5. 例子
\(^{136}\rm Xe\)是一个单满壳原子核，有4个价质子在核心\({}^{132}\rm Sn\)外。
最近实验学家通过（\(d\),\({}^3\rm He\)）反应测量了它的价质子在不同轨道上的占据率[2]，所以我们可以用BCS方程计算\({}^{136}\rm Xe\)的价质子占据率，并与实验数据比照。
我们考虑50-82大壳，加上一根\(0h_{9/2}\)轨道（其能量设为大壳间距，约\(8\rm MeV\)。
考虑对力\(G=0.18 \rm MeV\), \(^{136}\rm Xe\)，解得\(\lambda=-0.1854, \Delta=0.673049\)，\(v_j\)值与相应的占据值在表1中列出。
考虑到这个模型的简单性，我们的结果已经非常好了。
\(0h_{11/2}\)的重要性被低估了，其占据数为0.4而我们只得到0.174。
\(0h_{9/2}\)没有实验值，理论值接近于0是因为单极对力不足以导致跨壳激发。
| Orbits | \(\epsilon_j\)(MeV) | \(v_j\) | Occupancy(BCS) | Occupancy(EXP) |
|---|---|---|---|---|
| \(2s_{1/2}\) | 2.990 | 0.112 | 0.025 | 0.057(6) |
| \(1d_{5/2}\) | 2.690 | 0.123 | 0.061 | |
| \(1d_{3/2}\) | 0.963 | 0.277 | 0.461 | 0.52(3) (\(1d_{5/2} + 1d_{3/2}\)) |
| \(0g_{7/2}\) | 0.000 | 0.638 | 3.259 | 2.93(10) |
| \(0h_{11/2}\) | 2.760 | 0.120 | 0.174 | 0.40(4) |
| \(0h_{9/2}\) | 8.0 | 0.044 | 0.020 |
表1. 136Xe_occupancies \(^{136}\rm Xe\)的质子占据数。实验值从参考文献[2]中得到。
[1] R. D. Lawson, ``Theory of the nuclear shell model", Clarendon Press \(\cdot\) Oxford 1980, p339-p344.
[2] J. P. Entwisle, B. P. Kay, A. Tamii et al., Phys. Rev. C 93, 064312 (2016).
标签：psi,frac,对力,Delta,rm,BCS,单极,lambda 来源： https://www.cnblogs.com/luyi07/p/16580489.html