BdG 格点体系中Josephson键电流的格林函数推导

基本设定

考虑一个一般的格点BdG哈密顿量

$\hat H = \frac12 \sum_{ij} \hat\Psi_i^\dagger H_{ij} \hat\Psi_j .$

这里的$i,j$可以表示格点，也可以表示沿输运方向的切片；$\hat\Psi_i$是Nambu旋量，最简单的情形下为

$\hat\Psi_i = \begin{pmatrix} \hat c_i\\ \hat c_i^\dagger \end{pmatrix}.$

如果还有自旋、轨道、横向格点等自由度，它们都包含在$\hat\Psi_i$的内部指标中。$H_{ij}$是BdG矩阵块，满足厄米性$H_{ji}=H_{ij}^\dagger$。Nambu空间中的电荷矩阵为

$\tau_z = \begin{pmatrix} 1 & 0\\ 0 & -1 \end{pmatrix},$

更一般的表示是

$\tau_z = \begin{pmatrix} I & 0\\ 0 & -I \end{pmatrix}.$

局域电荷算符可写成

$\hat Q_m=-\frac{e}{2}\hat{\Psi}^\dagger_m\tau_z\hat{\Psi}_m$

这里的$1/2$来自 Nambu 表示的双计数。很多 BdG 线性响应文献会把这个$1/2$与 Nambu 求迹的双计数约定合并处理，最后电流公式常写成没有显式$1/2$的形式。下面推导重点是矩阵结构和格林函数关系；整体$1/2$因子取决于 Nambu 规范，但物理电流不变。

连续性方程

局域电荷变化由 Heisenberg 方程给出

$\frac{d\hat Q_m}{dt} = \frac{i}{\hbar} [\hat H,\hat Q_m].$

我们只关系从格点$m$到相邻格点$n$的键电流，因此只保留哈密顿量中连接$m$和$n$的部分

$\hat H_{mn} = \hat\Psi_m^\dagger H_{mn}\hat\Psi_n + \hat\Psi_n^\dagger H_{nm}\hat\Psi_m .$

这里先省略BdG总哈密顿量中的整体$1/2$，把它理解为已经包含在最终的Nambu规范中。利用费米子双线性算符的対易关系，对于任意矩阵$A,B$有

$[ \hat\Psi_i^\dagger A\hat\Psi_j, \hat\Psi_l^\dagger B\hat\Psi_l ] = \delta_{jl} \hat\Psi_i^\dagger A B\hat\Psi_l - \delta_{il} \hat\Psi_l^\dagger B A\hat\Psi_j .$

现在选择$l=m,B=\tau_z$，第一项为

$[ \hat\Psi_m^\dagger H_{mn}\hat\Psi_n, \hat\Psi_m^\dagger \tau_z \hat\Psi_m ] = - \hat\Psi_m^\dagger \tau_z H_{mn} \hat\Psi_n .$

第二项为

$[ \hat\Psi_n^\dagger H_{nm}\hat\Psi_m, \hat\Psi_m^\dagger \tau_z \hat\Psi_m ] = \hat\Psi_n^\dagger H_{nm}\tau_z \hat\Psi_m .$

因此，键$m\leftrightarrow n$对局域电荷变化的贡献为

$\left. \frac{d\hat Q_m}{dt} \right|_{mn} = -\frac{e i}{\hbar} \left[ - \hat\Psi_m^\dagger \tau_z H_{mn} \hat\Psi_n + \hat\Psi_n^\dagger H_{nm}\tau_z \hat\Psi_m \right].$

连续性方程写作

$\frac{d\hat Q_m}{dt} + \sum_n \hat I_{m\rightarrow n} = \hat S_m .$

其中$\hat{S}_m$是局域配对项带来的凝聚源项。如果我们在没有跨键配对的正常态跃迁键上取电流，则键电流定义为

$\hat I_{m\rightarrow n} = -\left. \frac{d\hat Q_m}{dt} \right|_{mn}.$

从而得到

$\hat I_{m\rightarrow n} = -\frac{e}{i\hbar} \left[ \hat\Psi_m^\dagger \tau_z H_{mn} \hat\Psi_n - \hat\Psi_n^\dagger H_{nm}\tau_z \hat\Psi_m \right].\label{eq:k1}$

如果$H_{mn}$是普通的BdG跃迁块

$H_{mn} = \begin{pmatrix} T_{mn} & 0\\ 0 & -T_{mn}^* \end{pmatrix},$

