Dyson 方程求解 G_{1N}

研究背景

在非平衡态Green函数理论中，Transmission 的公式为:

\[T(\omega) = Tr[\Gamma^{L}G^{R}\Gamma^{R}G^{R\dagger}]. \]

这里的 \(\Gamma^L,\Gamma^R, G^{R}\) 分别为:

\[\Gamma^L=\begin{bmatrix} \Gamma^L & 0 & \cdots & \cdots \\ 0 & 0 & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & \cdots & 0 \end{bmatrix}, \Gamma^R=\begin{bmatrix} 0 & 0 & \cdots & \cdots \\ 0 & 0 & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & \cdots & \Gamma^R \end{bmatrix}, G^R=\begin{bmatrix} G_{11} & G_{12} & \cdots & \cdots \\ G_{21} & G_{22} & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ G_{n1} & \cdots & \cdots & G_{nn} \end{bmatrix} \]

对上式求trace展开, 得

\[\Gamma^{L}G^{R}\Gamma^{R}G^{R\dagger} = \begin{bmatrix} \Gamma^LG_{11} & \Gamma^LG_{12} & \cdots & \Gamma^LG_{1n} \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & \cdots & 0 \end{bmatrix}* \begin{bmatrix} 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & \cdots & 0 \\ \Gamma^RG_{1n}^{\dagger} & \Gamma^RG_{2n}^{\dagger} & \cdots & \Gamma^RG_{nn}^{\dagger} \end{bmatrix}\\ =\begin{bmatrix} \Gamma^LG_{1n}\Gamma^RG_{1n}^{\dagger} & 0 & \cdots & 0 \\ 0 & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix} \]

所以，\(Tr[\Gamma^{R}G^{R}\Gamma^{L}G^{R\dagger}]= Tr[\Gamma^LG_{1n}\Gamma^RG_{1n}^{\dagger}]\), 此时求trace公式转化为求解 \(G_{1n}\). 一般而言，对哈密顿量求逆即可得到 \(G\), 但是对于特别大的矩阵，求逆过程特别费时. 下面讲介绍 Dyson 迭代方程来快速求解 \(G_{1n}\).

理论推导

器件体系的 Green 函数如下所示,

\[G^{R}(\omega)=\begin{bmatrix} \omega - H_{11} - \Sigma^L & -H_{12} & & & & &\\ -H_{21} & \omega - H_{22} & -H_{23} & & & &\\ & \cdots & \cdots &\cdots& & &\\ & & -H_{n-1,n-2} & \omega - H_{n-1,n-1} & -H_{n-1,n} \\ & & & -H_{n,n-1} & \omega - H_{n,n} - \Sigma^R \end{bmatrix}^{-1} \]

其中 \(\omega+i\eta\) 省略写成 \(\omega\), \(\Gamma^{R} = i*(\Sigma^R-\Sigma^{R\dagger}), \Gamma^{L} = i*(\Sigma^L-\Sigma^{L\dagger})\).

Dyson 迭代方程为:

求解 \(G_{1n}\) 的迭代关系如下:

\[G_{11}^{(1)} = [\omega + i\eta - H_{11}-\Sigma^L]^{-1}\\ G_{12}^{(2)} = G_{11}^{(1)}H_{12}G_{22}^{(2)}\\ G_{13}^{(3)} = G_{12}^{(2)}H_{23}G_{33}^{(3)}\\ \cdots\\ G_{1n}^{(n)} = G_{1,n-1}^{(n-1)}H_{n-1,n}G_{n,n}^{(n)}\\ \]

其中,

\[G_{22}^{(2)} = [\omega + i\eta - H_{22} - H_{12}^{\dagger}G_{11}^{(1)}H_{12}]^{-1}\\ G_{33}^{(3)} = [\omega + i\eta - H_{33} - H_{23}^{\dagger}G_{22}^{(2)}H_{23}]^{-1} \cdots\\ G_{n-1,n-1}^{(n-1)} = [\omega + i\eta - H_{n-1,n-1} - H_{n-2,n-1}^{\dagger}G_{n-2,n-2}^{(n-2)}H_{n-2,n-1}]^{-1}\\ G_{nn}^{(n)} = [\omega + i\eta - H_{n,n} - H_{n-1,n}^{\dagger}G_{n-1,n-1}^{(n-1)}H_{n-1,n}- \Sigma^R]^{-1} \]

参考:
[1]. A. MacKinnon. The Calculation of Transport Properties and Density of States of Disordered Solids. Phys. B - Condensed Matter 59, 385-390 (1985).
[2]. https://www.guanjihuan.com/archives/7650

代码验证

Dyson 优化算法验证

点击查看代码

