關(guān)夢雪1) 廉超1) 孟勝1)2)?

1)(中國科學(xué)院物理研究所,北京凝聚態(tài)物理國家研究中心,北京 100190)

2)(量子物質(zhì)科學(xué)協(xié)同創(chuàng)新中心,北京 100190)

(2018年3月19日收到;2018年4月17日收到修改稿)

1 引 言

實(shí)時(shí)演化的含時(shí)密度泛函理論[1,2](realtime time dependent density functional theory,rt-TDDFT)中,電子密度由時(shí)間依賴的Kohn-Sham(time-dependent Kohn-Sham,TDKS)方程數(shù)值積分得到,可以直接在非微擾的情況下提供電子波函數(shù)的時(shí)域演化動(dòng)力學(xué),結(jié)合原子核的運(yùn)動(dòng),能夠直觀地給出電子-原子核多體量子體系隨時(shí)間的演化路徑.因此,rt-TDDFT已成為計(jì)算模擬強(qiáng)場物理、超快物理過程的最有效的工具之一[3?8].

目前,rt-TDDFT的應(yīng)用范圍已經(jīng)從原子分子體系拓展到凝聚態(tài)物質(zhì)中,其數(shù)值實(shí)現(xiàn)大多基于平面波[9]或者實(shí)空間網(wǎng)格點(diǎn)方法[10,11].然而,在以上的兩種方法中,如果要研究能量尺度在10—100 eV的高能電子激發(fā),則需要包含很高能量的平面波基矢或者很密集的網(wǎng)格點(diǎn)才能準(zhǔn)確地描述i原子內(nèi)層電子的運(yùn)動(dòng).同時(shí),對于一個(gè)有Na個(gè)原子的體系,計(jì)算所需的基矢數(shù)目高達(dá)103Na—104Na,因此無法進(jìn)行有效的大規(guī)模模擬.這些因素都將顯著提高計(jì)算成本,并大大限制了rt-TDDFT理論的應(yīng)用范圍.

鑒于上述現(xiàn)狀,我們發(fā)展了一套基于數(shù)值原子軌道基的rt-TDDFT計(jì)算方法以及軟件TDAP(Time Dependent Ab initio Package)[12,13].相較于平面波和實(shí)空間展開,用基于數(shù)值化局域原子軌道的基矢展開波函數(shù),要求較少的基數(shù)目就可以很好地描述系統(tǒng)的電子態(tài)信息,大幅減小了計(jì)算量.另外,由于局域軌道只在有限空間分布,在此空間之外嚴(yán)格為零,所以對于大體系只有該空間內(nèi)的原子會(huì)有波函數(shù)的重疊,哈密頓量為準(zhǔn)線性,使得我們能夠獲得計(jì)算量隨尺寸成線性增長的好處.利用TDAP,能夠長時(shí)間地模擬大尺寸系統(tǒng)的動(dòng)力學(xué)演化過程,其計(jì)算成本較于平面波和實(shí)空間網(wǎng)格點(diǎn)方法大幅降低,并且能夠保持很高的精確性.與文獻(xiàn)中的常規(guī)方法[9?11]相比,其優(yōu)勢具體體現(xiàn)在以下四方面.

1)數(shù)值原子軌道基的使用極大地減小了計(jì)算量(計(jì)算量級為10Na),可以快速、準(zhǔn)確地模擬周期性體系或者包含很大真空層的超胞.

2)在每一電子和離子步中,電子密度都將自洽演化并通過平均場理論計(jì)算出激發(fā)態(tài)的實(shí)時(shí)軌跡,可以給出光激勵(lì)下的電子和原子核超快動(dòng)力學(xué)的微觀圖像.

3)原子核附近的電子狀態(tài)可以通過原子軌道的線性組合得到有效的處理,因此可以實(shí)現(xiàn)對內(nèi)層電子激發(fā)的理論模擬,研究較高能量范圍飛秒和亞飛秒時(shí)間尺度的電子動(dòng)力學(xué).

4)受限于計(jì)算效率,傳統(tǒng)的rt-TDDFT大多只能在動(dòng)量空間的單點(diǎn)上進(jìn)行計(jì)算,無法實(shí)現(xiàn)體系的動(dòng)量分辨.近期我們發(fā)展了動(dòng)量分辨的rt-TDDFT算法,可以用很小的單胞計(jì)算固體和表面,從而使計(jì)算成本降低了幾個(gè)數(shù)量級.此外,我們還發(fā)現(xiàn)體系激發(fā)狀態(tài)在動(dòng)量空間上的不同區(qū)域存在明顯差異,表明該方法有利于揭示深入細(xì)致的微觀機(jī)制.

