ナノエレクトロニクスなどの分野では複数の素材が接合されていたり、様々な端子がくっついていたり、あるいはデバイスの形状が電気特性に与える影響を考えることもある。こうした状況でデバイスの電気特性をシミュレーションするやり方がランダウアー・ビュティカー理論である。今回は再帰グリーン関数法を利用して計算するやり方を紹介したい。

この記事では物理学専攻の学部３年終了時点程度の知識を前提とする。

ランダウアー・ビュティカー理論

まず考えるシチュエーションを説明する。

考えるデバイスは図の緑で囲んだ部分であり、ハミルトニアン$H_c$で表現される。$H_C$は$N$個のサイトが１列に並んだものであり、それぞれのサイトはハミルトニアン$H_{i}$で表され(各$H_i$のサイズや中身はバラバラでもよい)、$V_{i,i+1}$で結合している($V_{i,i+1} = V^\dagger_{i+1,i}$)。
この左右両端がそれぞれハミルトニアン$H_{L/R}$で表される端子の端と重なり積分$V_{LC},V_{CR}$で結合している。
また、デバイス、端子の化学ポテンシャルは異なっていてもよくて、この化学ポテンシャルの違いがデバイスにかかる電圧になっている。これが電流などを流す原動力になっている。

全体のハミルトニアンはブロック行列として

H_{total} =\begin{pmatrix}
H_L & V_{L} & 0 \\
V_{L}^\dagger & H_C & V_{R} \\
0 & V_{R}^\dagger & H_R
\end{pmatrix}

と表される。
ハイゼンベルクの運動方程式と等価な式

(E-H_{total})G = I

ただし$G$は温度グリーン関数

G = \begin{pmatrix}
G_L & G_{LC} & G_{LR} \\
G_{CL} & G_C & G_{CR} \\
G_{RL} & G_{RC} & G_R
\end{pmatrix}

である。左辺を計算すると真ん中の列のブロック行列要素から

[ E - H_C - V\dagger_L(E-H_L)^{-1}V_L - V\dagger_R(E-H_R)^{-1}V_R]G_C = I

が導かれ、

\Sigma_{L/R}(E) = V\dagger_{L/R}(E-H_{L/R})^{-1}V_{L/R}

とすれば(この$\Sigma$を自己エネルギーという)、

G_C = (E-H_C-\Sigma_L(E)-\Sigma_R(E))^{-1}

となって$G_c$が表現できる。

ここで、左から右に流れる電流を計算するとMeir-Wingreen公式

\begin{align}
I &= \frac{e}{h}\int dE [f_L(E) - f_R(E) ] Tr[G_C^{(A)}(E)\Gamma_R(E)G_C^{(R)}(E)\Gamma_L(E)] \\
\Gamma_{L/R} &= -2 Im \Sigma^{(R)}_{L/R}
\end{align}

で表される(導出は割愛する)。$G,\Sigma$の上付き記号$(R),(A)$は遅延/先進グリーン関数であり、$i\omega_n \rightarrow \omega \pm i\delta$と置き換えたときの$\pm$の違いを表す。また、$f_L,f_r$は左右の端子の化学ポテンシャル$\mu_L,\mu_R$に対するフェルミ分布関数$\frac{1}{ e^{E-\mu } + 1 }$である。(左の端子、デバイス、右の端子で化学ポテンシャルが違っていてよいというのがポイントである)。

よって透過率

T(E) = Tr[G_C^{(A)}(E)\Gamma_R(E)G_C^{(R)}(E)\Gamma_L(E)]

を計算すればよい。

ここで、端子はシステムの両端のみと繋がっていて、

\begin{align}
H_C &= \begin{pmatrix}
H_1 & V_{12} & 0 & \cdots & 0 \\
V_{21} & H_2 & V_{23} & \cdots & 0 \\
\vdots & V_{32} & \ddots & & 0 \\
\vdots & & & \ddots & V_{N-1,N} \\
0 & \cdots & 0 & V_{N,N-1} & H_N
\end{pmatrix} \\

\Sigma_L &= \begin{pmatrix}
\Sigma_L & 0 & & 0 \\
0 & 0 & \cdots & 0 \\
\vdots &  & \ddots & \vdots \\
0 &\cdots & 0 & 0
\end{pmatrix} \\