clear all;
tic
W = 10;
L = 30; 
Ham = zeros(W*L,W*L);
SigmaL = zeros(W*L,W*L);
SigmaR = zeros(W*L,W*L);
GammaL = zeros(W*L,W*L);
GammaR = zeros(W*L,W*L);
t = 1;
eta = 1e-10;
H00 = zeros(W,W);
H01 = zeros(W,W);
for i = 1:W
    H00(i,i) = 4*t;
    H01(i,i) = t;
    if i < W
        H00(i,i+1) = t;
        H00(i+1,i) = t;
    end
end
for j = 1:L
    Ham(j*W-W+1:j*W,j*W-W+1:j*W) = H00;
    if j < L
        Ham(j*W-W+1:j*W,j*W+1:j*W+W) = H01;
        Ham(j*W+1:j*W+W,j*W-W+1:j*W) = H01';
    end
end

N_omega = 300;
omega = linspace(0,1,N_omega);
Cond = zeros(N_omega,1);
for i = 1:N_omega
    Gs = surface_green_function(H00,H01,omega(i));
    SigmaL(1:W,1:W) = H01*Gs*H01';
    SigmaR(W*L-W+1:W*L, W*L-W+1:W*L) = H01*Gs*H01';
    GammaL = 1j*(SigmaL-SigmaL');
    GammaR = 1j*(SigmaR-SigmaR');
    Gr = inv((omega(i) + 1j*eta)*eye(W*L) - Ham - SigmaL - SigmaR); 
    Ga = Gr';
    Cond(i) = trace(GammaL*Gr*GammaR*Ga);
end
figure;
plot(omega,real(Cond));
xlabel('Energy (t)');
ylabel('Conducatance (e^2/h)')
toc

tic
for i = 1:N_omega
    Gs = surface_green_function(H00,H01,omega(i));
    SigmaL1 = H01*Gs*H01';
    SigmaR1 = H01*Gs*H01';
    G11 = inv((omega(i) + 1j*eta)*eye(W)- H00 - H01*Gs*H01');
    Gjj_j = G11;G1j_j = G11;
    for j = 2:L
       %G_{N+1,N+1}^{N+1} = inv(omega - H_{N+1,N+1} - H01'*G_{N,N}^{N}*H01)
       Gj1j1_j1 = inv((omega(i) + 1j*eta)*eye(W)- H00 - H01'*Gjj_j*H01);
	   if j == L
	       Gj1j1_j1 = inv((omega(i) + 1j*eta)*eye(W) - H00 - H01'*Gjj_j*H01 - H01*Gs*H01');
		end
       %G_{i,N+1}^{N+1} = G_{i,N}^{N}*V_{N}*G_{N+1,N+1}^{N+1}.
       G1j1_j1 = G1j_j*H01*Gj1j1_j1;
       Gjj_j = Gj1j1_j1;
       G1j_j = G1j1_j1;
    end
    G1L_L = G1j_j;
    GammaL1 = 1j*(SigmaL1-SigmaL1');GammaR1 = 1j*(SigmaR1-SigmaR1');
    Cond(i) = trace(GammaL1*G1L_L*GammaR1*G1L_L');
end
figure;
plot(omega,real(Cond));
xlabel('Energy (t)');
ylabel('Conducatance (e^2/h)')
toc

% Output
% 历时 13.909239 秒。
% 历时 1.748965 秒。

生成结果:

surface_green_function.m 为

点击查看代码

function Gs = surface_green_function(H00,H01,omega)
    eta = 1e-10;
    [m, n] = size(H00); 
    % Initial
    alphai_1 = H01*((omega+1j*eta)*eye(m)-H00)^(-1)*H01;
    betai_1 = H01'*((omega+1j*eta)*eye(m)-H00)^(-1)*H01';
    ei_1s = H00 + H01*((omega+1j*eta)*eye(m)-H00)^(-1)*H01';
    ei_1 = H00 + H01*((omega+1j*eta)*eye(m)-H00)^(-1)*H01' + H01'*((omega+1j*eta)*eye(m)-H00)^(-1)*H01;

    for i = 1:100
        alphai = alphai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*alphai_1;
        betai = betai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*betai_1;
        ei = ei_1 + alphai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*betai_1 + betai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*alphai_1;
        eis = ei_1s + alphai_1 *((omega+1j*eta)*eye(m)-ei_1)^(-1)*betai_1;
        % exit the loop
        if sum(sum(abs(alphai))) <= 1e-8
            break
        end
        % iterative
        alphai_1 = alphai;
        betai_1 = betai;
        ei_1s = eis;
        ei_1 = ei;  
    end
    Gs = ((omega+1j*eta)*eye(m)-eis)^(-1);
end