因此,該方法能夠在廣泛的量子系統(tǒng)中精確、有效地處理各類超快動(dòng)力學(xué)過程,成為解釋并預(yù)測新奇量子現(xiàn)象的強(qiáng)大工具.本文簡單介紹TDAP軟件的核心算法,并結(jié)合一些實(shí)例,介紹軟件的主要特色和部分已經(jīng)實(shí)現(xiàn)的功能.

2 計(jì)算方法

rt-TDDFT算法的主要框架繼承于較早版本的TDAP[13],基于開源軟件包SIESTA[14,15].圖1描述了在給定離子步時(shí)對激發(fā)態(tài)模擬的流程圖.與早期版本略有不同的是,本文方法能夠進(jìn)行動(dòng)量分辨的時(shí)間演化,從而能夠適用于有限尺寸低維體系和周期性體系.下面均以這種一般性的動(dòng)量分辨rt-TDDFT算法為例.

2.1 哈密頓量及交迭矩陣

在周期性邊界條件下,晶格矢量用Rs(s=1,2,3,···)表示,單胞中的原子i的位置為bi. 每個(gè)原子都有一套數(shù)值原子軌道{ξiα}與之對應(yīng),其中原子態(tài)的軌道和角量子數(shù)均用α表示,用多重徑向基函數(shù)ζ展開[14].為了簡化方程形式,本文使用原子單位~=me=e=1.同時(shí),除非特別注明,只有顯式依賴于時(shí)間t的算符才以f(t)的形式標(biāo)記.

圖1 動(dòng)量分辨的TDDFT算法流程圖Fig.1.Flowchart of k-resolved TDDFT algorithm.

在每個(gè)k點(diǎn)的交迭矩陣Sk及哈密頓量Hk用原子軌道基表示為

哈密頓算符為

其中為動(dòng)能算符;及分別是原子I贗勢的局域及非局域部分,VH,VXC及Vext分別是Hartree、交換關(guān)聯(lián)及外部電場的勢能.目前,我們主要采用絕熱局域密度近似或者絕熱廣義梯度近似作為交換關(guān)聯(lián)泛函.關(guān)于的詳細(xì)計(jì)算方法已有文獻(xiàn)提及[15],在此不再贅述.

為了模擬光與各種材料之間的相互作用,時(shí)間依賴的電場E(t)以兩種形式體現(xiàn)在哈密頓量中:長度或者速度規(guī)范.在長度規(guī)范下,電場E(t)由標(biāo)量勢Vext引入:

為了防止哈密頓量的平移對稱性被電場E(t)破壞,可以在空間方向μ∈x,y,z上施加一個(gè)鋸齒狀的電場:

其中Lμ是沿著μ方向的原胞長度,并且ε→0.因此需要保證ρ(xμ)在?ε

為了彌補(bǔ)上述不足,通過度規(guī)變換

可以將TDKS方程用速度規(guī)范表達(dá),用于對無限大周期性體系的動(dòng)力學(xué)描述.此時(shí)哈密頓量為

此外,E(t)的波形在這兩種情形下均可任意調(diào)節(jié).例如,可采用高斯波包形式:

其中f為激光頻率,t0為峰值時(shí)間,?為附加相位.

2.2 電子態(tài)演化

有了隨時(shí)間演化的哈密頓量以及交迭矩陣后,利用TDKS方程就可以由上一個(gè)時(shí)間步的狀態(tài)數(shù)據(jù)得到新時(shí)刻的波函數(shù)unk(r,t):

其中unk(r,t)=ψnk(r,t)exp(?ik ·r)是布洛赫波函數(shù)ψnk(r,t)的空間周期性部分,而t′≈(t1+t2)/2,?t=t2?t1為時(shí)間的步長.通常,?t是個(gè)小量(<0.05 fs),離子位置bi在t1到t2的時(shí)間段內(nèi)幾乎不變,所以對Sk(t)能夠保證Sk(t′)=Sk(t2)這一假設(shè)的合理性.然而,由于電子態(tài)迅速演化,Hk(t)在不同時(shí)刻有很大差異,此時(shí)有兩個(gè)選擇:

