程 霄，謝 航

重慶大學(xué)物理學(xué)院，重慶，401331

目錄

I.理論部分 119

A.量子輸運簡介 119

B.量子輸運體系的說明 120

C.計算方法簡介 121

D.平衡態(tài)格林函數(shù) 122

E.非平衡格林函數(shù) 124

F.含時輸運的電流公式[3]125

G.單粒子的密度矩陣運動方程理論[27]126

H.復(fù)數(shù)吸收勢的理論[36,37]129

II.含時量子輸運的應(yīng)用 131

A.非對稱的瞬態(tài)電流[27]131

B.石墨烯納米帶中的瞬態(tài)電流[39]131

C.量子自旋霍爾絕緣體中的瞬態(tài)自旋極化電流[43]132

D.多體關(guān)聯(lián)作用體系的含時輸運[31,44]132

致 謝 134

參考文獻 134

含時量子輸運理論是量子輸運領(lǐng)域的一個比較復(fù)雜的問題。輸運問題研究開放體系，受到周圍環(huán)境中無限多個電子的影響，體系的態(tài)密度不同于孤立系統(tǒng)，呈現(xiàn)為連續(xù)的曲線。在單電子近似下，量子體系是多個單電子態(tài)從最低能到費米面附近能態(tài)的線性疊加而成。每個電子態(tài)遵循薛定諤方程；整個體系各個電子態(tài)構(gòu)成了密度矩陣，其演化遵從相應(yīng)的運動方程。在含時理論方面，由于有電極帶來的影響，格林函數(shù)或波函數(shù)的演化還包含有記憶項，也就是個非馬爾科夫過程。數(shù)學(xué)上，由于非平衡的特性，含時格林函數(shù)是定義在雙支實軸上的Keldysh 格林函數(shù)，并且由于電極影響，其運動方程是個微分積分方程。如何求解這種復(fù)雜的含時演化方程，特別是較大體系時如何快速地計算這種演化過程，是一個挑戰(zhàn)。近年來包括作者在內(nèi)的不少學(xué)者做了許多有益的嘗試。

本文以作者比較熟悉的運動方程方法為主，介紹含時量子輸運的基本理論。其中包括格林函數(shù)基礎(chǔ)，非平衡格林函數(shù)理論，電流計算公式和單粒子表象下的密度矩陣運動方程方法和復(fù)數(shù)吸收勢方法。對于其他一些的含時輸運方法，如級聯(lián)運動方程(HEOM)，波函數(shù)演化等方法，我們也有提及。最后我們介紹含時輸運的一些有趣的推廣和應(yīng)用，如次近鄰近似下的非對稱瞬態(tài)電流，石墨烯納米帶中的瞬態(tài)電流，含自旋和多體效應(yīng)的電子動力學(xué)演化等。

I.理論部分

A.量子輸運簡介

納米器件因為其廣闊的應(yīng)用前景，一直以來都是研究的熱點，對其輸運性質(zhì)（器件對電壓的響應(yīng)）的研究是其中很重要的一個方面。納米器件中的量子輸運不同于宏觀物體中由歐姆定律描述的輸運；也不同于原子量級的系統(tǒng)中由單個元胞的能帶性質(zhì)而體現(xiàn)的輸運，其被稱為介觀輸運。介觀，顧名思義就是介于宏觀和微觀之間。晶體中的載流子- 電子用波函數(shù)描述，電子的波函數(shù)除了有波幅信息外還有含有相位信息，正是由于相位信息，兩束相干的電子可以發(fā)生干涉。晶體中的電子在輸運時會發(fā)生定向運動，在此過程中除了和靜態(tài)晶格的彈性碰撞外，還會因為電子-電子和電子-聲子等相互作用發(fā)生非彈性碰撞，這些非彈性碰撞會改變電子的相位。類似電子平均（彈性碰撞）自由程，我們可以定義電子運動過程中的平均非彈性碰撞自由程，又被稱為相位相干長度lφ（對于Au來說，當(dāng)T=1K時，lφ≈1 μm）[1]。如果（費米面上的）電子的德布羅意波長為λB，那么器件的典型長度L 滿足λB?L ≤lφ就被視作介觀系統(tǒng)[2-4]。

所以，介觀系統(tǒng)的標(biāo)志是其電荷載流子（電子或空穴）的相位相干性在其傳輸過程中保持不變，這種相位相干性產(chǎn)生的干涉效應(yīng)可以反映在許多實驗測量中。相干性的一個例子是AB 效應(yīng)(Aharonov-Bohm effect)[5]，其中環(huán)形幾何結(jié)構(gòu)中兩個不同傳輸路徑的相干疊加導(dǎo)致振蕩的磁阻。另一個例子是普適電導(dǎo)漲落(universal conductance fluctuations)[1,6]，當(dāng)外部參數(shù)（磁場或溫度）改變時，樣品的電導(dǎo)率以e2/h（因此是普適性的）為尺度快速變化, 這也是由于體系包含多個導(dǎo)電的通道，由電子波動的干涉效應(yīng)引起的。

用量子力學(xué)處理量子輸運體系，通常需要將無限大的電極的影響，采用表面格林函數(shù)，即自能的方法，“折疊” 到有限的中心（器件）區(qū)域來計算。靜態(tài)的量子輸運，已經(jīng)有成熟的格林函數(shù)等理論，并且結(jié)合到了密度泛函等第一性原理的模型上，可以對實際體系做詳細的計算。而對于含時的量子輸運，目前還在發(fā)展當(dāng)中，有各種計算方法。由于計算量大的因素，大多數(shù)的理論還集中在較小的體系。這里我們著重介紹一種運動方程方法，有計算較大體系的潛力。

盡管輸運問題是用量子力學(xué)描述的，但是由于有較大的電極，整個輸運體系中的電子流動，可以想象為在兩個高低不同的水庫間連接著一個小小的儲水池，電流輸運對應(yīng)為其中水的流動。如果用有劃分的方法(partition scheme)，相當(dāng)于儲水池有兩個閘門，平時關(guān)閉，各自有不同水位（相當(dāng)于其化學(xué)勢）。而某個時刻突然閘門打開（相當(dāng)于電極-器件間的開關(guān)導(dǎo)通），這時高水位的水庫放水進入水池，水池另一端的水也突然流入低水位水庫，池中水開始流動，同時也引起一陣水波的蕩漾。而另外一種無劃分的方法(partition-free scheme)，則認為兩邊水庫和儲水池始終由兩個柔性管道相連。電極施加電壓時，相當(dāng)于虛擬地“抬高”及“降低”兩邊水庫的底部位置，引起其水位也相應(yīng)改變。這時也產(chǎn)生了高水位水庫通過水池向低水位水庫的放水的現(xiàn)象，當(dāng)中的水池根據(jù)兩邊進出水流大小，水位也發(fā)生漲落。