\Sigma_R &= \begin{pmatrix}
0 & 0 & \cdots & 0 \\
0 & \ddots &  & \vdots \\
\vdots &  & 0 & 0 \\
0 &\cdots & 0 & \Sigma_R
\end{pmatrix}
\end{align}

だから

T(E) = Tr[G_C^{(A)}(E)\Gamma_R(E)G_C^{(R)}(E)\Gamma_L(E)] \\
= Tr[G_{1N}^{(A)}(E)\Gamma_R(E)G_{N1}^{(R)}(E)\Gamma_L(E)]

である。さらに$G_{1N}^{(A)} = [G_{N1}^{(R)}]^\dagger$だから$T(E)$は結局$\Sigma^{(R)}_{L/R}$と$G_{N1}^{(R)}$のブロック成分さえ分かればよい。

また、$\Gamma_{L/R}$は半正定値行列行列なので、コレスキー分解可能であり$\Gamma_{L/R} = X^\dagger X$となる行列が存在して、$\Gamma^{1/2}_{L/R}=X$とする。
そこで散乱行列$S$を用いて

\begin{align}
S &= i \Gamma_R^{1/2} G_{N1}^{(R)} \Gamma_L^{1/2} \\
T(E) &= Tr[S^\dagger S]
\end{align}

となる。そしてスピン流に対する透過率(スピン透過率)は

T_S(E) =  Tr[S^\dagger \sigma_i S] \quad (i=x,y,z)

と表される。

再帰グリーン関数法

以下、$G_C \rightarrow G$と表記する(面倒なので...)。
グリーン関数は

\begin{align}
G &= ( E - H_C - \Sigma_L - \Sigma_R )^{-1} \\

&= \begin{pmatrix}
E-H_1-\Sigma_L & -V_{12} & 0 & \cdots & 0 \\
-V_{21} & E-H_2 & -V_{23} & \cdots & 0 \\
0 & V_{32} & \ddots & & 0 \\
\vdots & & & \ddots & -V_{N-1,N} \\
0 & \cdots & 0 & -V_{N,N-1} & E - H_N - \Sigma_R
\end{pmatrix}^{-1} \\

&= \begin{pmatrix}
G_{11} & \cdots & G_{1N} \\
\vdots & \ddots & \vdots \\
G_{N1} & \cdots & G_{NN}
\end{pmatrix} \\
\end{align}

$H_C$の各ブロックで一番大きいものが$m\times m$行列だったとすると、まともに$G_C = (E-H_C-\Sigma_L(E)-\Sigma_R(E))^{-1} $を計算すると$\sim O(m^3 N^3)$かかってしまう。
$N$は十分大きい数をとりたいが、その場合計算量が膨大になるので、逆行列$G_C$ の$N\times N$個すべてのブロックを計算するのではなく$G_{N1}$の１ブロックだけを計算することで計算量を$\sim O(m^3N)$に落とすことができる。$E - H_C - \Sigma_L - \Sigma_R$がブロック三重対角行列であることを利用して、$G$についての漸化式を作ることができるからである。

左(右)接続グリーン関数

個の漸化式をつくるためには左(右)接続グリーン関数$G^{L(n)},G^{R(n)}$を定義する(リードのグリーン関数と記号が似ているが、違うものなので注意)。これらはそれぞれ左(右)の端子とそれにつながる$n$個のサイトだけのシステムに対するグリーン関数である。つまり、

\begin{align}
G^{L(n)} &= \begin{pmatrix}
E-H_1-\Sigma_L & -V_{12} & 0 & \cdots & 0 \\
-V_{21} & E-H_2 & -V_{23} & \cdots & 0 \\
0 & V_{32} & \ddots & & 0 \\
\vdots & & & \ddots & -V_{n-1,n} \\
0 & \cdots & 0 & -V_{n,n-1} & E - H_n - \delta_{nN} \Sigma_R
\end{pmatrix}^{-1} \quad (1\leq n \leq N)\\

G^{R(n)} &= \begin{pmatrix}
E-H_n -\delta_{n1}\Sigma_L & -V_{n,n+1} & 0 & \cdots & 0 \\
-V_{n+1,n} & E-H_{n+1} & -V_{n+1,n+2} & \cdots & 0 \\
0 & V_{n+2,n+1} & \ddots & & 0 \\
\vdots & & & \ddots & -V_{N-1,N} \\
0 & \cdots & 0 & -V_{N,N-1} & E - H_N - \Sigma_R
\end{pmatrix}^{-1} \quad (1 \leq n \leq N) \\
\end{align}