目前已有文獻(xiàn)指出,哈密頓算符在0.2—0.5 fs幾乎線性變化,所以當(dāng)?t<0.05 fs時(shí)后者給出的結(jié)果更加準(zhǔn)確[17,18].另外,對于TDKS方程,當(dāng)采用ν表示時(shí),盡管unk(r)(t2)并不顯式地依賴其他的TDKS軌道un′k′(r)(t1)(n′?=n或者k′?=k),但是由于Hk由電荷密度決定,為所有被占據(jù)軌道的加權(quán)求和,所以能帶間的散射已經(jīng)被考慮[1].

數(shù)值實(shí)現(xiàn)上,可用一階Crank-Nicholson形式展開方程(10)中的演化算符

在技術(shù)層面上,方程(13)的計(jì)算速度可以通過兩方面加速:首先,采用ScaLAPACK做并行的矩陣求逆以及矩陣乘法;其次,為了縮短計(jì)算時(shí)間,的數(shù)值只在離子位置bi改變時(shí)才進(jìn)行更新,所以,當(dāng)離子位置固定時(shí),只在第一個(gè)離子步計(jì)算

2.3 電荷密度以及自洽計(jì)算

電荷密度矩陣ρiα,jβ(t2)可通過方程(10)計(jì)算出的unk(r,t2)得到:

其中,qn,k為第n個(gè)能帶在動(dòng)量空間k點(diǎn)的占據(jù)數(shù);cn,jβ,k(t2)為數(shù)值原子軌道基中各個(gè)基矢對unk(r,t2)的貢獻(xiàn),

在波函數(shù)演化過程中,采用電子密度自洽演化能夠顯著地提高計(jì)算結(jié)果的穩(wěn)定性[19].以電荷密度差作收斂判定標(biāo)準(zhǔn)為例,當(dāng)

且η≈10?4時(shí)可認(rèn)為達(dá)到收斂標(biāo)準(zhǔn).若未達(dá)到,則需要使用密度矩陣的線性組合來生成下一個(gè)循環(huán)的輸入量ρnext:

其中ρ,ρnew分別為輸入、輸出的密度矩陣;w為混合比例,通常為0.1—0.5.值得指出的是,TDAP與SIESTA軟件開發(fā)的所有收斂測試標(biāo)準(zhǔn)都兼容.

2.4 數(shù)據(jù)處理

有效數(shù)據(jù)的提取將在自洽過程結(jié)束后進(jìn)行,其中包括對總能量、Hellmann-Feynman力等物理量的計(jì)算.本文只介紹部分獨(dú)特的分析方法.

首先,本文方法可以計(jì)算不同能態(tài)之間躍遷的概率[20]:

其中vnk為絕熱基矢,

絕熱態(tài)nk上的電子占據(jù)數(shù)qnk由TDKS軌道上的投影得到:

其中nk,occ為在k點(diǎn)的占據(jù)態(tài).

其次,對于不同種類的模擬體系采用不同的響應(yīng)函數(shù),代表系統(tǒng)對外界刺激(比如光場)的響應(yīng).對于有限大小的結(jié)構(gòu),在有真空層的方向計(jì)算與時(shí)間相關(guān)的電偶極矩:

其中ρ(r)為電荷密度.然而,對于具有周期性的體系,偶極矩的定義失效,此時(shí),則采用含時(shí)的微觀電流作為響應(yīng)函數(shù):

3 計(jì)算實(shí)例

目前,我們已將基于原子軌道基的rt-TDDFT方法應(yīng)用到一些典型系統(tǒng)的超快動(dòng)力學(xué)研究中.研究重點(diǎn)在于低維體系中電子-光子、聲子-光子及電子-聲子的相互作用.下面給出應(yīng)用TDAP計(jì)算的幾個(gè)具體例子,包括石墨烯納米帶的光吸收性質(zhì)、超快激光誘導(dǎo)的硅晶體非熱熔化、MoS2/WS2界面間的電荷傳輸以及雙色光調(diào)制下的高次諧波產(chǎn)生.

3.1 石墨烯納米帶的光吸收譜及電子動(dòng)力學(xué)模擬