含時輸運體系有著許多有趣而復(fù)雜的現(xiàn)象。如果把其想象成一道落差巨大的瀑布，則可有“飛流直下三千尺，疑是銀河落九天”之“壯觀”的“電流”美景，如果把其想象成一條安靜的小溪，則可有“泉眼無聲惜細流，樹陰照水愛晴柔”的靜謐之境。

當(dāng)然，把量子態(tài)的電子完全比作水流，也有一定的局限性的。畢竟電子是一種波動的體系，這些波動的“流體”在器件中，會產(chǎn)生干涉，隧穿等現(xiàn)象。因此實際的輸運體系中，有著比前面提到的儲水池的水流更加復(fù)雜的“景象”，有一層層的量子“水流”在器件中流動，同時相干疊加在一起，不斷舞動。在新的納米材料，二維材料等體系中，由于材料本身的新的特性（如拓撲性質(zhì)，自旋極化等），也會帶來很多更加新奇的量子“水流”。這些現(xiàn)象，都是值得進一步研究的。

B.量子輸運體系的說明

我們通常用下圖來表示量子輸運體系[2,7]：

圖1.量子輸運系統(tǒng)的示意圖

我們研究的系統(tǒng)如圖1 所示，系統(tǒng)的總哈密頓量可以寫為：H =HL/R+HT+HC，其中HL/R描述了左/右電極(lead)的哈密頓量，HT是兩個電極和中心區(qū)域的隧穿耦合哈密頓量，HC是中心區(qū)域的哈密頓量。

電極部分一般是金屬，其中的電子通?？梢砸暈榉窍嗷プ饔貌⑶姨幱谄胶鉅顟B(tài)的電子“庫”中，由于電極的部分通常遠遠大于中心區(qū)域的納米體系，導(dǎo)電時通過的電量也遠小于電極中原來的電子數(shù)量，所以我們可以把電極想象為一個很大的水庫。其水位等水況不隨放水過程而改變。 因此電子庫（電極）擁有確定的能量分布和相應(yīng)的溫度，以及確定的其化學(xué)勢（相當(dāng)于水庫的水位）μα(α=L,R)。電極中電子服從費米-狄拉克分布：

左右電極的哈密頓量可以寫為：

中心區(qū)域的哈密頓量的形式取決于被研究體系的物理性質(zhì)和計算準(zhǔn)確度的需要。通常廣泛使用的有第一性原理（密度泛函） 模型和緊束縛模型。這里我們以緊束縛模型為例，其哈密頓量可以寫為：

其中tm,n和?n分布表示跳躍(hopping)能量和坐位(on-site)能量。〈···〉表示對相鄰格點（原子）求和。如果我們考慮電子-電子之間的庫倫相互作用，那么中心區(qū)域的哈密頓量通常需要加上哈伯德(Hubbard)哈密頓量[8]：

其中，U 是Hubbard 系數(shù)，通常由第一性計算擬合得到；是和α 相反的自旋，ni,α是第i 個格點上自旋為α 的粒子數(shù)算符。通常我們會采用平均場近似來計算。如下式，上（下）自旋的座位能取決于下（上）自旋粒子密度的平均值

這是一個自洽方程，通常只能通過迭代方法數(shù)值求解。這個Hubbard 自洽方程也可以從Hartree-Fock 理論中，近似推導(dǎo)出來。此外，在第一性原理模型的計算中，每個原子上的密度也滿足泊松方程，而從而影響該原子的電勢，而原子電勢就是上式中的坐位能?n，其也影響該原子位置的各個軌道的電子分布。因而泊松方程和中心區(qū)域的電子輸運方程也組成一個自洽體系，通常只能通過數(shù)值迭代來得到自洽解[9]。

在格林函數(shù)理論計算中，我們通常把總哈密頓量寫成如下的分塊矩陣形式，即左右電極和中心區(qū)域三部分：

通過引入自能的方式，將原本無限大的“電極+ 中心區(qū)域+ 電極”的體系，轉(zhuǎn)換為中心區(qū)域的有限大的體系，而電極的影響，通過表面格林函數(shù)的方式，反映在了自能函數(shù)里面[2]。

C.計算方法簡介

量子輸運可以分為靜態(tài)輸運和含時輸運，靜態(tài)輸運最有名的公式是郎道爾(Landauer)公式[10,11]

其中，fL(ε)和fR(ε)分別是左右電極的費米-狄拉克分布函數(shù)，T(ε)是通過中心區(qū)域的電子的透射率。從郎道爾公式出發(fā)，考慮零溫度左右兩端接頭施加對稱偏壓（即μL=EF+V/2, μR=EF-V/2），并考慮小偏壓的線性響應(yīng)的情況，我們可以得到電導(dǎo)率的計算公式G=(2e2/h)T，其中T 是費米能級EF處計算得到的電子透射率。

透射率反映了納米器件的幾何和原子特性，可以通過格林函數(shù)公式計算得到[2,7]

上式的Gr(E)和Ga(E)是推遲(retarded)和超前(advanced)格林函數(shù)，ΓL(E)和ΓR(E)分別是左邊和右邊的線寬函數(shù)，其值為自能函數(shù)的虛部乘以(-2)。tr[]表示對矩陣求跡。格林函數(shù)和線寬函數(shù)的定義和計算后面會說明。

另外，體系的局域能態(tài)密度(local density of states，LDOS)分布可以用下面公式計算

上式是針對平衡（無偏壓）情況的，如果有偏壓，應(yīng)用下面的公式計算[7]：

靜態(tài)的量子輸運的透射率計算，常見的方法有格林函數(shù)方法，波函數(shù)匹配法，傳輸（散射）矩陣方法，線性響應(yīng)理論等。具體可以參加相應(yīng)的參考文獻或書籍。