它与$\tau_z$时対易的$[H_{mn},\tau_z]=0$，这时候可以将$\tau_z$放到求迹的统一位置。

定义等时的单粒子密度矩阵

$\rho_{ij} = \langle \hat\Psi_j^\dagger \hat\Psi_i \rangle .$

利用恒等式

$\langle \hat\Psi_a^\dagger A \hat\Psi_b \rangle = \operatorname{Tr} \left[ A\rho_{ba} \right],$

键电流$\eqref{eq:k1}$的期望值可以表示为

$I_{m\rightarrow n} = -\frac{e}{i\hbar} \operatorname{Tr} \left[ \tau_z \left( \rho_{nm}H_{mn} - H_{nm}\rho_{mn} \right) \right].\label{eq:k2}$

这就是最核心的公式。后面将会看到，它比Matsubara形式和非平衡形式都更加基本，后面所有格林函数的表达，都是在解决如何计算$\rho_{ij}$？

Matsubara频率形式

平衡态中定义 Matsubara 格林函数

$G_{ij}(\tau) = - \langle T_\tau \hat\Psi_i(\tau) \hat\Psi_j^\dagger(0) \rangle .$

当$\tau\to 0^-$时，由于虚时序排列有

$G_{ij}(0^-) = \langle \hat\Psi_j^\dagger \hat\Psi_i \rangle = \rho_{ij}.$

另一方面

$G_{ij}(\tau) = \frac{1}{\beta} \sum_{\omega_n} e^{-i\omega_n\tau} G_{ij}(i\omega_n),\quad \omega_n=(2n+1)\pi/\beta .$

令$\tau=0^-$得到

$\rho_{ij} = \frac{1}{\beta} \sum_{\omega_n} G_{ij}(i\omega_n).$

因此

$\rho_{nm} = \frac{1}{\beta} \sum_{\omega_n} G_{nm}(i\omega_n),\qquad \rho_{mn} = \frac{1}{\beta} \sum_{\omega_n} G_{mn}(i\omega_n).$

代入密度矩阵形式的键电流公式$\eqref{eq:k2}$得到

$I_{m\rightarrow n} = -\frac{e}{i\hbar} \operatorname{Tr} \left[ \tau_z \left( \frac1\beta \sum_{\omega_n} G_{nm}(i\omega_n)H_{mn} - H_{nm} \frac1\beta \sum_{\omega_n} G_{mn}(i\omega_n) \right) \right].$

整理之后为

${\color{blue} I_{m\rightarrow n} = -\frac{e}{i\hbar\beta} \sum_{\omega_n} \operatorname{Tr} \left[ \tau_z \left( G_{nm}(i\omega_n)H_{mn} - H_{nm}G_{mn}(i\omega_n) \right) \right]. }$

非平衡小于格林函数形式

非平衡 Keldysh 方法中，核心对象是小于格林函数

$G^<_{ij}(t,t') = i \langle \hat\Psi_j^\dagger(t') \hat\Psi_i(t) \rangle .$

令$t’=t$得到

$G^<_{ij}(t,t) = i \langle \hat\Psi_j^\dagger(t) \hat\Psi_i(t) \rangle = i\rho_{ij}(t).$

因此

$\rho_{ij}(t) = -iG^<_{ij}(t,t).$

如果系统达到稳态，那么格林函数只依赖于时间差

$G^<_{ij}(t,t')=G^<_{ij}(t-t')$

Fourier变换为

$G^<_{ij}(t-t') = \int\frac{dE}{2\pi} e^{-iE(t-t')} G^<_{ij}(E).$

令$t=t’$得到

$G^<_{ij}(t,t) = \int \frac{dE}{2\pi} G^<_{ij}(E).$

因此密度矩阵为

$\boxed{ \rho_{ij} = -i \int \frac{dE}{2\pi} G^<_{ij}(E). }$

代入键电流公式$\eqref{eq:k2}$得到

$I_{m\rightarrow n} = -\frac{e}{i\hbar} \left(-i\right) \int\frac{dE}{2\pi} \operatorname{Tr} \left[ \tau_z \left( G^<_{nm}(E)H_{mn} - H_{nm}G^<_{mn}(E) \right) \right].$

其中

$\rho_{nm} = -i \int\frac{dE}{2\pi} G^<_{nm}(E),\qquad \rho_{mn} = -i \int\frac{dE}{2\pi} G^<_{mn}(E),$