首先介紹如何利用rt-TDDFT方法計(jì)算椅式石墨烯納米帶(armchair graphene nanoribbons,AGNRs)的光吸收,并根據(jù)動(dòng)量分辨的優(yōu)勢監(jiān)控其激發(fā)細(xì)節(jié).對于AGNRs[21,22],由于納米帶平面存在真空層,所以可以采用光場的長度規(guī)范.同時(shí),沿著μ∈(x,y,z)方向的光吸收強(qiáng)度可定義為介電函數(shù)的虛部:

其中αμ,υ代表在頻域偶極矩Pμ(ω)對電場Eυ(ω)的響應(yīng),

在rt-TDDFT的計(jì)算中,首先向體系施加電場來獲得偶極矩在時(shí)域上的變化,然后進(jìn)行傅里葉變換:

其中Eυ(t)可選作單位階躍函數(shù),

如圖2(a)所示,Eυ(t)沿著納米帶的寬度方向.圖2(b)顯示了不同寬度納米帶從紅外到紫外光波段(1—7 eV)的吸收光譜.我們觀察到在不同的能量范圍內(nèi)存在兩種不同的激發(fā)模式,其峰值分別位于可見光波段(1—4 eV)以及紫外波段(6—7 eV).隨著納米帶寬度的增加,這兩種模式都有明顯的紅移和收斂趨勢,這與已有的實(shí)驗(yàn)事實(shí)高度一致[23,24].對于紫外波段的激發(fā)模式,單調(diào)的紅移(隨著寬度從6.064 ?增加到60.36 ?,峰位由6.7 eV紅移至5.6 eV)可以外推到單層石墨烯的π—π躍遷(紅色虛線)[25,26],這說明了該演變對應(yīng)于一維納米帶向二維石墨烯能帶結(jié)構(gòu)的轉(zhuǎn)變.這與使用GW準(zhǔn)粒子近似但是不考慮激子效應(yīng)計(jì)算石墨烯光吸收譜得到的結(jié)果相符合[25],表明rt-TDDFT方法具有很高的精確性.

另一方面,對于可見光區(qū)域的激發(fā)模式,當(dāng)納米帶的寬度達(dá)到無窮大時(shí),尖銳的吸收峰將擴(kuò)展為一個(gè)平滑等高的吸收區(qū)域,對應(yīng)于單層石墨烯2.3%的吸收平臺(tái)區(qū)[27].同時(shí),以上現(xiàn)象在電場強(qiáng)度變化時(shí)仍保持不變,說明了兩種激發(fā)模式的普遍性.

為了得到更細(xì)致的動(dòng)力學(xué)信息,借助TDAP對動(dòng)量的分辨能力,圖3展示了寬度為6.064 ?的納米帶的不同吸收模式在激發(fā)過程中,能帶上不同k點(diǎn)的電子占據(jù)情況.對于可見光模式,電子激發(fā)主要集中在個(gè)別獨(dú)立的k點(diǎn)上.相比之下,高能吸收模式中,電子的激發(fā)趨向于集體性質(zhì)的k點(diǎn)分布.因此,可見光吸收模式有著個(gè)體激發(fā)的特征,而高能模式則包含更多集體激發(fā)的特點(diǎn).

圖2 (a)光場沿著寬度方向激發(fā)椅型納米帶;(b)不同寬度納米帶的光吸收譜Fig.2.(a)Schematic showing AGNRs under an external electric fi eld polarized long the ribbon width;(b)optical absorption spectra of AGNRs with di ff erent sizes.

圖3 兩種激發(fā)模式下的電子占據(jù) (a)可見光區(qū);(b)紫外光區(qū)Fig.3.Snapshot of electron population distribution of two modes:(a)Visible light(3.9 eV);(b)ultraviolet(6.73 eV)excitation.

3.2 超快激光誘導(dǎo)的單晶硅的非熱熔化