如果電極或中心區(qū)域的哈密頓量隨時間變化，就是含時輸運問題。這時這種系統(tǒng)的計算分析比較復(fù)雜。目前主要有幾類計算方法：基于非平衡格林函數(shù)(Nonequilibrium Green function, NEGF)的方法[14-17]； 基于波函數(shù)的波函數(shù)演化方法[18,19]；基于量子主方程(Quantum Master Equation, QME)的方法[20-23]等。

另外近年來人們對納米體系中的熱（聲子）傳輸也做了很多研究（如新加坡國大的Wang, Jian-Sheng，李保文等）。由于聲子是玻色子，和電子的格林函數(shù)理論（見后）類似，也可以用聲子的產(chǎn)生和湮滅算符定義各種聲子格林函數(shù)，從而建立聲子的NEGF 輸運理論[24,25]。這當(dāng)中的粒子分布函數(shù)也相應(yīng)地換為波色-愛因斯坦分布。另外對于電子-聲子的相互作用，也有不少理論。其中有在自洽波恩近似(self-consistent Born approximation, SCBA)框架下的低階展開等方法[25,26]。

本文中提到NEGF 理論框架下的密度矩陣運動方程(density matrix equation of motion, DM-EOM)方法[27]，是從Tanimura, 嚴以京, 陳冠華，鄭曉等人提出的級聯(lián)運動方程(Hierarchical Equation of Motio,HEOM)理論[15,16,28,29]進一步發(fā)展出來的。鄭曉，嚴以京等人最近也把HEOM 理論應(yīng)用到含多體相互作用的小體系，研究Kondo 共振等現(xiàn)象[30,31]。 不過要了解這種運動方程理論，我們必須先要了解其前置理論：格林函數(shù)技術(shù)。

D.平衡態(tài)格林函數(shù)

數(shù)學(xué)上，格林函數(shù)定義為微分算符L 的逆，即格林函數(shù)G 滿足[32]：

再結(jié)合一定的邊界條件就可以確定格林函數(shù)G 的具體形式。這樣通過格林函數(shù)我們就可以把微分方程轉(zhuǎn)化為積分方程，而積分方程可以用迭代法進行逼近求解?，F(xiàn)在我們考慮含時薛定諤方程：

那么相應(yīng)的格林函數(shù)G(r,r′；t,t′)滿足：

即在離散化的表象下

在能量域下，我們利用傅里葉變換，能夠得到其延遲格林函數(shù)為

上面是離散基矢的寫法，格林函數(shù)表現(xiàn)為一個方矩陣。公式中的η 為一個正的小虛數(shù)，之所以加上這個，是為了保證數(shù)值計算時的收斂?？梢钥闯鯣r(z)在復(fù)平面上的奇點是在上半平面上的，如果η 為負數(shù)，得到的是超前格林函數(shù)， Ga(z)在復(fù)平面上的奇點是在下半平面。

在多體系統(tǒng)中哈密頓量用二次量子化語言描述，延遲格林函數(shù)定義為[3,7]

其中, Ψ(rt)和Ψ?(r′t′)是來自rt 時空坐標(biāo)處的粒子湮滅算符和時空坐標(biāo)處的粒子生成算符；{} 是針對費米子的反對易算符：{A,B}=AB +BA。〈···〉是求平均，θ(t)是海維賽德(Heaviside)階躍函數(shù)，G的上標(biāo)r表示其滿足的一種邊界條件- 延遲邊界條件，即只有當(dāng)t ≥t′此函數(shù)才不為0。為了展示上式定義的格林函數(shù)其實是滿足原始定義式的，我們可以對其兩端應(yīng)用 i?t：

而Ψ(rt)滿足的薛定諤方程為

把這些代入上式整理后得到

類似的，我們還可以定義其它邊界條件下的格林函數(shù)[3]：

其中Ga是超前格林函數(shù)，只在t ≤t′有不為0 的值，G＜是小于(lesser)格林函數(shù)；而G＞是大于(greater)格林函數(shù)；GT是編時(time-order)或因果(causal)格林函數(shù)； T{···} 是編時算符，它總是把擁有更早時間（通常是較小t 值）的算符往右邊移動，即：

其中第二項的負號來自于費米產(chǎn)生/湮滅算符交換的換號特性。故編時格林函數(shù)表示為

所以，我們可以用大于小于格林函數(shù)來表示編時格林函數(shù)

現(xiàn)在我們考慮如何計算這些格林函數(shù)。首先考慮類似(2)式所示的自由粒子哈密頓量，我們可以很容易的計算出結(jié)果。由運動方程（我們用自然單位?=1，e=1）：

如果令初始時刻t0=0，可以得到：

所以，我們可以從定義中直接求出自由空間的格林函數(shù)

其中

是初始時刻t0=0 時的費米狄拉克分布。

更進一步，如果我們考慮總哈密頓量 H=H0+V(t)，其中H0是自由粒子哈密頓量，V(t)含時相互作用哈密頓量，我們可以把V(t)當(dāng)成微擾，通過微擾展開的辦法來計算格林函數(shù)。我們知道，薛定諤匯景是波函數(shù)隨時間演化但算符不隨時間演化，海森堡匯景是算符隨時間演化但波函數(shù)不隨時間演化。現(xiàn)在我們介紹一種相互作用匯景。在這個匯景中，波函數(shù)和算符都隨時間演化，但波函數(shù)隨時間演化只決定于微擾項V(t)，而算符隨時間演化只決定于H0，我們分開是為了方便后面的計算。具體來說，對于算符

其中，上標(biāo)I 代表相互作用匯景。對于波函數(shù)ψI，可以得到如下的方程

其中時間演化算符S(t,t′)為：

假設(shè)我們在相互作用匯景中，編時格林函數(shù)可以寫成如下形式

其中分母表示歸一化性質(zhì)。由于相互作用匯景下含時的場ψI(x1t1)可以用演化算符S 和初始態(tài)的場ψI(x0t0)聯(lián)系起來，當(dāng)在時間t0=-∞時，ψI(x0t0)= ψ0(x0t0)，表示一個非擾動下的態(tài)。任意的物理量A, 其平均值可以寫成

如果系統(tǒng)從初態(tài)ψ0到有擾動的時刻(t≈ 0)，體系的態(tài)變化很小，或擾動很慢（絕熱近似，即用exp(-?|t|)V(t)來代替V(t)，其中? 是無限小的正參數(shù)），這時近似地認為體系變化是可逆的。也就是說，如果在t=0 附近，時間翻轉(zhuǎn)，擾動反過來對稱地從t=0 退回到負無窮的時間處，體系可以演化回到原來的初態(tài)。這種情況就是平衡狀態(tài)。應(yīng)用這種時間對稱性，S(+∞,t)?S(t,-∞)。因此上式(40)可以寫成