ちなみに、$G^{L(N)}=G^{R(1)}=G$である。
また、行列のインデックスだが、元の$G$に合わせて$G^{L(n)}$は$1\cdots n$、$G^{R(n)}$は$n\cdots N$をとる。
左(右)接続グリーン関数の対角成分と$n,1$成分について次の漸化式が成り立つ(ブロック対角成分を$G_{nn} = G_n$と省略して表す)。

\begin{align}
&\begin{cases}
G^{L(1)}_1 = G^{L(1)}_{11} = (E-H_1-\Sigma_L)^{-1}\\
G^{L(n+1)}_{n+1} = (G^0_{n+1} - V_{n+1,n}G^{L(n)}_{n}V_{n,n+1})^{-1} \\
G^{L(n+1)}_{n+1,1} = G^{L(n+1)}_{n+1}V_{n+1,n}G^{L(n)}_{n,1}
\end{cases} \\

&\begin{cases}
G^{R(1)}_N = G^{L(1)}_{NN} = (E-H_N-\Sigma_R)^{-1}\\
G^{R(n-1)}_{n-1} = (G^0_{n-1} - V_{n-1,n}G^{R(n)}_{n}V_{n,n-1})^{-1} \\
G^{R(n-1)}_{N,n-1} = G^{R(n)}_{Nn}V_{n,n-1}G^{R(n-1)}_{n-1}
\end{cases}
\end{align}

さらに、元のグリーン関数$G$については漸化式

\begin{cases}
G_n = [ (G^0_n)^{-1} - \delta_{n>1} V_{n,n-1}G^{L(n-1)}_{n-1} V_{n-1,n} - \delta_{n<N} V_{n,n+1}G^{R(n+1)}_{n+1}V_{n+1,n} ] ^{-1} \\
G_{n1} = G_n V_{n,n-1} G^{L(n-1)}_{n-1,1}\\
G_{Nn} = G^{R(n+1)}_{N,n+1} V_{n+1,n} G_n
\end{cases}

よって$E \rightarrow E + i\delta$と置き換えることで $G^{(R)}_{N1}$（この$(R)$は遅延グリーン関数の記号であって、
右接続の記号ではない）を得て、さらに$G^{(A)}_{1N} = (G^{(R)}_{N1})^\dagger$となる。

端子の自己エネルギーの計算

$\Sigma_{L/R}$がどうなるかはここまでで議論していなかったが、実は端子のハミルトニアン$H_{L/R}$に対して再帰グリーン関数法を用いることで計算できる。
端子も$N \gg 1$の同種のサイトが１列につながったものと仮定して、ハミルトニアン

H = \begin{pmatrix}
H_0 & V & 0 & \cdots & 0 \\
V^\dagger & H_0 & V & \cdots & 0 \\
\vdots & V^\dagger & \ddots & & 0 \\
\vdots & & & \ddots & V \\
0 & \cdots & 0 & V^\dagger & H_0
\end{pmatrix}

で表す。そして再帰グリーン関数法を適応した結果、

\begin{align}
G^0 &= (E-H_0)^{-1} \\

&\begin{cases}
t_0 = G^0V^\dagger \\
\tilde{t}_0 = G^0V
\end{cases} \\

&\begin{cases}
t_i = ( 1 - t_{i-1}\tilde{t}_{i-1} - \tilde{t}_{i-1}t_{i-1} )^{-1} t^2_{i-1} \\
\tilde{t}_i = ( 1 - t_{i-1}\tilde{t}_{i-1} - \tilde{t}_{i-1}t_{i-1} )^{-1} \tilde{t}^2_{i-1} \\
\end{cases} \\

&\begin{cases}
T = [\sum_{n=1}^\infty(\prod_{m=0}^{n-1}\tilde{t}_m)]t_n \\
\tilde{T} = [\sum_{n=1}^\infty(\prod_{m=0}^{n-1}t_m)]\tilde{t}_n \\
\end{cases}
\end{align}

とすると、左端のサイトの成分$G_1$と右端の成分$G_N$は