激光激發(fā)能夠產(chǎn)生諸多超快現(xiàn)象以及新的物質(zhì)凝聚態(tài)[28].在皮秒時(shí)間尺度下,光激勵(lì)導(dǎo)致的熔化過程已經(jīng)在很多半導(dǎo)體[29?31]和二維材料[32,33]中被觀察到.然而,超快熔化的原子機(jī)制仍然處于爭議之中.由此,我們模擬了在實(shí)驗(yàn)條件下硅的激光熔化[33,34].基于含時(shí)密度泛函以及分子動(dòng)力學(xué)相結(jié)合的模擬結(jié)果表明,在沒有激光時(shí),由于熱擾動(dòng),晶格溫度在300 K左右振蕩.我們用激發(fā)到導(dǎo)帶上電子的百分比η來代表激光強(qiáng)度,光場持續(xù)時(shí)間為25 fs.圖4(a)和圖4(b)顯示了在η=10.16%以及η=0%兩種情況下由激光誘導(dǎo)的結(jié)構(gòu)差異.圖4(b)展現(xiàn)的結(jié)構(gòu)畸變顯然由光誘導(dǎo)并顯示出熔化的特征.圖4(c)—(f)為在光場結(jié)束后,Si—Si鍵的徑向分布函數(shù)在t=25,50,100,300 fs不同時(shí)刻的分布.可以看出第一個(gè)峰值逐漸向右移動(dòng),這意味著最近鄰距離的增加;同時(shí)所有的峰均變得更寬,表征晶體結(jié)構(gòu)正在被破壞.

為了進(jìn)行更定量的描述,我們采用Lindemann熔化標(biāo)準(zhǔn)(Lindemann criterion):當(dāng)硅原子平均平方位移的根ment,RMSD)大于臨界點(diǎn)Rc=0.35 ?時(shí)則認(rèn)為熔化發(fā)生[35].從圖5(a)可以看出:在沒有激光(η=0%)時(shí),RMSD的極大值僅為Rc的一半;然而,當(dāng)激光的強(qiáng)度為η=10.16%時(shí),RMSD在100 fs內(nèi)就到達(dá)Rc并且持續(xù)增大,顯示出很明顯的熔化行為.另一方面,RMSD與衍射強(qiáng)度I(t)通過Debye-Waller方程相關(guān)聯(lián):

其中,Q為被探測晶體平面的倒易晶格矢量;?u2(t)?為平均平方位移,即為RMSD的平方.在圖5(b)中,我們模擬(220)晶面的Iη=10.16%(t)以及Iη=0%(t),同時(shí)將已有實(shí)驗(yàn)數(shù)據(jù)列出作為比較.可以看出,模擬結(jié)果和實(shí)驗(yàn)數(shù)據(jù)的特征基本一致:Iη～11%(t)的數(shù)值在硅熔化后降低至0.2,然而Iη～0%(t)只在0.95附近振蕩.由于在整個(gè)模擬過程中晶格的溫度遠(yuǎn)遠(yuǎn)小于熱熔化溫度(圖5(c)),說明該過程為激光誘導(dǎo)的非熱熔化.非熱熔化可以歸因于電子從成鍵態(tài)到反鍵態(tài)之間的轉(zhuǎn)移,導(dǎo)致隨后熔化勢壘的降低.

圖4 硅的非熱熔化[34] (a),(b)不同光場強(qiáng)度下原子結(jié)構(gòu)的變化;(c)—(f)徑向分布函數(shù)的時(shí)間演化Fig.4.Nonthermal melting of silicon[34]:(a),(b)The atomic structures during laser melting;(c)–(f)time evolution of radial distrubution function under two laser intensity η.

圖5 (a)RMSD隨時(shí)間的變化;(b)硅(220)晶面衍射強(qiáng)度隨時(shí)間的變化;(c)晶體的溫度變化[34]Fig.5.(a)RMSD as a function of time;(b)simulated and experimental electron di ff raction intensity of(220)re fl ection as a function of time;(c)time evolution of ionic temperature[34].

3.3MoS2/WS2界面間的電荷傳輸