上式中T 表示編時算符。在進行格林函數(shù)計算時，可以用類似的計算

這時對應(yīng)的是平衡格林函數(shù)。

E.非平衡格林函數(shù)

在更普遍的情況下，擾動很快或較大時，系統(tǒng)被擾動到達了其他狀態(tài)，即使擾動拆除掉或逆時翻轉(zhuǎn)，系統(tǒng)也恢復(fù)不到原來狀態(tài)。這種情況下，時間對稱性不滿足。A 的平均值就要嚴格按照公式(40)進行計算。時間計算會有一個從負無窮到t(0 附近)， 再從t回到負無窮的過程[3]。

圖2.克爾德什圍道

如果我們定義如上圖所示的一類折疊的路徑，稱為克爾德什(Keldysh)圍道，那么我們可以把上式寫成如下的形式[7]

TC是克爾德什圍道上的編時算符，SC(∞,-∞)≡S-(-∞,∞)S+(∞,-∞)，表示沿著該圍道的傳播子。這樣對于格林函數(shù)G(t,t′)，根據(jù)t 和t′處于不同的分支（上分支C+或下分支C-），可以定義4 個不同的分量：

在用Dyson 公式來推導(dǎo)多個格林函數(shù)乘積的積分時，對非平衡體系，我們是在Keldysh 圍道上進行積分的。但通常，我們是通過真實時間軸（單一實軸）上的積分來計算克爾德什圍道上的積分。這個過程叫做解析延拓。根據(jù)Langreth 定理，對于形如：

的克爾德什圍道積分，解析延拓后得到[3]：

下面我們來推導(dǎo)第二個公式。按照上圖的Keldysh圍道，對于G＜(t,t′)，t 是定義在Keldysh 下支上的，而t′是定義在上支上的。這樣，可以把Keldysh 圍道寫成上下支兩部分，

其中第一項是在上支的積分，第二項是在下支的積分。我們進一步把寫為公式(26)的形式

則上支的積分部分寫為：

則下支的積分部分寫為：

因此，(50)和(51)之和，就得到公式(48)的積分部分。

F.含時輸運的電流公式[3]

我們這里用密度矩陣的劉唯爾運動方程理論來討論含時量子輸運問題。在量子力學(xué)中，劉維爾-凡諾依曼方程描述了密度矩陣的含時演化

在通常的輸運問題中，系統(tǒng)分為電極（含左(L)右(R)兩部分）和中間的器件(D)幾部分。我們用分塊矩陣的形式，把哈密頓量寫為(7)形式，密度矩陣寫成下面形式

相應(yīng)地，上面的運動方程就可以寫成下面的分塊矩陣形式：

其中的耗散項Qα(t)反映了電子通過電極而消耗出（或進入）器件區(qū)域。從上面可以看出，在含電極的開放體系當(dāng)中，獲得跨越器件和電極間的耦合密度矩陣σα,D(t)，是求解含電極影響的器件密度矩陣的演化的關(guān)鍵所在。在非平衡格林函數(shù)里面中，密度矩陣和小于格林函數(shù)有如下的聯(lián)系

這個可以從下面的定義得到

Gkα,m(t,t′)滿足如下的Dyson 方程

其中Gkα,m(t,t′), gka,k1α(t,τ)和Gl,m(τ,t′)分別是按Keldysh 圍道排序的器件- 電極耦合格林函數(shù)，電極格林函數(shù)和器件格林函數(shù)。kα和k1α表示電極α 中的波矢指標(biāo), l 和m 表示器件區(qū)的兩個軌道指標(biāo)。把上式替代σα,D(t)，代入開放的密度矩陣運動方程，可以得到

其中第一項中前面三個量的乘積，可以定義為自能項：

從而上面的耗散項可以寫為[14,15]：

注意上面的時間積分是按Keldysh 圍道C 來進行的。為了化簡為通常實軸上的積分，需要用到前面提到到Langreth 定理的解析延拓關(guān)系。具體到這里的問題，我們可以利用下面關(guān)系：

從而上面的Keldysh 圍道積分可以化簡為下面的表達式

這就是用體系的格林函數(shù)和自能表達的耗散項公式。從這個公式，可以進一步得到各個電極的電流表達式：

注意這個動態(tài)電流的表達式，各個電極可以是不同大小的。這和靜態(tài)電流的公式不一樣。靜態(tài)輸運時，由于體系中電荷是不變的，從電流連續(xù)性條件可以推導(dǎo)出，兩邊電極的電流始終是一樣的。另外，在交流電壓時，可以把式(62)通過傅里葉變換，轉(zhuǎn)換為頻率域，也能得到類似的公式[33]。

G.單粒子的密度矩陣運動方程理論[27]

我們已經(jīng)知道了瞬時電流的表達式，但其中含時大于小于格林函數(shù)和大于小于含時自能的不是那么容易計算，我們想要找到一種簡便的計算方法。前面我們用運動方程技術(shù)求解了格林函數(shù)，可以預(yù)期的是，只要總哈密頓量H 中不含有如(5)式的相互作用項（包含4 個產(chǎn)生/湮滅算符），那么這種技術(shù)產(chǎn)生的方程組就是封閉的，叫做級聯(lián)運動方程。

從上面的耗散項可以看出

其中，我們定義了一階輔助密度矩陣φα(t)

這是我們下一次用運動方程求解的目標(biāo)。如果我們求解到φα(t)，就可以通過(63)式得到Qα(t)，進而得到電流。

φα(t)的運動方程，需要針對(64)式的兩邊對時間求導(dǎo)。在對(64)式右邊時間求導(dǎo)時，求導(dǎo)算符可以作用到三種地方：格林函數(shù)，自能，以及積分號上。對于前兩種，我們需要計算出和的表達式。對于中心區(qū)域無相互作用的系統(tǒng)，從時域的Dyson 方程[7,15]（采用Keldysh 路徑），可以得到：

其中，

其中，f(ε)是費米分布函數(shù)，A(ε)被稱為譜函數(shù)(spectral function)：

其中，Γα(ε)被稱為線寬(level-width)函數(shù)，其和自能的虛部有關(guān)