化简最终得到

${\color{blue} I_{m\rightarrow n} = \frac{e}{\hbar} \int\frac{dE}{2\pi} \operatorname{Tr} \left[ \tau_z \left( G^<_{nm}(E)H_{mn} - H_{nm}G^<_{mn}(E) \right) \right]. }\label{eq:k3}$

这就是非平衡格林函数的键电流公式。

上式中只是将电流写成了$G^<$，真正非平衡物理来自于$G^<$如果计算。在稳定的非平衡格林函数中，有Keldysh方程

$\boxed{ G^<(E) = G^r(E)\Sigma^<(E)G^a(E). }$

其中

$G^r(E) = \left[ E+i0^+ - H - \Sigma^r(E) \right]^{-1},\qquad G^a(E) = \left[ G^r(E) \right]^\dagger .$

这里的$\Sigma^<(E)$描述外部reservoir 向系统注入粒子的分布信息。若系统连接多个 reservoir，则

$\Sigma^<(E) = \sum_\alpha \Sigma_\alpha^<(E).$

对处于局域平衡的第$\alpha$个reservoir

$\Sigma_\alpha^<(E) = i f_\alpha(E)\Gamma_\alpha(E),$

其中

$\Gamma_\alpha(E) = i \left[ \Sigma_\alpha^r(E) - \Sigma_\alpha^a(E) \right].$

所以小于格林函数为

${ G^<(E) = G^r(E) \left[ \sum_\alpha i f_\alpha(E)\Gamma_\alpha(E) \right] G^a(E). }$

代入非平衡键电流公式$\eqref{eq:k3}$中得到

${\color{magenta} I_{m\rightarrow n} = \frac{e}{\hbar} \int \frac{dE}{2\pi} \operatorname{Tr} \left\{ \tau_z \left[ \left( G^r\Sigma^<G^a \right)_{nm} H_{mn} - H_{nm} \left( G^r\Sigma^<G^a \right)_{mn} \right] \right\}. }$

这就是非平衡形式，它的核心不是把$G$替换成了$G^r$，而是将平衡态密度矩阵

$\rho = \frac1\beta \sum_{\omega_n}G(i\omega_n)$

推广为

$\rho = -i \int \frac{dE}{2\pi} G^<(E),$

并且用$G^<
=
G^r\Sigma^<G^a$来反映非平衡分布。

平衡态极限

如果所有 reservoir 处于同一温度和同一化学势，则

$f_\alpha(E)=f(E).$

此时

$\Sigma^<(E) = i f(E)\Gamma(E),\qquad \Gamma(E) = i [ \Sigma^r(E)-\Sigma^a(E) ].$

因此

$\Sigma^<(E) = f(E) [ \Sigma^a(E)-\Sigma^r(E) ].$

通过Dyson方程可得恒等式

$G^a(E)-G^r(E) = G^r(E) [ \Sigma^a(E)-\Sigma^r(E)]G^a(E)$

所以

$G^<(E) = G^r(E)\Sigma^<(E)G^a(E) = f(E) [ G^a(E)-G^r(E) ].$

于是

$\rho_{ij} = -i \int \frac{dE}{2\pi} G^<_{ij}(E) = -i \int \frac{dE}{2\pi} f(E) [ G^a_{ij}(E)-G^r_{ij}(E) ].$

这就是平衡实频率格林函数的表达式。

另一方面，Matsubara形式给出

$\rho_{ij} = \frac1\beta \sum_{\omega_n} G_{ij}(i\omega_n).$

两者是等价的，因为谱表示给出

$G(z) = \int \frac{dE}{2\pi} \frac{ A(E) }{ z-E },\qquad A(E) = i [ G^r(E)-G^a(E) ].$

代入Matsbara求和

$\frac1\beta \sum_{\omega_n} G(i\omega_n) = \int \frac{dE}{2\pi} A(E) \left[ \frac1\beta \sum_{\omega_n} \frac{1}{i\omega_n-E} \right],\quad \frac1\beta \sum_{\omega_n} \frac{1}{i\omega_n-E} = f(E).$

从而得到

${\color{teal} \frac1\beta \sum_{\omega_n} G(i\omega_n) = \int \frac{dE}{2\pi} f(E) i [ G^r(E)-G^a(E) ]. }$

因此，Matsubara频率求和与平衡实频谱积分实际上是等价的。