平均场超导理论中的格林函数、序参量与粒子数方程

模型与平均场近似

考虑标准的吸引相互作用模型

$H=\sum_{\mathbf{k},\sigma}\epsilon_{\mathbf{k}}\, c_{\mathbf{k}\sigma}^\dagger c_{\mathbf{k}\sigma} -\frac{U}{N}\sum_{\mathbf{k},\mathbf{k}'} c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow}, \qquad U>0,$

其中$N$是晶格点数或者$\mathbf{k}$点的数量。我们在巨正则系综中进行研究

$K \equiv H-\mu \hat N = \sum_{\mathbf{k},\sigma}\xi_{\mathbf{k}}\, c_{\mathbf{k}\sigma}^\dagger c_{\mathbf{k}\sigma} -\frac{U}{N}\sum_{\mathbf{k},\mathbf{k}'} c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow},$

其中$\xi_{\mathbf{k}}=\epsilon_{\mathbf{k}}-\mu$。接下来对相互作用进行平均场分解

$c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow} \approx \langle c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger\rangle c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow} + c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger \langle c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow}\rangle - \langle c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger\rangle \langle c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow}\rangle.$

定义超导序参量

$\Delta \equiv -\frac{U}{N}\sum_{\mathbf{k}} \langle c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}\rangle.\label{eq:k3}$

于是得到平均场下的哈密顿量

$K_{\rm MF} = \sum_{\mathbf{k}} \Big[ \xi_{\mathbf{k}} \big( c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow} + c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}\downarrow} \big) + \Delta\, c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger + \Delta^*\, c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow} \Big] +\frac{N|\Delta|^2}{U}.$

最后一项$\frac{N|\Delta|^2}{U}$是平均场分解带来的常数项，在计算自由能以及自洽方程的时候是需要保留的。

Nambu表象与矩阵形式

定义Nambu旋量

$\Psi_{\mathbf{k}} = \begin{pmatrix} c_{\mathbf{k}\uparrow}\\[3pt] c_{-\mathbf{k}\downarrow}^\dagger \end{pmatrix}, \qquad \Psi_{\mathbf{k}}^\dagger = \begin{pmatrix} c_{\mathbf{k}\uparrow}^\dagger & c_{-\mathbf{k}\downarrow} \end{pmatrix}.$

并引入Nambu空间中的Pauli矩阵

$\tau_x= \begin{pmatrix} 0&1\\ 1&0 \end{pmatrix}, \quad \tau_y= \begin{pmatrix} 0&-i\\ i&0 \end{pmatrix}, \quad \tau_z= \begin{pmatrix} 1&0\\ 0&-1 \end{pmatrix}.$

接下来看两个双线性组合，首先直接计算

$\Psi_{\mathbf{k}}^\dagger \tau_z \Psi_{\mathbf{k}} = c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow} - c_{-\mathbf{k}\downarrow}c_{-\mathbf{k}\downarrow}^\dagger\label{eq:k1}$

利用费米子反对易关系$c c^\dagger = 1-c^\dagger c$得到

$\Psi_{\mathbf{k}}^\dagger \tau_z \Psi_{\mathbf{k}} = c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow} + c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}\downarrow} -1.$

因此

$c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow} + c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}\downarrow} = 1+\Psi_{\mathbf{k}}^\dagger \tau_z \Psi_{\mathbf{k}}.\label{eq:k2}$

这里说明在Nambu表象中$\tau_z$测量的是”粒子-空穴差”$\eqref{eq:k1}$，而不是简单的粒子数，真正的粒子数还会多出一个常数$\eqref{eq:k2}$。对于非对角项则有

$\Psi_{\mathbf{k}}^\dagger \tau_x \Psi_{\mathbf{k}} = c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger + c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow},\\ \Psi_{\mathbf{k}}^\dagger \tau_y \Psi_{\mathbf{k}} = -i\,c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger +i\,c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}.$

将序参量表示为$\Delta=\Delta_x+i\Delta_y$，则可以得到

$\Delta\, c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger +\Delta^*\, c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow} = \Psi^\dagger_{\mathbf{k}}(\Delta_x\tau_x-\Delta_y\tau_y)\Psi_{\mathbf{k}}$

因此平均场哈密顿量表示为