它的矩陣元為：

現(xiàn)在來看電極小于格林函數(shù)的表達式。對于穩(wěn)態(tài)而言，電極上的施加的電壓不隨時間變化，即（2.2）式上的在位能εkα不隨時間變化，這時的g＜如(69)式所示。對于含時的情況，我們在兩端電極上施加隨時間變化的電壓Δα(t)，在位能也隨時間變化εkα(t)=其中，是穩(wěn)態(tài)的時候的在位能，這時g＜表示為

在這里及以下敘述中，我們用上標(biāo)～表示穩(wěn)態(tài)相關(guān)的量。對于Σ＜也有

從而能量分量（譜形式）的小于（大于）自能可以寫成

利用(65)、(80)式我們可以計算得到：

φαα′(ε,ε,′t)被稱為二階輔助密度矩陣的能量分量，是我們再下一次用運動方程求解的目標(biāo)，其表達式為：

φαα′(ε,ε′,t)的運動方程，可以從定義出發(fā)，用上面同樣的辦法得到。經(jīng)過一系列計算，

最后，把(54)式改寫成能量分量的形式，即：

公式(82)、(84)和(85)就是我們需要的DM-EOM的中心方程，其中的未知數(shù)為 {σC(t), φα(ε,t),φαα′(ε,ε′,t)}。

在實際計算中，上面是包括能量積分的微分方程組，直接求解，計算量太大（比如可以采用Gauss-Legendre quadrature 積分公式，化為近似求和來求解[16]）。我們可以進一步把上述DM-EOM 中的積分化為求和來簡化計算。從以上的推導(dǎo)中我們看到φα(ε,t)定義中的積分號來自于傅里葉變換，我們可以把其中的積分，通過留數(shù)定理改寫為求和。具體來說，從小于自能函數(shù)出發(fā)，

現(xiàn)在要利用留數(shù)定理計算這個積分，我們需要知道fα(z)Γα(z)的極點。其中，fα(z)的極點用Pad′e展開來得到（極點在虛軸上，見圖3 黑色點）[34]：

圖3.DM-EOM 方法中小于自能函數(shù)的圍道積分和極點分布

而Γα(z)可以通過用洛倫茲函數(shù)擬合來得到（極點見圖3 紅色點）：

這里用了Nd個洛倫茲函數(shù)來對線寬函數(shù)進行最小二乘法的擬合。其中Ωd，Wd和ηd分別是第d 個洛倫茲函數(shù)的中心，寬度和幅度，是針對線寬矩陣各個矩陣元的調(diào)節(jié)矩陣。由于線寬函數(shù)每個非零矩陣元都是能量的函數(shù)，如果單獨對每個獨立的非零矩陣元進行洛倫茲擬合，再把所有擬合的洛倫茲函數(shù)合并起來，那么Nd非常巨大，這樣導(dǎo)致對稍微大一點的體系，計算量很大。所以，我們通常選取一組比較合適的公共的洛倫茲函數(shù)集合，把每個非零矩陣元隨能量變化的曲線合并在一起作為目標(biāo)函數(shù)，再調(diào)節(jié)這組公共的洛倫茲函數(shù)組的所有參數(shù)，進行全局性的擬合[35]。這樣可以大大減少洛倫茲函數(shù)的個數(shù)。而且實際計算時，不同的初始擬合條件，可能會得到不同的擬合結(jié)果。但是只要擬合得足夠好，最終計算結(jié)果都是一致和收斂的。另外，除了純數(shù)值的洛倫茲擬合方法外，還可以用另外一種復(fù)數(shù)吸收勢的方法對線寬函數(shù)進行洛倫茲展開。這個方法見后面介紹。

對fα(z)和Γα(z)進行了展開和擬合后，可以看出它們都是有奇點的解析函數(shù)。因此，我們可以用留數(shù)定理，把實軸上的積分化為某個圍道積分，從而用留數(shù)和代替該積分。復(fù)變函數(shù)積分的圍道如圖3 所示，除了實軸部分，還有上下平面兩個半圓。之所以取兩種圍道，是因為這樣的取法可以保證當(dāng)圍道半徑R →∞時，沿著半圓路徑的積分收斂為0。從在一階輔助密度矩陣和二階輔助密度矩陣定義中，兩個時間變量分別有滿足τ - t ＜0和t2-t ＜0 的條件。而為了保證R →∞時，在虛軸上，ε=iR exp[-iε(τ -t)]= exp[-R|τ -t|]收斂，不發(fā)散，我們需要對兩種情況規(guī)定圍道在上半復(fù)平面或下半復(fù)平面。最后，利用留數(shù)定理，我們得到

由上式可以看出，小于自能函數(shù)（穩(wěn)態(tài)）可以展開為Nk項。Nk越大，線寬函數(shù)越精確，但是計算成本也越大。進一步，利用(89)式，一階和二階輔助密度也可以由能量積分轉(zhuǎn)換為離散項求和形式：

最后，我們可以得到如下的離散版本的DM-EOM 方程：

對于上述的微分方程組，通?？梢杂盟碾A的龍格-庫塔(Runge-Kutta)方法求解。

圖4.復(fù)試吸收勢的示意圖；右：輸運計算時的左電極-器件-右電極分區(qū)的哈密頓矩陣，中間部分是CAP 器件的較小的電極-器件耦合區(qū)域。

H.復(fù)數(shù)吸收勢的理論[36,37]

在傳統(tǒng)的NEGF 理論框架中，一個方向上為半無窮大的電極的計算，是用迭代的格林函數(shù)(recursive Green's function)方法來進行的。這種方法需要對電極單元的哈密頓量進行多次迭代，對于較大的第一性原理的體系，有相當(dāng)?shù)挠嬎懔?。在含時輸運方面，對無窮大的電極的線寬函數(shù)的處理，也是很耗時的。對于我們前面提到的DM-EOM 理論，需要用多個洛倫茲函數(shù)的矩陣來擬合電極的線寬函數(shù)矩陣。這種擬合算法所需要的初始參數(shù)，通常是建立在多次擬合的經(jīng)驗上的，而且并沒有唯一性。而不同的擬合結(jié)果，其中包含的微小誤差積累起來，嚴格來說，可能對較長時間的量子演化產(chǎn)生不同的，甚至錯誤的影響。

另一方面，人們從電磁場的含時計算-時域有限差分(DFTD)算法的吸收邊界條件（如完美匹配層，PML）得到啟發(fā)，嘗試也發(fā)明一種能夠很好地吸收入射電子波的邊界區(qū)域。如果這種邊界區(qū)域?qū)Ω鱾€能量電子波都很較好地吸收，那么就能起到模擬無窮大電極的作用（波進入無窮大區(qū)域，不會有任何反射）。在早期的David E.Manolopoulos 的工作中[38]，其提出在有限大小區(qū)域內(nèi)，加上純虛數(shù)的勢場，即所謂的復(fù)吸收勢(complex absorbing potential,CAP)。最優(yōu)化這種勢場，可以使得其對相當(dāng)大范圍內(nèi)的電子波的反射率非常小，起到近似完全吸收。

