基于原子軌道基的實時密度泛函理論:方法及應用?

關夢雪1) 廉超1) 孟勝1)2)?

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

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

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

1 引 言

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

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

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

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

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

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

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

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

2 計算方法

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

2.1 哈密頓量及交迭矩陣

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

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

在每個k點的交迭矩陣Sk及哈密頓量Hk用原子軌道基表示為

哈密頓算符為

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

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

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

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

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

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

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

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

2.2 電子態(tài)演化

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

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

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

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

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

2.3 電荷密度以及自洽計算

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

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

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

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

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

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

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

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

其中vnk為絕熱基矢,

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

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

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

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

3 計算實例

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

3.1 石墨烯納米帶的光吸收譜及電子動力學模擬

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

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

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

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

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

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

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

圖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 超快激光誘導的單晶硅的非熱熔化

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

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

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

圖4 硅的非熱熔化[34] (a),(b)不同光場強度下原子結(jié)構的變化;(c)—(f)徑向分布函數(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隨時間的變化;(b)硅(220)晶面衍射強度隨時間的變化;(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界面間的電荷傳輸

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

圖6 MoS2/WS2雙層之間的超快電荷傳輸[38] (a)MoS2/WS2雙層結(jié)構的側(cè)視圖;(b),(c)兩種MoS2/WS2結(jié)構下電荷傳輸過程,內(nèi)插圖為對應的原子結(jié)構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)生

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

如圖7(a),采用雙色光來調(diào)控單層MoS2中電子的動力學過程.光與MoS2之間的相互作用通過求解速度規(guī)范下的TDKS方程(方程(6))得到.雙色激光場表示為高斯波包包絡的正弦函數(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)為兩束光的電場峰值大小、相位及頻率;兩者之間的相對強度比用r2表示,r2=E2/E1.在我們的模擬中,激光的極化方向沿著MoS2的zigzag方向并且設定E1=0.07 V/?,ω1=1/3Egap,ω2=1/2Egap(Eg=1.94 eV),φi=0,脈沖峰值時間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(t)做傅里葉變換得到:可以很明顯地看出,通過簡單的控制第二束光的強度,在適當?shù)碾p色場中,高階諧波的強度在E>10 eV的能量范圍內(nèi)可以提高10—100倍.與此對比,輸入光的總強度在峰值處的變化僅為40%(0.0525 V/?).同時,平臺區(qū)的截斷能也得到了擴展,提供了生成孤立的阿秒脈沖的可能[45,46].以上現(xiàn)象歸因于第二束激光脈沖的參與可以獲得更大的矢勢.

4 結(jié) 論

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

[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