\begin{cases}
G_{11} = (E-H_0 - VT)^{-1} \\
G_{NN} = (E-H_0 - V^\dagger\tilde{T})^{-1} \\
\end{cases}

となる。
すると、システムの右端と右の端子の左側、システムの左端と左の端子の右側が結合しているので、

\Sigma_L = V_{CL} G_{NN} V_{LC} \\
\Sigma_R = V_{CR} G_{11} V_{RC} \\

となる。

実装

端子の自己エネルギーの計算

もっとも簡単なケースとして次の半無限の三重対角ハミルトニアンで表される端子を考える。

H = \begin{pmatrix}
E_0 & t & 0 & \cdots & 0 \\
t & E_0 & t & \cdots & 0 \\
\vdots & t & \ddots & & 0 \\
\vdots & & & \ddots & \ddots \\
0 & \cdots & 0 & \ddots & \ddots
\end{pmatrix}

このときのグリーン関数の端の成分は解析的に求めることができて

\begin{align}
G^{(R)}_{11}(E) &= \begin{cases}
\frac{ \Delta - \sqrt{\Delta^2-1} }{t} & (\Delta \geq -1) \\
\frac{ \Delta - i\sqrt{1 - \Delta^2} }{t} & (|\Delta| \leq 1) \\
\frac{ \Delta + \sqrt{\Delta^2-1} }{t} & (\Delta \geq 1)
\end{cases} \\
\Delta &= \frac{E-E_0}{2t}
\end{align}

となる。
参考資料:

これを十分大きな長さ$N$の端子について再帰グリーン関数法で求めた場合と比べてみよう。
サイトポテンシャル$E_0 = 0$、ホッピング$t=1$とした。
まず、解析解は

lead.py

import numpy as np
from numpy import linalg as LA

def G_lead_end0(E,E0,t):
    D = (E-E0)/(2*t)
    if(D>=1): return complex((D-np.sqrt(D**2-1))/t,0)
    if(abs(D)<1): return complex(D,-np.sqrt(1-D**2))/t
    return complex((D+np.sqrt(D**2-1))/t,0)

で計算される。

また、端子の再帰グリーン関数による計算は

lead.py

def self_energy_in_G_lead(E,H0,V,VL,VR,n_iter,delta=1e-3):
    dagger = lambda M: np.conj(M.T)
    dim = len(H0)
    I = np.identity(dim).astype("complex")
    Zero = np.zeros((dim,dim)).astype("complex")
    G0 = LA.inv( complex(E,delta)*I - H0 )
    A = G0.dot(dagger(V)); B = G0.dot(V)
    # multiplication of A[0]*...*A[n]
    AA = I.copy(); BB = I.copy()
    Asum = Zero.copy(); Bsum = Zero.copy()
    for i in  range(n_iter):
        Asum += BB.dot(A); Bsum += AA.dot(B);  
        AA = AA.dot(A); BB = BB.dot(B)
        W = LA.inv( I - A.dot(B) - B.dot(A) )
        A = W.dot(A.dot(A)); B = np.dot( W, B.dot(B) )
    GL_inv = complex(E,delta)*I - H0 - V.dot(Asum)
    GR_inv = complex(E,delta)*I - H0 - dagger(V).dot(Bsum)
    return dagger(VL).dot(LA.inv(GL_inv)).dot(VL), VR.dot(LA.inv(GR_inv)).dot(dagger(VR))

で計算される。

実行すると、解析解は

lead.py

E_arr = np.linspace(-10,10,2**7)
Gr = np.array([ G_lead_end0(E,0,1).real for E in E_arr ])
Gi = np.array([ G_lead_end0(E,0,1).imag for E in E_arr ])

数値解は端子の長さを$N=2^{10}$、遅延グリーン関数で$i\omega_n\rightarrow \omega + i\delta$としたときの微小量$\delta=1\times 10^{-3}$とすると、

lead.py

H0 = np.array([ [0,0],[0,0] ]).astype("complex")
V = np.array([ [1,0],[0,1] ]).astype("complex")
VL=VR = np.array([ [1,0],[0,1] ]).astype("complex")
G_recursive = np.array([ self_energy_in_G_lead(E,H0,V,VL,VR,2**10,delta=1e-3)[0][0,0] for E in  E_arr ])

これ図示すると