用有限大小的區(qū)域來代替真正無窮大的電極區(qū)，其哈密頓量的“坐位能” 需要滿足如下CAP 公式：

上面的x1和x2是CAP 區(qū)域的開始和結(jié)束位置，Δx=x2-x1是CAP 區(qū)域的長度；c 是一個常數(shù)，其對最終吸收結(jié)果不太敏感，我們這里取c=1.0。當(dāng)CAP 方法用于較復(fù)雜的原子系統(tǒng)時，我們采用一系列的x 方向的重復(fù)單元來構(gòu)成CAP 區(qū)域，如圖4（左）所示（上部是傳統(tǒng)的兩邊無限大的電極，下部是用有限大的復(fù)吸收區(qū)域代替電極區(qū)）。在靜態(tài)輸運計算時，其格林函數(shù)可以寫為[37]

這里的ψk(r)和Φk(r)是來自孤立的器件體系的本征態(tài)。由于考慮了純虛數(shù)的復(fù)數(shù)吸收勢，器件部分總體的哈密頓量是非厄米的。所以這種非厄米體系的格林函數(shù)計算需要用到兩類本征態(tài)：

其中H0是沒有CAP 的哈密頓量。上面兩式的本征值是共軛的，但是這兩式的本征態(tài)之間并非滿足共軛關(guān)系，也不是正交關(guān)系；并且ψk(r)和Φk(r)自身也沒有歸一性。它們之間滿足所謂的‘雙正交歸一' 關(guān)系。在離散的原子軌道基{φm(r)}表象下，上式可以寫為：

圖5.(a)考慮次近鄰近似下的一維原子鏈體系的幾個線寬函數(shù)的準(zhǔn)確值（用格林函數(shù)迭代計算）；(b)用洛倫茲函數(shù)擬合的該線寬函數(shù)的結(jié)果；(c)對N=4 的該體系，t ＞0 時施加對稱偏壓后的瞬態(tài)電流

注意這里的本征值是個復(fù)數(shù)。從這個復(fù)數(shù)展開式，可以很自然地得到E 在實數(shù)時的洛倫茲函數(shù)。為了表示這個，我們把上式的分子和εk寫為實數(shù)形式：以及εk=Ωk+iWk，那么上式可以改寫為

從這里看出，格林函數(shù)已經(jīng)自然地展開為洛倫茲函數(shù)級數(shù)和的形式了，其中實部和虛部對應(yīng)不同的展開系數(shù)。這里的洛倫茲函數(shù)和普通的形式有些不同，相應(yīng)地，其留數(shù)定理展開時的公式也要做相應(yīng)地修改。具體公式可以參考文獻[36]。

在實際DM-EOM 計算中，洛倫茲展開來自線寬函數(shù)（自能），也就是和半無窮電極的表面格林函數(shù)相聯(lián)系的。我們可以只上面計算格林函數(shù)的通式中只取一個電極的復(fù)吸收勢矩陣，從而得到等效于半無窮的體系的格林函數(shù)。并且，在較大體系，特別是第一性原理的計算時，由于CAP 的區(qū)域較大，展開式中包括很多的極點。在算出整個（以左電極為例）后，我們可以只取中的靠近器件區(qū)域的部分來進一步計算，其對應(yīng)的格林函數(shù)和耦合矩陣分別是和HD,Ls(HLs,D)，見上面圖4（右）。由于通常電極和器件區(qū)的耦合只限于靠近電極器件邊界的區(qū)域，所以這種化簡是有效的，并且因為洛倫茲展開項減少了，可以很大地節(jié)省計算時間。這時自能表示為

上式寫為下標(biāo)的形式，就是

圖6.石墨烯納米帶的電極-器件的劃分圖（左）以及瞬態(tài)電流圖（右）。電流的單位是V σ0，其中，σ0 =e2/h 是量子電導(dǎo)單位，時間的單位是tv =?/v。圖中的兩條虛線表示最小電導(dǎo)率和直流電導(dǎo)率。

其中的分子是

II.含時量子輸運的應(yīng)用

在這部分，我們主要例舉幾個只有在含時量子輸運情況下才能出現(xiàn)的物理現(xiàn)象。這幾個例子也顯示了含時量子輸運的重要性和獨特的魅力。同時，它們各自采用了幾種不同的含時量子輸運的方法，也說明了這一領(lǐng)域目前也還在不斷發(fā)展中。

A.非對稱的瞬態(tài)電流[27]

在多量子點體系中，通常情況下，采用最近鄰的緊束縛近似的一維原子線模型，其能量密度分布關(guān)于E=0 點是對稱的。這時，如果用這種原子線做電極和器件，兩邊對稱地加上電壓，其瞬態(tài)電流每時每刻都是對稱的，直到達到穩(wěn)態(tài)。而在DFT 等模型中，我們通常需要考慮次近鄰等躍遷。比如我們采用ti,i±1=-2.0 eV, ti,i±2=-0.4 eV, 這時線寬函數(shù)就不是對稱分布了，見圖5(a)。我們利用本文提到的DM-EOM 方法來計算含時電流。我們采用22 個洛倫茲函數(shù)來擬合這個線寬函數(shù)矩陣，得到的擬合曲線如圖5(b)。進一步采用上面提到的DM-EOM 方法計算對稱偏壓(±1.0 V)時的左右瞬時電流，見圖5(c)。仔細分析該圖，我們發(fā)現(xiàn)，在開始時刻左右電極的電流并不是對稱的（比較上面的左電極電流和反號的右電極電流，并不重合）。