$K_{\rm MF} = \sum_{\mathbf{k}} \Psi_{\mathbf{k}}^\dagger \Big( \xi_{\mathbf{k}}\tau_z+\Delta_x\tau_x-\Delta_y\tau_y \Big) \Psi_{\mathbf{k}} +\sum_{\mathbf{k}}\xi_{\mathbf{k}} +\frac{N|\Delta|^2}{U}.$

通过一个全局的规范变换，可以将$\Delta$变为实数，此时$\Delta_y=0$。因此

$K_{\rm MF} = \sum_{\mathbf{k}} \Psi_{\mathbf{k}}^\dagger \big( \xi_{\mathbf{k}}\tau_z+\Delta\tau_x \big) \Psi_{\mathbf{k}} +\sum_{\mathbf{k}}\xi_{\mathbf{k}} +\frac{N\Delta^2}{U}.$

定义BdG矩阵

$\mathcal H_{\rm BdG}(\mathbf{k})=\xi_{\mathbf{k}}\tau_z+\Delta\tau_x.$

得到

$K_{\rm MF} = \sum_{\mathbf{k}}\Psi_{\mathbf{k}}^\dagger \mathcal H_{\rm BdG}(\mathbf{k})\Psi_{\mathbf{k}} +\sum_{\mathbf{k}}\xi_{\mathbf{k}} +\frac{N\Delta^2}{U}.$

对$\mathcal H_{\rm BdG}$对角化即可得到Bogoliubov 准粒子能谱

$\mathcal H_{\rm BdG}^2(\mathbf{k}) = (\xi_{\mathbf{k}}^2+\Delta^2)\,\mathbb 1\Longrightarrow E_{\mathbf{k}}^{(\pm)}=\pm E_{\mathbf{k}}, \qquad E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^2+\Delta^2}.$

这就是标准的 BCS 准粒子色散关系。超导打开的能隙本质上体现在$\Delta$的出现。

矩阵格林函数

定义 Matsubara 格林函数

$\hat G(\mathbf{k},\tau) = -\langle T_\tau \Psi_{\mathbf{k}}(\tau)\Psi_{\mathbf{k}}^\dagger(0)\rangle.$

在频率空间中有

$\hat G(\mathbf{k},i\omega_n) = \Big[i\omega_n-\mathcal H_{\rm BdG}(\mathbf{k})\Big]^{-1} = \Big[i\omega_n-\xi_{\mathbf{k}}\tau_z-\Delta\tau_x]^{-1}$

结合关系式

$(a\mathbb 1-\mathbf b\cdot\boldsymbol\tau)^{-1} = \frac{a\mathbb 1+\mathbf b\cdot\boldsymbol\tau}{a^2-|\mathbf b|^2},$

可以得到格林函数为

$\hat G(\mathbf{k},i\omega_n) = \frac{i\omega_n+\xi_{\mathbf{k}}\tau_z+\Delta\tau_x} {(i\omega_n)^2-\xi_{\mathbf{k}}^2-\Delta^2}\qquad 或者\qquad \hat G(\mathbf{k},i\omega_n) = -\frac{i\omega_n+\xi_{\mathbf{k}}\tau_z+\Delta\tau_x} {\omega_n^2+E_{\mathbf{k}}^2}.$

将其表示为矩阵形式

$\hat G(\mathbf{k},i\omega_n)= \begin{pmatrix} G(\mathbf{k},i\omega_n) & F(\mathbf{k},i\omega_n)\\ \bar F(\mathbf{k},i\omega_n) & \tilde G(\mathbf{k},i\omega_n) \end{pmatrix},$

则

$G(\mathbf{k},i\omega_n) = -\frac{i\omega_n+\xi_{\mathbf{k}}}{\omega_n^2+E_{\mathbf{k}}^2},\qquad \tilde G(\mathbf{k},i\omega_n) = -\frac{i\omega_n-\xi_{\mathbf{k}}}{\omega_n^2+E_{\mathbf{k}}^2},\\ F(\mathbf{k},i\omega_n) = -\frac{\Delta}{\omega_n^2+E_{\mathbf{k}}^2}, \qquad \bar F(\mathbf{k},i\omega_n) = -\frac{\Delta}{\omega_n^2+E_{\mathbf{k}}^2}.\label{eq:k4}$