激光誘導(dǎo)的層間超快的電荷傳輸在決定電荷分離的速度和效率方面起著非常重要的作用,為光電器件的應(yīng)用提供了新的平臺(tái)[36,37].如圖6(a)所示,以MoS2/WS2雙層體系為例,展示該方法如何被用來監(jiān)控這種范德瓦耳斯異質(zhì)結(jié)構(gòu)的電子動(dòng)力學(xué)[38].通過選擇光子的能量恰好等于MoS2的帶隙大小,電子將從特定的態(tài)上激發(fā)并在MoS2留下空穴,隨后空穴會(huì)進(jìn)一步傳遞到WS2.通過同時(shí)求解TDKS方程以及離子的運(yùn)動(dòng),可以得到界面間電子的動(dòng)力學(xué)演化.結(jié)果表明,通過施加扭轉(zhuǎn)、平移、改變層間距等結(jié)構(gòu)調(diào)制方法,能夠非常有效地調(diào)節(jié)電荷傳輸?shù)膭?dòng)力學(xué)過程.在圖6(b)和圖6(c)中,AB1-2H與AA1-3R均是MoS2/WS2雙層體系的穩(wěn)定結(jié)構(gòu),且層間距離相同(6.3 ?),但兩層間旋轉(zhuǎn)180?[39].統(tǒng)計(jì)受光場激發(fā)后WS2軌道空穴(電子)密度的時(shí)間演化χ,并用指數(shù)函數(shù)χ=a+b×exp(?t/τ)擬合[40],其中τ為電荷傳輸?shù)膲勖?如圖6(b)和圖6(c)所示,τAB1及τAA1分別在100 fs及1000 fs的時(shí)間尺度.進(jìn)一步的分析表明,層間的電荷轉(zhuǎn)移只與某些特定層間態(tài)之間的耦合相關(guān),在該情形下,為布里淵區(qū)的K點(diǎn)的|?2?及|?1?態(tài).用偶極躍遷矩陣元M=??2|?Z|?1?來衡量這兩個(gè)態(tài)之間的耦合強(qiáng)度,其中?Z為沿著垂直MX2平面方向的位置算符.發(fā)現(xiàn)MAB1?MAA1且1/τ指數(shù)依賴于M.根據(jù)以上發(fā)現(xiàn)就能夠從微觀機(jī)制上理解兩種結(jié)構(gòu)間電荷傳輸動(dòng)力學(xué)過程的差異,并且為控制電荷傳輸?shù)牧孔觿?dòng)力學(xué)過程提供了新的途徑,從而促進(jìn)二維異質(zhì)結(jié)構(gòu)在新奇光電器件中的應(yīng)用.

圖6 MoS2/WS2雙層之間的超快電荷傳輸[38] (a)MoS2/WS2雙層結(jié)構(gòu)的側(cè)視圖;(b),(c)兩種MoS2/WS2結(jié)構(gòu)下電荷傳輸過程,內(nèi)插圖為對應(yīng)的原子結(jié)構(gòu)Fig.6.Ultrafast charge transfer in MoS2/WS2bilayer[38]:(a)Side view of MoS2/WS2bilayer;(b),(c)hole transfer dynamics for MoS2/WS2in two stacking modes,the insets show the schematic atomic structure.

3.4 雙色光調(diào)制的高次諧波產(chǎn)生

強(qiáng)激光輻照下,固體能夠產(chǎn)生高度非線性的電子和光學(xué)行為.在非線性光學(xué)中,一個(gè)最基本和重要的方面即是在凝聚態(tài)系統(tǒng)中的高次諧波產(chǎn)生(high harmonic generation,HHG)[41?44]. 通過HHG可以重構(gòu)電子運(yùn)動(dòng)的微觀圖像并產(chǎn)生更高能的超短脈沖.然而,HHG的應(yīng)用卻一直受限于其較低的產(chǎn)率,所以如何提高脈沖的強(qiáng)度已成為當(dāng)今超快科學(xué)領(lǐng)域中非常重要的課題.

如圖7(a),采用雙色光來調(diào)控單層MoS2中電子的動(dòng)力學(xué)過程.光與MoS2之間的相互作用通過求解速度規(guī)范下的TDKS方程(方程(6))得到.雙色激光場表示為高斯波包包絡(luò)的正弦函數(shù):

圖7 (a)模擬體系示意圖;(b)不同r2下的電場波形和(c)矢勢波形Fig.7.(a)Schematic showing high harmonic generation in monolayer MoS2under two-color fi elds;(b)the electric fi eld and(c)the vector potential of four representive laser pulses with various r2.

其中,Ei,φi及ωi(i=1,2)為兩束光的電場峰值大小、相位及頻率;兩者之間的相對強(qiáng)度比用r2表示,r2=E2/E1.在我們的模擬中,激光的極化方向沿著MoS2的zigzag方向并且設(shè)定E1=0.07 V/?,ω1=1/3Egap,ω2=1/2Egap(Eg=1.94 eV),φi=0,脈沖峰值時(shí)間t0=30 fs,半高寬σ=8.5 fs.在圖7(b)和圖7(c)給出了不同r2下的激光波形.

圖8 不同r2下的高次諧波產(chǎn)生譜Fig.8.HHG spectra under two-color fi elds with di ff erent relative intensity ratio r2between two laser pulses.