雖然電壓是對稱施加的，但是由于電極的能態(tài)分布曲線，相對E=0 是上下不對稱的。加偏壓時，左電極的占據(jù)態(tài)電子瞬間進入器件，而器件中的電子瞬間也進入右電極的非占據(jù)態(tài)。這是這種能量分布上的不對稱，導(dǎo)致了開始的瞬時進入和流出的電子數(shù)量不相同，從而左右瞬間的電流不對稱。當(dāng)然，由于器件不可能一直積累電荷，最終兩邊的穩(wěn)態(tài)電流是相同的。

這種非對稱的瞬態(tài)電流現(xiàn)象，在很多實際體系中都是存在的。而通常只能用這種含時的量子輸運理論，才能得到詳細的研究。

B.石墨烯納米帶中的瞬態(tài)電流[39]

Enrico Perfetto，Gianluca Stefanucci 等人采用了一種“長電極” 的策略來進行石墨烯的含時量子輸運計算。這種策略和前面提到的格林函數(shù)方法等“嵌入”(embedding)式處理完全不同。有點類似波函數(shù)演化方法。但是這里并沒有應(yīng)用吸收（透明）邊界條件來處理器件邊緣處的波函數(shù)。作者根據(jù)大部分電子波的傳播是在費米速度vF)以下的特點，采用足夠長的左右電極來模擬無限長的電極，只要計算時間小于Tmax=2Lr/vF，有限大的電極中流出的電子的狀況和無限大的電極的情況是類似的。所以可以用簡單的有限的體系的含時演化來模擬開放體系的量子演化。這種對輸運問題的處理倒類似真實實驗的情況。

考慮寬度為250 個原胞(W≈106 nm)的石墨烯納米帶，把它分為3 個部分：左右電極(Nc=2000)和長度為100 個原胞(Na=100，L≈25 nm)的中心區(qū)域。在左右電極上施加V =8×10-4V 的電壓，中心區(qū)域上的電場為E=(VL-VR)/L，可以計算得到如圖6 所示的瞬態(tài)電流。當(dāng)電壓開始施加的時候，電流迅速上升，然后振蕩著趨近一個暫時穩(wěn)態(tài)電流Imin，此電流是石墨烯在狄拉克點附近具有相對論性質(zhì)的能譜這一事實的體現(xiàn)，與理論計算的石墨烯納米帶的最小電導(dǎo)率精確一致[40,41]：

這個過程的持續(xù)時間為τc=L/vF，其角頻率為ω1=（其中，v =2.7 eV 是石墨烯的跳躍積分），此頻率與石墨烯光學(xué)響應(yīng)中的共振頻率一致，是石墨烯的塊(bulk)性質(zhì)[42]。當(dāng)電子進入電極后，因為電子與沒有電場的電極的相互作用，電流最終會回落到的穩(wěn)態(tài)值IDC（與通過郎道爾計算的結(jié)果一致），這個過程也是振幅衰減的振蕩過程，其頻率ω2=2πvF/W（其中vF是電子的費米速度），而?ω2正是橫向能量子帶的能量間隔。

C.量子自旋霍爾絕緣體中的瞬態(tài)自旋極化電流[43]

香港大學(xué)王健教授組近來用含時輸運理論研究了一些拓撲材料的一些含時自旋極化電流。他們采用了量子自旋霍爾絕緣體- 鍺烯納米帶作為例子。其緊束縛模型如下：

其中，第一項是近鄰跳躍項，h0=1.3 eV 是跳躍積分；第二項描述強度為λSO=0.1 eV 的固有軌道自旋耦合，只含次近鄰〈〈i,j〉〉跳躍項，其中，vij=±1，σz是泡利矩陣的z 分量；第三項描述垂直平面的電場的作用。其中，l=0.4是屈曲(buckling)高度，ξi=±1，i=AB。在此電場的作用下，處于z 軸上不同位置的AB 兩套格點分別具有不同的在位能，這樣電場就破壞了原來的AB 格點的對稱性，系統(tǒng)出現(xiàn)能隙，系統(tǒng)從量子自旋霍爾絕緣相轉(zhuǎn)變?yōu)槠胀ǖ模芟叮┙^緣相。

圖7.鍺烯納米帶上的瞬態(tài)自旋電流Iσ(t)(σ =↑,↓)。(a)Ez =0，(b)lEz =0.1h(Ez ＞0)，(c)Ez =0.1h(Ez ＜0)。內(nèi)圖：瞬態(tài)自旋極化電流P(t)=I↑(t)-I↓(t)。虛線表示通過郎道爾公式計算得到的穩(wěn)態(tài)值。

他們把寬度W 為12 個原胞的鍺烯納米帶分成左右電極和長度L 為16 個原胞的中心區(qū)域3 個部分，左電極上施加的0.2 V 的瞬態(tài)電壓，右電極上的電壓為0 V。不同Ez時的瞬態(tài)自旋電流Iσ(t)（σ =↑,↓）和瞬態(tài)自旋極化電流P(t)=I↑(t)-I↓(t)由含復(fù)吸收勢的NEGF-DFT 方法進行計算[37]，結(jié)果如圖7 所示?？偟膩碚f，瞬態(tài)電流可分為3 個階段：開始階段的強烈振蕩階段，緊接著的穩(wěn)步上升階段和最后的穩(wěn)態(tài)階段。在圖7(a)中，Ez=0，中心區(qū)域處于量子自旋霍爾絕緣相，瞬態(tài)的電流對于自旋來說是完全簡并的。在圖7(b)(c)中，Ez/= 0，中心區(qū)域處于普通（能隙）絕緣相，這時因為上下自旋的電子的散射矩陣的相位不同，導(dǎo)致了其形成的上下自旋電流達到穩(wěn)態(tài)電流所需的特征時間τσ（σ=↑,↓ 不同，具體來說，在圖7(b)中，Ez＞0，τ↑＞τ↓，在圖7(c)中，Ez＜0，τ↑＜τ↓，這種時間延遲產(chǎn)生了瞬態(tài)自旋極化電流P(t)=I↑(t)-I↓(t)，如圖7 內(nèi)圖所示。

D.多體關(guān)聯(lián)作用體系的含時輸運[31,44]

在多體關(guān)聯(lián)體系中，哈伯爾(Hubbard)模型和安德森(Anderson)）雜質(zhì)模型經(jīng)常用到。對單量子點而言，其哈密頓量是