lead.py

from matplotlib import pyplot as plt

plt.plot(E_arr,Gr,label="Re(G)")
plt.plot(E_arr,Gi,label="Im(G)")
plt.plot(E_arr,G_recursive.real,label="Re(G) recursive",marker="x",ms=4,lw=0)
plt.plot(E_arr,G_recursive.imag,label="Im(G) recursive",marker="x",ms=4,lw=0)
plt.legend()

となり。下図のように一致している（厳密には$\delta = 1\times 10^{-3}$の影響の分だけずれている。）。

システムのグリーン関数

各サイトのハミルトニアン

H = \begin{pmatrix}
 H_1 & H_2 & \cdots & H_n
\end{pmatrix}

とサイト間の遷移

V = \begin{pmatrix}
 V_{12} & V_{23} & \cdots & V_{n-1,n}
\end{pmatrix}

と端子の自己エネルギー$\Sigma_{L/R}$を与えて、グリーン関数のうちで左端と右端のサイト間のブロック成分$G_{N1}$を再帰グリーン関数法で求めるクラスを作る

recursive.py

import numpy as np
from numpy import linalg as LA

class RecursiveGreenFunction:
    def __init__(self,N,H_list,V_list,delta=1e-3):
        self.N =N #num. of site
        self.delta = delta # minute parameter of reterted Green func.
        self.H = H_list #len = N
        self.V = V_list #len= N-1, list of V[i-1] = V[i-1][i]
        
    def dagger(self,M):
        return np.conj(M.T)
    
    #B*AB
    def mat_Bd_A_B(self,A,B):
        return self.dagger(B).dot(A).dot(B)
    
    #BAB*
    def mat_B_A_Bd(self,A,B):
        return B.dot(A).dot( self.dagger(B) )
    
    # G0^{-1}[i] = (E+i*delta)I - H[i]
    def G0_inv(self,E,i):
        size = len(self.H[i])
        diag_E = complex(E,self.delta)*np.identity(size).astype("complex")
        return diag_E - self.H[i]

    # calc. left-connected Green func. (i=0,1,...,N-1)
    def L_connected_green_func(self,E,self_EL,self_ER):
        # contain self energy
        GL_list = [ LA.inv( self.G0_inv(E,0) - self_EL) ]
        for i in range(1,self.N):
            GL_inv = self.G0_inv(E,i) - self.mat_Bd_A_B(GL_list[-1],self.V[i-1])
            # contain self energy
            if(i==self.N-1): GL_inv -= self_ER
            GL_list.append( LA.inv( GL_inv ) )
        return GL_list

    # calc. right-connected Green func. (i=N-1,N-2,...,0)
    def R_connected_green_func(self,E,self_EL,self_ER):
        # contain self energy
        GR_list = [ LA.inv( self.G0_inv(E,self.N-1) -self_ER) ]
        for i in range(self.N-2,-1,-1):
            GR_inv = self.G0_inv(E,i) -self.mat_B_A_Bd(GR_list[-1],self.V[i])
            # contain self energy
            if(i==0): GR_inv -= self_EL
            GR_list.append( LA.inv( GR_inv ) )
        return GR_list[::-1] #reorder into 0,1,...,N-1

    # calc. G[i][i]
    def main_green_func(self,i,E,self_EL,self_ER,GL_list,GR_list):
        G_inv = self.G0_inv(E,i)
        # contain self energy
        if(i==0): G_inv -= self_EL
        if(i==self.N-1): G_inv -= self_ER
        # left-connect
        if(i>0): G_inv -= self.mat_Bd_A_B(GL_list[i-1],self.V[i-1]) 
        # right-connect
        if(i<self.N-1): G_inv -= self.mat_B_A_Bd(GR_list[i+1],self.V[i])
        return LA.inv( G_inv )


    def green_edge(self,E,self_EL,self_ER):
        #list of left(right)-connected Green functions
        GL_list = self.L_connected_green_func(E,self_EL,self_ER)
        GR_list = self.R_connected_green_func(E,self_EL,self_ER)
        # calc. G[0][0]
        G = self.main_green_func(0,E,self_EL,self_ER,GL_list,GR_list)
        Gn0 = G.copy()
        #put GL[0][0]
        GL_n0 = GL_list[0].copy()
        for i in range(1,self.N):
            #calc. G[i][i]
            G = self.main_green_func(i,E,self_EL,self_ER,GL_list,GR_list)
            #V[i][i-1]
            Vd = self.dagger(self.V[i-1])
            #G[i][0]= G[i][i]V*[i-1]GL[i-1][0]
            Gn0 = G.dot(Vd).dot(GL_n0)
            #GL[i][0] = GL[i]V*[i-1]GL[i-1][0]
            GL_n0 = GL_list[i].dot(Vd).dot(GL_n0)
        return Gn0