其中$G$是正常态电子格林函数，$\tilde{G}$是正常态空穴格林函数，$F,\bar{F}$是反常格林函数，它非零正是超导配对存在的体现。结合序参量的定义$\eqref{eq:k3}$以及反常平均值与格林函数之间的关系

$\langle c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}\rangle = T\sum_n F(\mathbf{k},i\omega_n)$

因此

$\Delta =-\frac{UT}{N}\sum_{\mathbf{k},n}F(\mathbf{k},i\omega_n).$

代入方程$\eqref{eq:k4}$得到

$\Delta=\frac{UT}{N}\sum_{\mathbf{k},n}\frac{\Delta}{\omega_n^2+E_{\mathbf{k}}^2}.$

考虑序参量$\Delta\neq 0$，两边同时约掉从而有

$\frac{1}{U}=\frac{T}{N}\sum_{\mathbf{k},n} \frac{1}{\omega_n^2+E_{\mathbf{k}}^2}.$

结合松原频率求和

$T\sum_n \frac{1}{\omega_n^2+E^2} = \frac{1}{2E}\tanh\frac{\beta E}{2},\label{eq:k5}$

从而得到有限温能隙方程

$\boxed{ \frac{1}{U} = \frac{1}{N}\sum_{\mathbf{k}} \frac{1}{2E_{\mathbf{k}}} \tanh\frac{\beta E_{\mathbf{k}}}{2} }.$

考虑零温极限

$T\to 0,\qquad \tanh\frac{\beta E_{\mathbf{k}}}{2}\to 1,$

从而得到标准的BCS自洽能隙方程

$\boxed{ \frac{1}{U} = \frac{1}{N}\sum_{\mathbf{k}}\frac{1}{2E_{\mathbf{k}}} }.$

粒子数摸底定义为

$n = \frac{1}{N}\sum_{\mathbf{k},\sigma} \langle c_{\mathbf{k}\sigma}^\dagger c_{\mathbf{k}\sigma}\rangle.$

将$\mathbf{k}\uparrow$和$-\mathbf{k}\downarrow$配成一组可以写成

$n = \frac{1}{N}\sum_{\mathbf{k}} \left\langle c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow} + c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}\downarrow} \right\rangle.$

结合等式$\eqref{eq:k2}$

$c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow} + c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}\downarrow} = 1+\Psi_{\mathbf{k}}^\dagger \tau_z \Psi_{\mathbf{k}},\nonumber$

可得

$n = 1+\frac{1}{N}\sum_{\mathbf{k}} \langle \Psi_{\mathbf{k}}^\dagger \tau_z \Psi_{\mathbf{k}}\rangle.$

因此在Nambu表象中，下分量是空穴算符$c^\dagger$，因此粒子数算符不是简单的$\Psi^\dagger\Psi$，而是要通过$\tau_z$来投影粒子与空穴，并且自动带出一个常数。

进一步将双线性平均值写成 Matsubara 格林函数

$\langle \Psi_{\mathbf{k}}^\dagger \tau_z \Psi_{\mathbf{k}}\rangle = T\sum_n {\rm Tr}\big[\tau_z \hat G(\mathbf{k},i\omega_n)\big].$

所以有

$n = 1+\frac{T}{N}\sum_{\mathbf{k},n} {\rm Tr}\big[\tau_z \hat G(\mathbf{k},i\omega_n)\big].$

代入格林函数

$\hat G(\mathbf{k},i\omega_n) = -\frac{i\omega_n+\xi_{\mathbf{k}}\tau_z+\Delta\tau_x} {\omega_n^2+E_{\mathbf{k}}^2},\nonumber$

计算

${\rm Tr}\big[\tau_z \hat G(\mathbf{k},i\omega_n)\big]=-\frac{ {\rm Tr}(i\omega_n\tau_z+\xi_{\mathbf{k}}\tau_z^2+\Delta\tau_z\tau_x)}{\omega_n^2+E_{\mathbf{k}}^2}.$

结合关系式

${\rm Tr}\,\tau_z=0,\qquad {\rm Tr}(\tau_z\tau_x)=0,\qquad {\rm Tr}\,\tau_z^2=2,$

得到

${\rm Tr}\big[\tau_z \hat G(\mathbf{k},i\omega_n)\big] = -\frac{2\xi_{\mathbf{k}}}{\omega_n^2+E_{\mathbf{k}}^2}.$