圖8.含相互作用的量子點體系在瞬態(tài)直流偏壓下的情況。左：量子點接入開放體系以及其能級分布圖。右：在不同溫度（見圖中各種顏色曲線）時，器件處于EOR、MVR 和KR 狀態(tài)時瞬時電流[(a),(c),(e)]和譜函數(shù)分布[(b),(d),(f)]。

其中上式就是安德森雜質(zhì)模型，表示一個量子點中兩種自旋的電子所具有的坐位能（勢能）和相互排斥能。第二項也就是哈伯爾項（見公式5）。前面的平均場近似（公式6）通常適用于弱相互作用的情況。

嚴格處理這類多體含時問題，理論和計算上都有一定困難。目前有不少理論可以用到，如含時的NEGF+GW 方法，含時密度矩陣重整化群(density matrix renormalization group, DMRG)等等。這里我們舉的兩個例子采用的是級聯(lián)運動方程(HEOM)方法[29]。該方法是從影響泛函理論出發(fā)推導(dǎo)的，理論上是一組無窮階的微分方程組。當(dāng)處理非相互作用粒子體系時，該方程組自動在二階上閉合，等價于本文中提到的DM-EOM 理論。由于高階HEOM 的計算量非常巨大，目前大部分的多體計算都集中在單量子點體系，并且通常在4～5 階進行截斷。

圖8 是魏建華教授等人最近的工作[44]。他們考察了含相互作用的量子點(quantum dot, QD)- 電極體系加上臺階電壓VSD時的情況。其中左邊的(a)圖中QD 器件連接兩個電極，門電壓(gate volatge)可調(diào)節(jié)量子點的坐位能。左邊的(b)圖描述了電極，量子點的能量狀況。當(dāng)門電壓控制?g降低到-1.5 meV時，電極中電子進入QD 較多，系統(tǒng)處于近藤效應(yīng)區(qū)域(Kondo regime, KR)，此時電子的庫倫相互作用很強；而當(dāng)?g處于0 meV 附近時，QD 部分填充，系統(tǒng)處于混合價帶區(qū)域(mixed valence regime, MVR)；當(dāng)?g處于較高的能量0.6 meV 處時，系統(tǒng)處于空軌道區(qū)域(empty orbital regime, EOR)。圖8 右邊幾幅圖描述了這三組區(qū)域時，不同溫度的瞬態(tài)電流的變化和相應(yīng)的譜函數(shù)圖?？梢钥闯觯诮賲^(qū)時，只有當(dāng)溫度kBT 小于0.07 meV（大致的Kondo 溫度）時，QD 中有近藤振蕩，這種多體振蕩增強了左右電極電流通過QD 的相干隧穿，導(dǎo)致瞬態(tài)電流產(chǎn)生很多振蕩。對應(yīng)地，能譜上也出現(xiàn)了個尖銳的Kondo 峰。而溫度高于Kondo 溫度時，這種共振導(dǎo)致的電流振蕩消失了，只有單粒子體系中的“過沖”電流導(dǎo)致的非線性響應(yīng)。同時，仔細觀察近藤區(qū)的不同溫度的穩(wěn)態(tài)電流，可以看出低溫(kBT =0.03 meV)時的電流反而更小，這和在MVR 和EOR 情況下是相反的（高溫時穩(wěn)態(tài)電流更小）。這也正是近藤效應(yīng)的體現(xiàn)：非常低溫時相互作用系統(tǒng)的電阻，會由于存在多體的近藤效應(yīng)反而會升高。

圖9 中，鄭曉，嚴以京和M.Di Ventra 研究了近藤量子點在交流電極環(huán)境下，由強關(guān)聯(lián)引起的記憶效應(yīng)[31]。 兩個電極隨時變化為：VL(t)=-VR(t)=V0sin(Ωt)，電極的線寬函數(shù)是洛倫茲形式，其耦合常數(shù)是Δ。當(dāng)QD 的坐位能?g處于-2Δ、-6Δ 和-10Δ 時，對應(yīng)的Kondo 溫度分別是0.44Δ、0.025Δ和0.0013Δ。 從上圖可以看出當(dāng)溫度較高時（如10Δ），電流-電壓呈線性關(guān)系，動態(tài)電流也隨時間呈正弦變化；而降到較低溫度時，電流-電壓曲線顯示出較復(fù)雜的有兩次交叉的滯后回線。這說明電流與電壓變化并不是同步的，電流對過去一段時間的電壓有記憶效應(yīng)。實際上，當(dāng)交流頻率Ω 很高時，電流-電壓由于“量子電感”效應(yīng)[45]而呈現(xiàn)出有固定延遲相位的橢圓形滯后回線，Ω 非常低時，滯后曲線會消失。而當(dāng)Ω 在圖9 中的中間值時，Kondo 多體效應(yīng)也會導(dǎo)致電流滯后（或記憶）效應(yīng)。同時，他們也發(fā)現(xiàn)更低的?g和更大的偏壓V0也會壓制近藤效應(yīng)。因為本文篇幅所限，具體細節(jié)可以參考相關(guān)的文獻。

圖9.含相互作用的量子點近藤體系在交流偏壓下的情況。(a)動態(tài)的電流-電壓曲線；小圖是T=0.1 Δ 時曲線的局部放大圖；(b)實時的電流曲線。這里不同顏色代表不同的溫度。其他參數(shù)為：W=20 Δ，eV0=Δ， Ω=0.3 Δ，其中W 是洛倫茲函數(shù)的寬度，V0 是交流電壓幅度，Ω交流電的頻率對應(yīng)的能量。

致 謝

本文作者感謝香港科技大學(xué)的沈平教授在介觀電子輸運方面的悉心指導(dǎo)；感謝香港大學(xué)的陳冠華教授，王健教授，加拿大麥吉爾大學(xué)的郭鴻教授等在含時量子輸運方面的指導(dǎo)和幫助。作者還特別感謝中國科學(xué)技術(shù)大學(xué)的嚴以京教授和鄭曉教授，上海電力大學(xué)的蔣鋒副教授，以及北京計算科學(xué)研究中心的任志勇研究員，美國洛斯阿拉莫斯國家實驗室的張余研究員等人在級聯(lián)運動方程，格林函數(shù)方面的幫助。最后，作者感謝重慶大學(xué)“百人計劃” 啟動基金和國家自然科學(xué)基金理論物理平臺(No.11847301)的支持。