圖8所示的HHG譜可以通過方程(7)及(22)計(jì)算出的電流J(t)做傅里葉變換得到:可以很明顯地看出,通過簡單的控制第二束光的強(qiáng)度,在適當(dāng)?shù)碾p色場中,高階諧波的強(qiáng)度在E>10 eV的能量范圍內(nèi)可以提高10—100倍.與此對比,輸入光的總強(qiáng)度在峰值處的變化僅為40%(0.0525 V/?).同時(shí),平臺(tái)區(qū)的截?cái)嗄芤驳玫搅藬U(kuò)展,提供了生成孤立的阿秒脈沖的可能[45,46].以上現(xiàn)象歸因于第二束激光脈沖的參與可以獲得更大的矢勢.

4 結(jié) 論

本文介紹了一種自主發(fā)展的基于數(shù)值原子軌道基的rt-TDDFT模擬方法,用于對凝聚物質(zhì)激發(fā)態(tài)的動(dòng)力學(xué)過程進(jìn)行大規(guī)模、精確的模擬.電磁場的長度和速度規(guī)范的有效使用,保證了該方法在處理固體、界面以及二維材料中各類超快過程時(shí)的靈活性和可靠性.作為實(shí)例,把該方法應(yīng)用于幾個(gè)典型的體系,如石墨烯納米帶的光吸收性質(zhì)、單晶硅的非熱熔化、MoS2/WS2界面間的電荷傳輸、雙色光調(diào)制的高次諧波產(chǎn)生等,它們均代表著超快光學(xué)領(lǐng)域中的研究熱點(diǎn).通過模擬可以得到與實(shí)驗(yàn)事實(shí)和其他高精度計(jì)算方法高度一致的結(jié)果,并且能從微觀電子結(jié)構(gòu)的角度解釋其工作機(jī)理.我們希望該方法可以擴(kuò)展到更廣闊的研究領(lǐng)域,成為物理學(xué)家、化學(xué)家和材料學(xué)家研究工作中的有效工具.

[1]Runge E,Gross E K 1984 Phys.Rev.Lett.52 997

[2]Ullrich C A 2011 Time-Dependent Density-Functional Theory:Concepts And Applications(Oxford:Oxford University Press)

[3]Sato S,Yabana K,Shinohara Y,Otobe T,Lee K M,Bertsch G 2015 Phys.Rev.B 92 205413

[4]Takimoto Y,Vila F,Rehr J 2007 J.Chem.Phys.127 154114

[5]Sánchez-Portal D,Hernandez E 2002 Phys.Rev.B 66 235415

[6]Lopata K,Govind N 2011 J.Chem.Theory Comput.7 1344

[7]Yabana K,Sugiyama T,Shinohara Y,Otobe T,Bertsch G 2012 Phys.Rev.B 85 045134

[8]Castro A,Werschnik J,Gross E K 2012 Phys.Rev.Lett.109 153603

[9]Yost D C,Yao Y,Kanai Y 2017 Phys.Rev.B 96 115134

[10]Andrade X,Strubbe D,de Giovannini U,Larsen A H,Oliveira M J,Alberdi-Rodriguez J,Varas A,Theophilou I,Helbig N,Verstraete M J 2015 Phys.Chem.Chem.Phys.17 31371

[11]Sato S A,Yabana K 2014 J.Adv.Simulat.Sci.Eng.1 98

[12]Meng S,Kaxiras E 2008 J.Chem.Phys.129 054110

[13]Ma W,Zhang J,Yan L,Jiao Y,Gao Y,Meng S 2016 Comp.Mater.Sci.112 478

[14]Soler J M,Artacho E,Gale J D,García A,Junquera J,Ordejón P,Sánchez-Portal D 2002 J.Phys.Condens.Matter 14 2745

[15]Ordejón P,Artacho E,Soler J M 1996 Phys.Rev.B 53 R10441

[16]Yabana K,Nakatsukasa T,Iwata J I,Bertsch G 2006 Phys.Status Solidi(b)243 1121

[17]Wang Z,Li S S,Wang L W 2015 Phys.Rev.Lett.114 063004

[18]Ren J,Vukmirovi? N,Wang L W 2013 Phys.Rev.B 87 205117

[19]Ren J,Kaxiras E,Meng S 2010 Mol.Phys.108 1829