これは直接逆行列を求めた場合と同じなので、システム全体のハミルトニアンから直接グリーン関数を求めた場合と比較する。
簡単のために、全てのサイトが同じハミルトニアン$H_0$で、全ての遷移が同じホッピング$V$で表されるとすると、ハミルトニアンは次の関数で得られる。

recursive.py

def hamiltonian_real(n,H0,V_in,periodic=True):
    h_u = np.diag(np.ones(n-1),k=1) 
    h_v = np.diag(np.ones(n-1),k=-1)
    if(periodic):
        h_u[n-1][0] = 1
        h_v[0][n-1] = 1
    return np.kron(np.identity(n),H0) + np.kron(h_u,V_in) + np.kron(h_v,np.conj(V_in).T)

次のセットアップで計算する。

recursive.py

#energy mesh
E_list = np.linspace(-10,10,64)

#site 
n = 2**9
#states of each site
m = 2

#site Hamiltonian
H0 = np.array([
    [ 0.2, 0.0 ],
    [ 0.0, -0.2 ]
]).astype("complex128")


#hopping between 2 sites
V = np.array([
    [ 1.0, 0.0 ],
    [ 0.0, 1.0 ]
]).astype("complex128")

ちなみに、これはz方向の強磁性体を表している。

これを再帰グリーン関数法で計算すると

recursive.py

#calc. by recursive Green function method
H_list = [ H0 for i in range(n) ]
V_list = [ V for i in range(n-1)]
self_E = np.zeros((m,m)).astype("complex128")
RGF = RecursiveGreenFunction(n,H_list,V_list,self_E,self_E)
G_list = np.array([ RGF.green_edge(E) for E in E_list ])

直接、逆行列を計算すると

recursive.py

#exact inverse-matrix calc.
Hr = hamiltonian_real(n,H0,V,periodic=False)
G0 = lambda E,delta,H: LA.inv( complex(E,delta)*np.identity(len(H)).astype("complex") - H )
G_ref_list = np.array([ G0(E,RGF.delta,Hr)[m*(n-1):m*n,0:m] for E in E_list ])

各要素の誤差の最大値を計算すると

recursive.py

#error
abs(G_ref_list-G_list).max()

1.7991094746387624e-14

$G_{N1}$の各成分の実部をプロットすると、

recursive.py

from matplotlib import pyplot as plt

name = lambda i,j: "G[" + str(i) + "," + str(j) + "]"

for i in range(m):
    for j in range(m):
        plt.plot(E_list,G_ref_list[:,i,j].real,label="Re"+name(i,j))
        plt.plot(E_list,G_list[:,i,j].real,lw=0,marker="x",ms=3)
plt.xlabel("energy")
plt.ylabel("Re[G]")
plt.legend(loc="upper right")
plt.savefig("G_rec_re.png")

虚部をプロットすると

recursive.py

for i in range(m):
    for j in range(m):
        plt.plot(E_list,G_ref_list[:,i,j].imag,label="Im"+name(i,j))
        plt.plot(E_list,G_list[:,i,j].imag,lw=0,marker="x",ms=3)
plt.xlabel("energy")
plt.ylabel("Im[G]")
plt.legend(loc="upper right")
plt.savefig("G_rec_im.png")

である。実線は直接逆行列計算をした結果、xは再帰グリーン関数法で計算した結果である。
以上の図から、再帰グリーン関数法でちゃんと所望の成分が計算できている。

透過率の計算

電子透過率、スピン透過率は散乱行列$S$で単位行列やスピン行列を挟めばよくて

transmittance.py

def S_matrix(Gn1,self_energy_L,self_energy_R):
    dagger = lambda M: np.conj(M.T)
    sq_VgL = scipy.linalg.sqrtm(-2*self_energy_L.imag)
    sq_VgR = scipy.linalg.sqrtm(-2*self_energy_R.imag)
    return sq_VgR.dot( Gn1.dot(sq_VgL) )