所以粒子数为

$n = 1-\frac{2T}{N}\sum_{\mathbf{k},n} \frac{\xi_{\mathbf{k}}}{\omega_n^2+E_{\mathbf{k}}^2}.$

同样结合方程$\eqref{eq:k5}$，得到粒子数为

$n = 1-\frac{1}{N}\sum_{\mathbf{k}} \frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh\frac{\beta E_{\mathbf{k}}}{2} .$

零温下为

$n = \frac{1}{N}\sum_{\mathbf{k}} \left( 1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}} \right) .$

从而就得到了超导平均场平衡态中的粒子数期望值，也就是超导态的平均粒子数密度。

在Bogoliubov变换中

$u_{\mathbf{k}}^2=\frac12\left(1+\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\right), \qquad v_{\mathbf{k}}^2=\frac12\left(1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\right).$

零温下的占据数满足

$n = \frac{2}{N}\sum_{\mathbf{k}}v_{\mathbf{k}}^2,$

因为每个动量$\mathbf{k}$上两个自旋的总占据正好是$2v_{\mathbf{k}}^2$，代入之后得到

$n = \frac{1}{N}\sum_{\mathbf{k}} \left(1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\right).$

因此，Nambu 矩阵格林函数方法与通常的 Bogoliubov 波函数方法是完全等价的，只是前者更适合推广到多轨道、多带以及更一般的矩阵问题。至此，我们就得到了平均场超导理论中的两个核心自洽方程

$\frac{1}{U} = \frac{1}{N}\sum_{\mathbf{k}} \frac{1}{2E_{\mathbf{k}}} \tanh\frac{\beta E_{\mathbf{k}}}{2},\\ n = 1-\frac{1}{N}\sum_{\mathbf{k}} \frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh\frac{\beta E_{\mathbf{k}}}{2}, \qquad E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^2+\Delta^2}.$

其中第一个方程决定序参量$\Delta$，第二个方程则根据填充数决定化学势$\mu$，两个方程联立则给出相互作用$U$、温度$T$以及粒子数密度$n$条件下的平衡态。

热力学方法推导

BdG方程对角化之后，准粒子能量为$E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^2+\Delta^2}$，平均场热力学势可以写成

$\Omega_{\rm MF} = \frac{N\Delta^2}{U} +\sum_{\mathbf{k}}(\xi_{\mathbf{k}}-E_{\mathbf{k}}) -\frac{2}{\beta}\sum_{\mathbf{k}}\ln\!\big(1+e^{-\beta E_{\mathbf{k}}}\big).$

由平衡态条件

$\frac{\partial \Omega_{\rm MF}}{\partial \Delta}=0$

可得

$\frac{2N\Delta}{U} -\sum_{\mathbf{k}}\frac{\Delta}{E_{\mathbf{k}}} \bigl[1-2f(E_{\mathbf{k}})\bigr]=0.$

结合

$1-2f(E)=\tanh\frac{\beta E}{2},$

得到自洽方程

$\frac1U = \frac1N\sum_{\mathbf{k}} \frac{1}{2E_{\mathbf{k}}} \tanh\frac{\beta E_{\mathbf{k}}}{2}.$

粒子数的定义为

$n=-\frac1N\frac{\partial \Omega_{\rm MF}}{\partial \mu}$

计算可得

$n = 1-\frac1N\sum_{\mathbf{k}} \frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh\frac{\beta E_{\mathbf{k}}}{2}.$

这和上面从矩阵格林函数得到的结果完全一致。

这里的粒子数方程给出的不是正常态粒子数，而是超导态自洽解对应的平均粒子数。而且它不是“粒子数本征值”，而是期望值，因为在BCS平均场中配对项
$\Delta\, c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger+\text{h.c.}$
破坏了$U(1)$规范对称性，所以平均场波函数一般不是粒子数算符$\hat{N}$的本征态。也就是超导基态满足
$\hat N |\Psi_{\rm BCS}\rangle \neq N_0 |\Psi_{\rm BCS}\rangle,$
但是有
$\langle \Psi_{\rm BCS}|\hat N|\Psi_{\rm BCS}\rangle = \text{某个固定值}.$
所以粒子数方程求出来的是超导态中的平均粒子数，不是严格投影到固定粒子数后的本征值。