[20]Rohringer N,Peter S,Burgd?rfer J 2006 Phys.Rev.A 74 042512

[21]Son Y W,Cohen M L,Louie S G 2006 Phys.Rev.Lett.97 216803

[22]Yang L,Park C H,Son Y W,Cohen M L,Louie S G 2007 Phys.Rev.Lett.99 186801

[23]Yan H,Li X,Chandra B,Tulevski G,Wu Y,Freitag M,Zhu W,Avouris P,Xia F 2012 Nat.Nanotech.7 330

[24]Ju L,Geng B,Horng J,Girit C,Martin M,Hao Z,Bechtel H A,Liang X,Zettl A,Shen Y R 2011 Nat.Nanotech.6 630

[25]Yang L,Deslippe J,Park C H,Cohen M L,Louie S G 2009 Phys.Rev.Lett.103 186802

[26]Trevisanutto P E,Holzmann M,C?té M,Olevano V 2010 Phys.Rev.B 81 121405

[27]Gomez C V,Pisarra M,Gravina M,Pitarke J M,Sindona A 2016 Phys.Rev.Lett.117 116801

[28]Ostrikov K K,Beg F,Ng A 2016 Rev.Mod.Phys.88 011001

[29]Shank C,Yen R,Hirlimann C 1983 Phys.Rev.Lett.50 454

[30]Harb M,Ernstorfer R,Hebeisen C T,Sciaini G,Peng W,Dartigalongue T,Eriksson M A,Lagally M G,Kruglik S G,Miller R D 2008 Phys.Rev.Lett.100 155504

[31]Sokolowski-Tinten K,Blome C,Dietrich C,Tarasevitch A,von Hoegen M H,von der Linde D,Cavalleri A,Squier J,Kammler M 2001 Phys.Rev.Lett.87 225701

[32]Porer M,Leierseder U,Ménard J M,Dachraoui H,Mouchliadis L,Perakis I,Heinzmann U,Demsar J,Rossnagel K,Huber R 2014 Nat.Mater.13 857

[33]Hellmann S,Beye M,Sohrt C,Rohwer T,Sorgenfrei F,Redlin H,Kall?ne M,Marczynski-Bühlow M,Hennies F,Bauer M 2010 Phys.Rev.Lett.105 187401

[34]Lian C,Zhang S,Meng S 2016 Phys.Rev.B 94 184310[35]Zijlstra E S,Kalitsov A,Zier T,Garcia M E 2013 Phys.Rev.X 3 011005

[36]Falke S M,Rozzi C A,Brida D,Maiuri M,Amato M,Sommer E,de Sio A,Rubio A,Cerullo G,Molinari E 2014 Science 344 1001

[37]Rozzi C A,Falke S M,Spallanzani N,Rubio A,Molinari E,Brida D,Maiuri M,Cerullo G,Schramm H,Christoffers J 2013 Nat.Commun.4 1602

[38]Zhang J,Hong H,Lian C,Ma W,Xu X,Zhou X,Fu H,Liu K,Meng S 2017 Adv.Sci.4 1700086

[39]van der Zande A M,Kunstmann J,Chernikov A,Chenet D A,You Y,Zhang X,Huang P Y,Berkelbach T C,Wang L,Zhang F 2014 Nano Lett.14 3869

[40]Long R,Prezhdo O V 2016 Nano Lett.16 1996

[41]Ndabashimiye G,Ghimire S,Wu M,Browne D A,Schafer K J,Gaarde M B,Reis D A 2016 Nature 534 520

[42]Luu T T,Garg M,Kruchinin S Y,Moulet A,Hassan M T,Goulielmakis E 2015 Nature 521 498

[43]Vampa G,Hammond T J,Thire N,Schmidt B E,Legare F,McDonald C R,Brabec T,Corkum P B 2015 Nature 522 462

[44]Liu H,Li Y,You Y S,Ghimire S,Heinz T F,Reis D A 2016 Nat.Phys.13 262

[45]Li J B,Zhang X,Yue S J,Wu H M,Hu B T,Du H C 2017 Opt.Express 25 18603

[46]Shiner A D,Trallero-Herrero C,Kajumba N,Bandulet H C,Comtois D,Legare F,Giguere M,Kie ff er J C,Corkum P B,Villeneuve D M 2009 Phys.Rev.Lett.103 073902