def transmittance(X,S):
    dagger = lambda M: np.conj(M.T)
    return np.trace( dagger(S).dot(X).dot(S) ).real

とすればよい。

全てまとめると、以下のようになる。

landauer.py

import numpy as np
from numpy import linalg as LA
import scipy.linalg

def self_energy_lead(E,H0,V,VL,VR,n_iter,delta=1e-3):
    dagger = lambda M: np.conj(M.T)
    dim = len(H0)
    I = np.identity(dim).astype("complex")
    Zero = np.zeros((dim,dim)).astype("complex")
    G0 = LA.inv( complex(E,delta)*I - H0 )
    A = G0.dot(dagger(V)); B = G0.dot(V)
    # multiplication of A[0]*...*A[n]
    AA = I.copy(); BB = I.copy()
    Asum = Zero.copy(); Bsum = Zero.copy()
    for i in  range(n_iter):
        Asum += BB.dot(A); Bsum += AA.dot(B);  
        AA = AA.dot(A); BB = BB.dot(B)
        W = LA.inv( I - A.dot(B) - B.dot(A) )
        A = W.dot(A.dot(A)); B = np.dot( W, B.dot(B) );
    GR = LA.inv( complex(E,delta)*I - H0 - V.dot(Asum) )
    GL = LA.inv( complex(E,delta)*I - H0 - dagger(V).dot(Bsum) )
    return dagger(VL).dot(GL).dot(VL), VR.dot(GR).dot( dagger(VR) )

class RecursiveGreenFunction:
    def __init__(self,N,H_list,V_list,delta=1e-3):
        self.N =N #num. of site
        self.delta = delta # minute parameter of reterted Green func.
        self.H = H_list #len = N
        self.V = V_list #len= N-1, list of V[i-1] = V[i-1][i]
        
    def dagger(self,M):
        return np.conj(M.T)
    
    #B*AB
    def mat_Bd_A_B(self,A,B):
        return self.dagger(B).dot(A).dot(B)
    
    #BAB*
    def mat_B_A_Bd(self,A,B):
        return B.dot(A).dot( self.dagger(B) )
    
    # G0^{-1}[i] = (E+i*delta)I - H[i]
    def G0_inv(self,E,i):
        size = len(self.H[i])
        diag_E = complex(E,self.delta)*np.identity(size).astype("complex")
        return diag_E - self.H[i]

    # calc. left-connected Green func. (i=0,1,...,N-1)
    def L_connected_green_func(self,E,self_EL,self_ER):
        # contain self energy
        GL_list = [ LA.inv( self.G0_inv(E,0) - self_EL) ]
        for i in range(1,self.N):
            GL_inv = self.G0_inv(E,i) - self.mat_Bd_A_B(GL_list[-1],self.V[i-1])
            # contain self energy
            if(i==self.N-1): GL_inv -= self_ER
            GL_list.append( LA.inv( GL_inv ) )
        return GL_list

    # calc. right-connected Green func. (i=N-1,N-2,...,0)
    def R_connected_green_func(self,E,self_EL,self_ER):
        # contain self energy
        GR_list = [ LA.inv( self.G0_inv(E,self.N-1) -self_ER) ]
        for i in range(self.N-2,-1,-1):
            GR_inv = self.G0_inv(E,i) -self.mat_B_A_Bd(GR_list[-1],self.V[i])
            # contain self energy
            if(i==0): GR_inv -= self_EL
            GR_list.append( LA.inv( GR_inv ) )
        return GR_list[::-1] #reorder into 0,1,...,N-1

    # calc. G[i][i]
    def main_green_func(self,i,E,self_EL,self_ER,GL_list,GR_list):
        G_inv = self.G0_inv(E,i)
        # contain self energy
        if(i==0): G_inv -= self_EL
        if(i==self.N-1): G_inv -= self_ER
        # left-connect
        if(i>0): G_inv -= self.mat_Bd_A_B(GL_list[i-1],self.V[i-1]) 
        # right-connect
        if(i<self.N-1): G_inv -= self.mat_B_A_Bd(GR_list[i+1],self.V[i])
        return LA.inv( G_inv )


    def green_edge(self,E,self_EL,self_ER):
        #list of left(right)-connected Green functions
        GL_list = self.L_connected_green_func(E,self_EL,self_ER)
        GR_list = self.R_connected_green_func(E,self_EL,self_ER)
        # calc. G[0][0]
        G = self.main_green_func(0,E,self_EL,self_ER,GL_list,GR_list)
        Gn0 = G.copy()
        #put GL[0][0]
        GL_n0 = GL_list[0].copy()
        for i in range(1,self.N):
            #calc. G[i][i]
            G = self.main_green_func(i,E,self_EL,self_ER,GL_list,GR_list)
            #V[i][i-1]
            Vd = self.dagger(self.V[i-1])
            #G[i][0]= G[i][i]V*[i-1]GL[i-1][0]
            Gn0 = G.dot(Vd).dot(GL_n0)
            #GL[i][0] = GL[i]V*[i-1]GL[i-1][0]
            GL_n0 = GL_list[i].dot(Vd).dot(GL_n0)
        return Gn0

