表面格林函数求解
背景
在凝聚态物理中，特别是输运计算时，经常会遇到很多求解表面格林函数的问题。本文将简单介绍表面格林函数的求解。
考虑一个半无限长的体系，其哈密顿量可以写成：
\[H=\begin{bmatrix}
H_{00} & H_{01} & & & & &\\
H_{01} & H_{00} & H_{01} & & & &\\
& H_{01} & H_{00} & H_{01} & & &\\
& & \cdots & \cdots &\cdots& &\\
& & & \cdots & \cdots &\cdots&
\end{bmatrix}
\]
其中，\(H_{00}\)表示单层的哈密顿量，\(H_{01}\) 表示层间的相互作用哈密顿量，这里只考虑最近邻。
其推迟格林函数为 \((\omega + i\eta-H)G^{R} = I\), 可写成:
\[G^{R}=\begin{bmatrix}
\omega + i\eta - H_{00} & -H_{01} & & & & &\\
-H_{01} & \omega + i\eta - H_{00} & -H_{01} & & & &\\
& -H_{01} & \omega + i\eta - H_{00} & -H_{01} & & &\\
& & \cdots & \cdots &\cdots& &\\
& & & \cdots & \cdots &\cdots&
\end{bmatrix} ^{-1}
\]
注意，我们目的是求解上式中的 \(G^{R}(1,1)\), 即表面格林函数。显然，对于无限大推迟函数矩阵求逆不现实，需要考虑其他的方法。
简单解析解
我们先考虑最简单的一维无限长链，\(H_{00}\) 和 \(H_{01}\) 均为一维，设为 $H_{00}=\varepsilon, H_{01}=t $, 即推迟哈密顿量为:
\[G^{R}=\begin{bmatrix}
\omega + i\eta - \varepsilon & -t & & & & &\\
-t & \omega + i\eta - \varepsilon & -t & & & &\\
& -t & \omega + i\eta - \varepsilon & -t & & &\\
& & \cdots & \cdots &\cdots& &\\
& & & \cdots & \cdots &\cdots&
\end{bmatrix} ^{-1}
\]
由矩阵知识可知，逆矩阵等于伴随矩阵除以原矩阵的行列式，具体细节为:
\[D_{N}=\begin{vmatrix}
\omega + i\eta - \varepsilon & -t & & & & &\\
-t & \omega + i\eta - \varepsilon & -t & & & &\\
& -t & \omega + i\eta - \varepsilon & -t & & &\\
& & \cdots & \cdots &\cdots& &\\
& & & \cdots & \cdots &\cdots&
\end{vmatrix}_{N\times N}
\]
那么根据逆矩阵元 \(G_{ij}\) 等于伴随矩阵行列式\(|M_{ij}|\) 除以原矩阵行列式 \(|G|\), 即 \(G^{R}(1,1) = D_{N-1}/D_{N}\), 又因为对于上述行列式 \(D_{N}\)存在如下关系:
\[D_{N} = (\omega+i\eta-\varepsilon)D_{N-1} - t^2 D_{N-2}
\]
可写成: \(D_{N}/D_{N-1} = (\omega+i\eta-\varepsilon) - t^2 D_{N-2}/D_{N-1}\).
对于无限大的矩阵, 近似假设 \(D_{N}/D_{N-1} = D_{N-1}/D_{N-2}\), 那么上式转化为求解一元二次方程 \(t^2x^2-(\omega+i\eta-\varepsilon)x + 1=0\), 其解为:
\[x =G^{R}(1,1)= \frac{(\omega+i\eta-\varepsilon) \pm \sqrt{(\omega+i\eta-\varepsilon)^2 - 4t^2}}{2t^2}.
\]
相关测试代码如下
点击查看代码
clear all
N = 1000;
H00 = 0.0;
H01 = 0.5;
w = 0.5;
eta = 1e-10;
Ham = zeros(N,N);
for i = 1:N
Ham(i,i) = H00;
if i < N
Ham(i,i+1) = H01;
Ham(i+1,i) = H01;
end
end
Grinv = (w + 1j*eta)*eye(N) - Ham;
Gr = inv(Grinv);
Gr(1,1)
% 解析解
analy = ((w+1j*eta-H00) - 1j*sqrt(4*H01^2-(w+1j*eta-H00)^2))/(2*H01^2)
%analy1 = ((w+1j*eta-H00) + sqrt((w+1j*eta-H00)^2 - 4*H01^2))/(2*H01^2)
%analy2 = ((w+1j*eta-H00) - sqrt((w+1j*eta-H00)^2 - 4*H01^2))/(2*H01^2)
% Output Result Examples:
% H00 = 0; H01 = 0.5; w = 1.1; Gr(1,1) = 1.2835; analy = 1.2835 - 0.0000i
% H00 = 0; H01 = 0.5; w = 0.5; Gr(1,1) = 2.0000; analy = 1.0000 - 1.7321i
快速迭代求解
迭代求解表面格林函数方法，详细参考这篇论文 M P Lopez Sancho et al. J. Phys. F: Met. Phys. 15(1985) 851-858. http://iopscience.iop.org/0305-4608/15/4/009.
\[\begin{pmatrix}
\omega-H_{00} & -H_{01} & & & & \\
-H_{01}^{\dagger} & \omega-H_{00} & -H_{01} & & & \\
& -H_{01}^{\dagger} & \omega- H_{00} & -H_{01} & & \\
& & \cdots & \cdots &\cdots& \\
\end{pmatrix}
\begin{pmatrix}
G_{00} & G_{01} & G_{02} & \cdots\\
G_{10} & G_{11} & G_{12} & \cdots\\
G_{20} & G_{21} & G_{22} & \cdots \\
\vdots & \vdots & \vdots & \ddots\\
\end{pmatrix} = I
\]
有着如下关系:
\[(\omega-H_{00})G_{00} = I + H_{01}G_{10}\\
(\omega-H_{00})G_{10} = H_{01}^{\dagger}G_{00} + H_{01}G_{20}\\
\vdots\\
(\omega-H_{00})G_{n0} = H_{01}^{\dagger}G_{n-1,0} + H_{01}G_{n+1,0}
\]
最后可得如下关系:
\[G_{2n,0} = (\omega-H_{00})^{-1}(H_{01}^{\dagger}G_{2n-1,0} + H_{01}G_{2n+1,0})
\]
当 \(n>0\) 时：
\[[\omega-H_{00}-H_{01}(\omega-H_{00})^{-1}H_{01}^{\dagger} - H_{01}^{\dagger}(\omega-H_{00})^{-1}H_{01}]G_{2n,0}\\
=H_{01}^{\dagger}(\omega-H_{00})^{-1}H_{01}^{\dagger}G_{2n-2,0} + H_{01}(\omega-H_{00})^{-1}H_{01}G_{2n+2,0}
\]
当 \(n=0\) 时:
\[[\omega-H_{00}-H_{01}(\omega-H_{00})^{-1}H_{01}^{\dagger}]G_{00}\\
= I + H_{01}(\omega-H_{00})^{-1}H_{01}G_{20}
\]
将上面两式的因子重新标记为:
\[\alpha_1 = H_{01}(\omega-H_{00})^{-1}H_{01}\\
\beta_1 = H_{01}^{\dagger}(\omega-H_{00})^{-1}H_{01}^{\dagger}\\
\varepsilon^{s}_{1} = H_{00} + H_{01}(\omega-H_{00})^{-1}H_{01}^{\dagger}\\
\varepsilon_{1} = H_{00} + H_{01}(\omega-H_{00})^{-1}H_{01}^{\dagger} + H_{01}^{\dagger}(\omega-H_{00})^{-1}H_{01}
\]
那么可以得到:
\[(\omega-\varepsilon_{1}^{s})G_{00} = I + \alpha_1G_{20}\\
(\omega-\varepsilon_1)G_{2n,0} = \beta_1G_{2n-2,0} + \alpha_1G_{2n+2,0}
\]
这些参数的迭代关系为:
\[\alpha_i = \alpha_{i-1}(\omega-\varepsilon_{i-1})^{-1}\alpha_{i-1}\\
\beta_i = \beta_{i-1}(\omega-\varepsilon_{i-1})^{-1}\beta_{i-1}\\
\varepsilon_{i} = \varepsilon_{i-1} + \alpha_{i-1}(\omega-\varepsilon_{i-1})^{-1}\beta_{i-1} + \beta_{i-1}(\omega-\varepsilon_{i-1})^{-1}\alpha_{i-1}\\
\varepsilon^{s}_{i} = \varepsilon^{s}_{i-1} + \alpha_{i-1}(\omega-\varepsilon_{i-1})^{-1}\beta_{i-1}
\]
那么:
\[(\omega-\varepsilon_{i}^{s})G_{00} = I + \alpha_iG_{2^i,0}\\
(\omega-\varepsilon_{i})G_{2^in,0} = \beta_iG_{2^i(n-1),0} + \alpha_iG_{2^i(n+1),0}\\
\]
不断迭代后， \(\alpha\) 和 \(\beta\) 将逐渐减小, 取截断可以得到:
\[(\omega-\varepsilon^{s}_{final})G_{00} = I\\
G_{00} = (\omega-\varepsilon^{s}_{final})^{-1}.
\]
WannierTools 中的 Surfgreen.f90 核心代码如下:
快速迭代测试代码如下:
点击查看代码
% Quick iterative method
% Ref: M P Lopez Sancho et al. J. Phys. F: Met. Phys. 15(1985) 851-858.
clear all
N = 100;
H00 = 0.0;
H01 = 0.5;
w = 1.1;
eta = 1e-10;
Ham = zeros(N,N);
% Initial Data
alpha1 = H01*(w+1j*eta-H00)^(-1)*H01;
beta1 = H01*(w+1j*eta-H00)^(-1)*H01;
e1s = H00 + H01*(w+1j*eta-H00)^(-1)*H01;
e1 = H00 + H01*(w+1j*eta-H00)^(-1)*H01 + H01*(w+1j*eta-H00)^(-1)*H01;
alphai_1 = alpha1;
betai_1 = beta1;
ei_1s = e1s;
ei_1 = e1;
for i = 1:100
alphai = alphai_1*(w+1j*eta-ei_1)^(-1)*alphai_1;
betai = alphai_1*(w+1j*eta-ei_1)^(-1)*betai_1;
ei = ei_1 + alphai_1*(w+1j*eta-ei_1)^(-1)*betai_1 + betai_1*(w+1j*eta-ei_1)^(-1)*alphai_1;
eis = ei_1s + alphai_1 *(w+1j*eta-ei_1)^(-1)*betai_1;
alphai_1 = alphai;
betai_1 = betai;
ei_1s = eis;
ei_1 = ei;
end
Gs = (w+1j*eta-eis)^(-1)
% Output Result Examples:
% H00 = 0; H01 = 0.5; w = 1.1; Gs = 1.2835 - 0.0000i
% H00 = 0; H01 = 0.5; w = 0.5; Gs = 1.0000 - 1.7321i
与Quantum Transport 书中的P214 结果对比，点击查看代码
clear all;
Em_1 = 0;
Em_nx = 0.5;
N_E = 1200;
E = linspace(-1,3,N_E);
Sigma1 = zeros(N_E,1);
Sigma2 = zeros(N_E,1);
Sigma3 = zeros(N_E,1);
Sigma4 = zeros(N_E,1);
tx = 0.5;
for i_E = 1:N_E
GsL = surface_green_function(Em_1+2*tx,-tx,E(i_E),1e-10);
GsR = surface_green_function(Em_nx+2*tx,-tx,E(i_E),1e-10);
Sigma1(i_E) = tx*GsL*tx;
Sigma2(i_E) = tx*GsR*tx;
E_ratio1 = 1 - ( E(i_E) - Em_1 )/(2*tx);
ka1 = acos( E_ratio1 +2*(E_ratio1 < -2) +2*(E_ratio1 < -4) );
Sigma3(i_E) = -tx*exp(1j*ka1);
E_ratio2 = 1 - ( E(i_E) - Em_nx )/(2*tx);
ka2 = acos( E_ratio2 +2*(E_ratio2 < -2) +2*(E_ratio2 < -4) );
Sigma4(i_E) = -tx*exp(1j*ka2);
end
figure;
plot(E,real(Sigma1),'k.-',E,imag(Sigma1),'r.-')
figure;
plot(E,real(Sigma3),'k.-',E,imag(Sigma3),'r.-')
figure;
plot(E,Sigma2,'k.-')
figure;
plot(E,Sigma4,'r.-')
浙公网安备 33010602011771号