def S_matrix(Gn1,self_energy_L,self_energy_R):
    dagger = lambda M: np.conj(M.T)
    sq_VgL = scipy.linalg.sqrtm(-2*self_energy_L.imag)
    sq_VgR = scipy.linalg.sqrtm(-2*self_energy_R.imag)
    return sq_VgR.dot( Gn1.dot(sq_VgL) )

def transmittance(X,S):
    dagger = lambda M: np.conj(M.T)
    return np.trace( dagger(S).dot(X).dot(S) ).real

delta = 1e-3

#energy mesh
E_arr = np.linspace(-10,10,64)

##### lead #####
n_lead = 2**10
H0_lead = np.array([ [0,0],[0,0] ]).astype("complex")
V_lead = np.array([ [1,0],[0,1] ]).astype("complex")

tL=tR=1
VL = tL*np.array([ [1.,0.],[0.,1.] ]).astype("complex128")
VR = tR*np.array([ [1.,0.],[0.,1.] ]).astype("complex128")

##### system #####
#site 
n = 2**8
#states of each site
m = 2

#site Hamiltonian
H0 = np.array([
    [ 0.2, 0.0 ],
    [ 0.0, -0.2 ]
]).astype("complex128")


#hopping between 2 sites
V = np.array([
    [ 1.0, 0.0 ],
    [ 0.0, 1.0 ]
]).astype("complex128")

H_list = [ H0 for i in range(n) ]
V_list = [ V for i in range(n-1)]

RGF = RecursiveGreenFunction(n,H_list,V_list)

self_E_list = [ self_energy_lead(E,H0_lead,V_lead,VL,VR,n_lead,delta=delta) for E in  E_arr ]
# right-end of system connected with left-end of right lead (visa versa)
Gn0_arr = np.array([ RGF.green_edge(E,SE[1],SE[0]) for E,SE in zip(E_arr,self_E_list) ])
S_matrix_arr = np.array([ S_matrix(Gn0,SE[1],SE[0]) for Gn0, SE in zip(Gn0_arr,self_E_list) ])

up = np.array([
    [1.,0],
    [0.,0.]
])

down = np.array([
    [0.,0],
    [0.,1.]
])

Tup_arr = np.array([ transmittance(up,S) for S in S_matrix_arr ])
Tdown_arr = np.array([ transmittance(down,S) for S in S_matrix_arr ])

from matplotlib import pyplot as plt
plt.plot( E_arr, Tup_arr+Tdown_arr, label="T" )
plt.plot( E_arr, Tup_arr-Tdown_arr, label="Ts" )
plt.xlabel("energy")
plt.ylabel("transmittance")
plt.legend()
plt.savefig("landauer_test.png")

計算結果は以下の通りである。up/downスピンの電子の透過率$T_\uparrow,T_\downarrow$は

であるから、電子の透過率$T＝T_\uparrow+T_\downarrow$とスピン透過率$T_S=T_\uparrow-T_\downarrow$は以下のようになる。

参考文献

D. A. Ryndyk , MRS Bulletin 42, 74–75 (2017); Theory of Quantum Transport at Nanoscale: An Introduction
https://link.springer.com/article/10.1007/s10825-013-0458-7
https://iopscience.iop.org/article/10.1088/0305-4608/15/4/009