基于有限元法的光子并矢格林函數(shù)重整化及其在自發(fā)輻射率和能級(jí)移動(dòng)研究中的應(yīng)用?

趙運(yùn)進(jìn) 田錳 黃勇剛? 王小云 楊紅 米賢武

1)(吉首大學(xué)物理與機(jī)電工程學(xué)院,吉首 416000)

2)(懷化學(xué)院電氣與信息工程學(xué)院,懷化 418000)

(2018年5月6日收到；2018年7月19日收到修改稿)

1 引 言

根據(jù)量子電動(dòng)力學(xué)[1?3],自發(fā)輻射和能級(jí)移動(dòng)來(lái)源于量子點(diǎn)與電磁環(huán)境的相互作用.近年來(lái),人們?cè)O(shè)計(jì)了各種人工微納結(jié)構(gòu)來(lái)調(diào)控它們,如光子晶體[4?13]、金屬納米結(jié)構(gòu)[14?21]等. 對(duì)于表面等離激元金屬納米結(jié)構(gòu),在金屬和電介質(zhì)界面處可產(chǎn)生表面等離激元共振,能將光場(chǎng)壓縮到納米尺度,極大地改變了附近量子點(diǎn)的自發(fā)輻射和能級(jí)移動(dòng)特性,譬如:單分子熒光增強(qiáng)[14]、表面等離激元增強(qiáng)發(fā)光二極管[15]、單分子探測(cè)[16,17]、等離激元納米激光器[18,19]、拉曼散射增強(qiáng)[20,21]等.對(duì)于量子點(diǎn)距離金屬表面很近的情況,相互作用的強(qiáng)度甚至能夠超過(guò)輻射和非輻射的損耗,使系統(tǒng)進(jìn)入強(qiáng)耦合區(qū)域,導(dǎo)致了很多新奇的量子現(xiàn)象[22?25].

理論上,自發(fā)輻射率及能級(jí)移動(dòng)可用并矢格林函數(shù)來(lái)表達(dá)[26,27](即Γ(r；ω)=2d2?[n.G(r；r；ω).n]/~ε0；?(r；ω)= ?d2?[n.G(r；r；ω).n]/~ε0,其中,Γ(r；ω)為自發(fā)輻射率,?(r；ω)為能級(jí)移動(dòng)).當(dāng)源點(diǎn)與場(chǎng)點(diǎn)的位置相同時(shí),格林函數(shù)的實(shí)部是一個(gè)發(fā)散的量,數(shù)值計(jì)算會(huì)遇到困難.為處理這樣發(fā)散的非物理的能級(jí)移動(dòng)問(wèn)題,一種方法是將格林函數(shù)分解為均勻部分(發(fā)散,可解析處理[28?34])和散射部分,將均勻空間的貢獻(xiàn)歸結(jié)于躍遷頻率,而只考慮散射場(chǎng)的貢獻(xiàn),用散射格林函數(shù)來(lái)代替總的格林函數(shù),該方法是最自然的重整化方式.2018年,我們利用基于有限元的COMSOL Multiphysics軟件,通過(guò)計(jì)算有納米結(jié)構(gòu)和均勻空間兩種情況下偶極子輻射場(chǎng)的差,準(zhǔn)確獲得了任意微納結(jié)構(gòu)中的散射格林函數(shù)[26].值得注意的是點(diǎn)電偶極子輻射場(chǎng)依賴(lài)于人為設(shè)定的網(wǎng)格大小,因此,在兩次計(jì)算輻射場(chǎng)時(shí),只有采用相同的網(wǎng)格,才能消除發(fā)散的部分.

數(shù)值技術(shù)上,不僅有限元方法會(huì)導(dǎo)致點(diǎn)電偶極子輻射場(chǎng)依賴(lài)于網(wǎng)格大小,時(shí)域有限差分方法也有類(lèi)似的問(wèn)題.2012年,van Vlack等[35]利用時(shí)域有限差分法研究耗散系統(tǒng)中的自發(fā)輻射問(wèn)題,證實(shí)該方法得到的偶極子輻射場(chǎng)是點(diǎn)電偶極子所在的網(wǎng)格內(nèi)的平均場(chǎng).物理上,該發(fā)散來(lái)源于點(diǎn)電偶極子模型,考慮量子點(diǎn)的實(shí)際大小,在量子點(diǎn)體積內(nèi)取平均的重整化方式被廣泛采用[34?36].時(shí)域有限差分法采用直角矩形網(wǎng)格,其差分格式會(huì)不可避免地導(dǎo)致臺(tái)階近似,相比之下,有限元法可以應(yīng)用在更復(fù)雜和更精細(xì)的納米結(jié)構(gòu)中,其網(wǎng)格劃分方式更加靈活,能更好地?cái)M合曲面結(jié)構(gòu),可精確地仿真電磁場(chǎng)的急劇變化.另外,我們已經(jīng)證實(shí)[26],利用有限元可求解任意微納結(jié)構(gòu)中的散射格林函數(shù)問(wèn)題.基于以上考慮,如果有限元方法能準(zhǔn)確處理均勻空間格林函數(shù)的重整化問(wèn)題,則可利用該方法準(zhǔn)確處理任意微納結(jié)構(gòu)中格林函數(shù)的重整化問(wèn)題.不同于散射格林函數(shù)方法,該方法只需一次仿真,能極大地減少計(jì)算量.

本文介紹了自發(fā)輻射和能級(jí)移動(dòng)相關(guān)的理論和格林函數(shù)重整化及其有限元計(jì)算方法；在均勻空間和納米球系統(tǒng)中利用基于有限元的COMSOL Multiphysics軟件,證實(shí)該算法的可應(yīng)用性和準(zhǔn)確性；將該方法應(yīng)用到球-平板納米腔系統(tǒng),研究自發(fā)輻射和能級(jí)移動(dòng)問(wèn)題.

2 理論及算法驗(yàn)證

2.1 自發(fā)輻射率和能級(jí)移動(dòng)的并矢格林函數(shù)表示

根據(jù)量子電動(dòng)力學(xué)[33],耗散吸收介質(zhì)中二能級(jí)原子與電磁場(chǎng)相互作用,電偶極近似下Hamiltonian可寫(xiě)為:

其中,H0是未相互作用部分,包含電磁場(chǎng)部分和原子部分；HI是原子與場(chǎng)之間的相互作用；d=?g|d?|e?=dn為原子的躍遷偶極子矩陣元,其大小和方向分別為d和n；(r,ω)和(r,ω)表示電磁場(chǎng)和介質(zhì)的產(chǎn)生和湮滅場(chǎng)算符,滿(mǎn)和r0為原子的躍遷頻率和位置；|e?和|g?分別表示原子的激發(fā)態(tài)和基態(tài).考慮局域響應(yīng)介質(zhì)[22?25],當(dāng)相對(duì)介電函數(shù)ε(r,ω)=εR(r′,ω)+iεI(r′,ω)的實(shí)部和虛部滿(mǎn)足Kramers-Kronig關(guān)系,相對(duì)磁導(dǎo)率為1時(shí),電場(chǎng)算符可以表示為:

其中,ε0為真空介電常數(shù),光子并矢格林函數(shù)G(r,r′,ω)滿(mǎn)足[26,27]:

I為單位并矢.

假設(shè)原子初始時(shí)刻處于激發(fā)態(tài),由于旋波近似,系統(tǒng)的粒子數(shù)守恒,系統(tǒng)相關(guān)的態(tài)|a?和|b(r,ω)?分別為:|a?=|e,0?和|b(r,ω)?=|g?f??(r,ω)|0?.則系統(tǒng)狀態(tài)可寫(xiě)為

為求解a(t)= ?a|U(t)|a?,可利用投影算符方法[37,38]. 利用推遲超前格林函數(shù)G±(ω)=其中,G(z)=(z?H/~)?1,則時(shí)間演化算符為

定義投影算符P=|a??a|及其補(bǔ)空間投影算符Q=I?P,其中I為恒等算子,根據(jù)恒等式(z?H/~)G(z)=I,可得:

由方程(4)可得格林函數(shù)矩陣元Gaa(z)=?a|G(z)|a?=1/(z?ω0?Raa(z)),其中Raa(z)=?a|R(z)|a?.依據(jù)格林函數(shù)理論,格林函數(shù)的奇點(diǎn),即方程z?ω0?Raa(z)=0的解為系統(tǒng)的本征值,其實(shí)部對(duì)應(yīng)著系統(tǒng)能級(jí),虛部對(duì)應(yīng)著衰減快慢.將電場(chǎng)的表達(dá)式((2)式)代入(1)式,利用(5)式,則能移算符矩陣元為:R±aa(ω) ≡ limη→0+?a|R(z=ω±iη)|a?= ?(r0,ω)±Γ(r0,ω)/2.經(jīng)過(guò)進(jìn)一步的推導(dǎo),可得:

其中,?和?分別表示虛部和實(shí)部,d和n仍然為(1)式中躍遷偶極矩陣元d=dn的大小和單位方向矢量.

在弱耦合情況下,即Γ(r；ω)和?(r；ω)遠(yuǎn)小于原子的躍遷頻率ω0,且是頻率的緩變函數(shù)時(shí),Γ(r；ω)和?(r；ω)可近似為Γ=Γ(r；ω0)和?=?(r；ω0),此時(shí),時(shí)間演化算符矩陣元為,

此時(shí),Γ和?表示自發(fā)輻射率和能級(jí)移動(dòng).本文稱(chēng)Γ(r；ω)和?(r；ω)為自發(fā)輻射率和能級(jí)移動(dòng).

值得注意的是,能級(jí)移動(dòng)公式(7)中第二項(xiàng)非常小,本文主要考慮第一項(xiàng),即

此時(shí),自發(fā)輻射率和能級(jí)移動(dòng)的表達(dá)式與文獻(xiàn)[26,27]一致.(6)式和(8)式是本節(jié)的主要結(jié)果,下節(jié)將說(shuō)明其中的格林函數(shù)需用重整化的表達(dá)式代替.

2.2 并矢格林函數(shù)的重整化及其有限元計(jì)算

2.2.1 并矢格林函數(shù)的重整化及其與自發(fā)輻射率和能級(jí)移動(dòng)的關(guān)系

本小節(jié)主要目的是闡明自發(fā)輻射率((6)式)中的并矢格林函數(shù)需用其重整化的表達(dá)式,而能級(jí)移動(dòng)((8)式)中的并矢格林函數(shù)需用重整化的散射格林函數(shù)代替,下面我們具體說(shuō)明.

依電磁理論,并矢格林函數(shù)可用點(diǎn)電偶極子的輻射場(chǎng)表示.r0處點(diǎn)電偶極子P=dδ(r?r0)在r處的場(chǎng)為

為解決以上發(fā)散問(wèn)題,可采用實(shí)腔模型,即用一個(gè)無(wú)耗的小體積?V(理論上通常為球形或者小立方體)表示量子點(diǎn)所在的區(qū)域,其相對(duì)介電常數(shù)為εB,量子點(diǎn)的躍遷偶極矩位于該小體積的中心r0,用重整化的電場(chǎng),即?V中平均局域場(chǎng)表示量子點(diǎn)所感受到的場(chǎng)[35].

類(lèi)似于文獻(xiàn)[36]的處理,設(shè)r0處?V內(nèi)的極化強(qiáng)度為P(r；ω),則?V內(nèi)r處的電場(chǎng)可用散射格林函數(shù)Gscatt(r,r′；ω)和均勻空間格林函數(shù)Ghom(r,r′；ω)表達(dá)為

其中Vδ為r′處一個(gè)無(wú)限小的體積元,且對(duì)于立方型或者球形體積?V[28],

由于散射格林函數(shù)是空間的緩變函數(shù),(11)式中Gscatt(r,r′；ω) ≈ Gscatt(r,r0；ω),利用此關(guān)系,并將(11)式代入(10)式,則電場(chǎng)在?V內(nèi)的平均值,即重整化電場(chǎng)變?yōu)?/p>

對(duì)于半徑為a的球形體積?V=4πa3/3,波數(shù)

定義重整化的格林函數(shù):

則(13)式,即重整化的電場(chǎng)也可寫(xiě)為

由(11)式可知,偶極子的輻射場(chǎng)為散射場(chǎng)和均勻空間輻射場(chǎng)的和,其中第一項(xiàng)代表散射場(chǎng),第二和第三項(xiàng)代表均勻空間的輻射場(chǎng).重整化后,即(13)式中的第一項(xiàng)為重整化的散射場(chǎng),可寫(xiě)為

第二和第三項(xiàng)為重整化的均勻空間輻射場(chǎng),可寫(xiě)為

其中,重整化均勻空間格林函數(shù)為

由(17)—(19)式,利用d=dn,重整化的電場(chǎng)與重整化的格林函數(shù)分量只相差系數(shù),因此,本文不特別指明的情況下,重整化的電場(chǎng)和重整化的格林函數(shù)具有相同的意義.

對(duì)比電場(chǎng)的重整化形式((17)式)和未重整化的形式((9)式),由于原子感受到的場(chǎng)需由重整化的場(chǎng)替換(此時(shí)作用于點(diǎn)電偶極子的輻射反作用場(chǎng)為ERR=i?[Greg(r0,r0；ω)].d/ε0,該場(chǎng)阻礙偶極子的振蕩,詳細(xì)解釋請(qǐng)參閱文獻(xiàn)[34]),則自發(fā)輻射率中的格林函數(shù)應(yīng)當(dāng)采用重整化的格林函數(shù).此時(shí),(6)式變?yōu)?/p>

與文獻(xiàn)[35]一致.

對(duì)于能級(jí)移動(dòng),通常將均勻空間的貢獻(xiàn)歸結(jié)于量子點(diǎn)的躍遷頻率(均勻空間中的能級(jí)移動(dòng)可解析處理,詳細(xì)處理方法見(jiàn)文獻(xiàn)[3]),本文只考慮散射體的存在對(duì)能級(jí)移動(dòng)的影響.此時(shí),(8)式變?yōu)?/p>

(21)和(22)式是計(jì)算自發(fā)輻射率和能級(jí)移動(dòng)的一般表達(dá).

2.2.2 有限元方法及其在重整化并矢格林函數(shù)計(jì)算中的應(yīng)用

由(17)式可知,Greg(r0,r0；ω)可通過(guò)重整化的偶極子輻射場(chǎng)Ereg(r0；ω)得到.具體的計(jì)算流程如下:首先求解r0處以頻率ω沿n方向振蕩的點(diǎn)電偶極子P=dnδ(r?r0)所產(chǎn)生的輻射場(chǎng)E(r；ω),隨后求該輻射場(chǎng)E(r；ω)在?V內(nèi)的平均值,即計(jì)算重整化的電場(chǎng)Ereg(r0；ω)((10)式),進(jìn)一步,由(17)式可得重整化格林函數(shù)

由(16)式可得重整化的散射格林函數(shù)

其中,對(duì)于球形?V,M和L的表達(dá)式見(jiàn)(15)和(12)式.

由以上可知,計(jì)算Greg(r0,r0；ω)和Gr0,r0；ω)的關(guān)鍵是計(jì)算點(diǎn)電偶極子的輻射場(chǎng)及其在?V內(nèi)的平均值.因此,原則上能計(jì)算點(diǎn)電偶極子重整化輻射場(chǎng)的數(shù)值方法均可用于計(jì)算重整化的格林函數(shù)和散射格林函數(shù).本文擬用有限元法計(jì)算,下面簡(jiǎn)要介紹有限元算法以及利用基于有限元的COMSOL Multiphysics軟件實(shí)現(xiàn)以上計(jì)算過(guò)程.

有限元方法是一種基于變分法的離散化數(shù)值技術(shù),在電磁場(chǎng)及納米光子學(xué)領(lǐng)域具有廣泛的應(yīng)用[39,40].文獻(xiàn)[41]詳細(xì)介紹了有限元方法計(jì)算線電流源的輻射場(chǎng),為方便,我們簡(jiǎn)要概括其要點(diǎn).依頻率域麥克斯韋方程,電流元J(r,ω)產(chǎn)生的輻射場(chǎng)E(r,ω)滿(mǎn)足:

為仿真無(wú)限大的空間,需將空間截?cái)?在邊界處通常采用吸收邊界條件,使得電磁波入射到邊界處無(wú)反射.為此,通常在計(jì)算邊界外設(shè)定吸收區(qū)域(完全匹配層PML),通過(guò)設(shè)置特定的各向異性吸收材料或者采用復(fù)坐標(biāo)系拉伸形式的波方程,使得電磁波在該區(qū)域無(wú)反射地衰減.

依有限元理論,對(duì)于任意試探函數(shù)F(r,ω),則方程(23)的解等價(jià)于求以下泛函的極值點(diǎn),

由以上可知,有限元方法解的是以上積分方程,能處理數(shù)學(xué)上的點(diǎn)源情況,即J(r,ω)=nδ(r?r0).

為求解方程(24),將空間離散成連續(xù)的子域,然后將解E(r,ω)和試探函數(shù)F(r,ω)用子域上的局域矢量基函數(shù)Nj(r)做展開(kāi),如:

則方程(24)轉(zhuǎn)換成線性方程組:

求解以上線性方程組,可得電場(chǎng)展開(kāi)系數(shù)Ek,進(jìn)而得到原方程的解于Nj(r)可選擇高階函數(shù)形式,能較高精度地處理復(fù)雜的幾何結(jié)構(gòu)[39?41],使得有限元方法特別適用于納米光子學(xué)領(lǐng)域.

利用基于有限元的COMSOL Multiphysics[42]軟件的射頻模塊,以上利用有限元求解點(diǎn)電偶極子輻射場(chǎng)及其平均值的過(guò)程可方便實(shí)現(xiàn).該軟件提供了點(diǎn)電偶極子模塊,考慮長(zhǎng)度為la(la? λ)的線電流I,其偶極矩為d=?iIla/ω(注意,軟件采用時(shí)諧因子eiωt,由J= ?P/?t可得極化強(qiáng)度P(r,ω)=J(r,ω)/iω,進(jìn)而得到

綜上所述,利用有限元計(jì)算重整化格林函數(shù)和重整化散射格林函數(shù)的計(jì)算步驟分為以下四步:1)設(shè)點(diǎn)電偶極子沿n方向,利用類(lèi)似于文獻(xiàn)[43]的方法計(jì)算輻射場(chǎng)；2)將偶極子的輻射場(chǎng)在以點(diǎn)偶極子為中心,半徑為a的小球體內(nèi)做積分得到重整化的輻射場(chǎng)；3)利用(17)式得到重整化的格林函數(shù)分量；4)利用(16)式將重整化格林函數(shù)減去解析的重整化均勻空間格林函數(shù),得到重整化的散射格林函數(shù)分量.

2.2.3 算法驗(yàn)證——均勻空間及金屬納米球系統(tǒng)

為證實(shí)以上重整化格林函數(shù)的有限元計(jì)算方法的可靠性,首先研究重整化格林函數(shù)有解析解的真空情況.依2.2.2小節(jié)計(jì)算步驟可知,計(jì)算重整化格林函數(shù)的關(guān)鍵在于計(jì)算點(diǎn)電偶極子的輻射場(chǎng).為簡(jiǎn)便,設(shè)點(diǎn)電偶極子處于坐標(biāo)原點(diǎn),其極化方向沿著z方向,考慮對(duì)稱(chēng)性,計(jì)算區(qū)域可選為球的四分之一(盡管此處情況可利用二維旋轉(zhuǎn)軸對(duì)稱(chēng)性的求解方法,為不失一般性,本文仍然用三維模型計(jì)算；另外由于對(duì)稱(chēng)性,實(shí)際的偶極矩為軟件設(shè)定的偶極矩的四倍,因此在利用軟件計(jì)算時(shí),點(diǎn)電偶極子的強(qiáng)度設(shè)為Ila=0.25 A.m).詳細(xì)計(jì)算過(guò)程如下:1)建模及求解域剖分,建立如圖1所示的計(jì)算空間,設(shè)中心處點(diǎn)電偶極子振蕩頻率為ω=2.75 eV(此處及以后用ω表示~ω),中間藍(lán)色網(wǎng)格密集區(qū)域?yàn)橹卣男∏蝮w,其半徑為a(采用四面體網(wǎng)格剖分,網(wǎng)格尺寸設(shè)定見(jiàn)表1計(jì)算說(shuō)明),中間空氣層的厚度為波長(zhǎng)λ(采用四面體網(wǎng)格剖分,最大網(wǎng)格尺寸為λ/10),最外層為PML吸收邊界層,其厚度為0.25λ(采用10層掃掠網(wǎng)格),用以截?cái)酂o(wú)限大的空間；2)設(shè)定材料參數(shù)及邊界條件,根據(jù)情況設(shè)定求解域中的材料介電函數(shù),在xz和yz截面上采用完美磁導(dǎo)體邊界條件.通過(guò)COMSOL Mutiphysics求解,可得到點(diǎn)電偶極子輻射場(chǎng)在空間的分布,將其在偶極子源附近小球體內(nèi)做平均,可得重整化的輻射場(chǎng).

圖1 真空重整化格林函數(shù)有限元計(jì)算模型圖 利用對(duì)稱(chēng)性,只仿真四分之一個(gè)球,xz和yz截面上采用完美磁導(dǎo)體邊界條件,點(diǎn)電偶極子源處于球心處,中間網(wǎng)格加密藍(lán)色區(qū)域?yàn)橹卣瘏^(qū)域,在該區(qū)域中,我們考慮兩種不同的網(wǎng)格尺寸Fig.1. Models for calculating the renormalized Green’s function in vacuum.The simulation domain is reduced to one quarter for symmetry consideration.The point dipole source located at the center is thought to be polarized along the z-axis.Perfect match layer(PML)is used to truncate the in finite space and perfect magnetic conductor conditions are used in the xz and yz planes.There is a small sphere region which is used to renormalizing the radiation field.

在偶極子源附近,即重整化區(qū)域內(nèi),采用兩種不同的網(wǎng)格尺寸:一般網(wǎng)格精度(最大網(wǎng)格尺寸為0.4 nm)和精細(xì)網(wǎng)格精度(最大網(wǎng)格尺寸為0.05 nm).偶極子處的輻射場(chǎng)的z方向分量Ez實(shí)部和虛部如表1的第1列和第2列所示,發(fā)現(xiàn)源點(diǎn)處的偶極子輻射場(chǎng)虛部強(qiáng)烈地依賴(lài)于網(wǎng)格,源點(diǎn)處網(wǎng)格尺寸變小時(shí),Ez虛部由?5.7563×1024V/m變?yōu)?7.3745×1026V/m,有兩個(gè)數(shù)量級(jí)的變化(見(jiàn)表1第2列數(shù)據(jù)),Ez實(shí)部幾乎不變.通過(guò)仿真,當(dāng)我們進(jìn)一步減小源點(diǎn)附近的網(wǎng)格尺寸時(shí),輻射場(chǎng)虛部變得越來(lái)越大,而其實(shí)部變化很小,暗含著源點(diǎn)處輻射場(chǎng)的虛部發(fā)散；當(dāng)網(wǎng)格過(guò)于密集(例如,小于0.001 nm)時(shí),不能得到收斂的結(jié)果.這主要是由于輻射場(chǎng)的虛部對(duì)應(yīng)格林函數(shù)的實(shí)部(見(jiàn)(9)式,其中偶極矩為d=iIla/ω),而當(dāng)源點(diǎn)和場(chǎng)點(diǎn)在同一位置時(shí),均勻空間格林函數(shù)的實(shí)部是發(fā)散的.綜上可知,有限元計(jì)算得到的點(diǎn)電偶極子的輻射場(chǎng)依賴(lài)于計(jì)算過(guò)程中設(shè)定的網(wǎng)格尺寸,其值不具有確定的意義.

考慮重整化,即對(duì)以上電場(chǎng)在以偶極子為中心、半徑為a的小球體內(nèi)做平均,當(dāng)a=1 nm,Ezreg(r0；ω)的計(jì)算結(jié)果如表1的第3和第4列所示,發(fā)現(xiàn)重整化的輻射場(chǎng)實(shí)部與未重整化情況幾乎沒(méi)有差別,更重要的是其虛部不再隨網(wǎng)格尺寸的改變而變化,獲得收斂的結(jié)果(見(jiàn)表1第4列).

下面,考慮不同的半徑a(0.05—1 nm),重整化電場(chǎng)實(shí)部和虛部分別如圖2(a)和圖2(b)所示,黑色實(shí)線表示解析解((19)式,其中G(r0,r0；ω)((20)式)中的M和L見(jiàn)(15)和(12)式),紅色空心圓圈表示有限元結(jié)果((10)式,其中E(r,ω)由有限元軟件計(jì)算得到).我們發(fā)現(xiàn):有限元法的數(shù)值結(jié)果與解析解符合得很好.為展示有限元結(jié)果與解析解的差別,在插圖中展示了相對(duì)誤差,即用數(shù)值解與解析解的差除以解析解,發(fā)現(xiàn)輻射場(chǎng)實(shí)部的相對(duì)誤差小于0.004%,虛部的相對(duì)誤差在0.3%左右,進(jìn)一步證實(shí)有限元方法能處理均勻空間中并矢格林函數(shù)的重整化問(wèn)題.

表1 真空中點(diǎn)電偶極子的輻射場(chǎng):1,2列為不重整化的有限元結(jié)果；3,4列為重整化的結(jié)果(依(10)式計(jì)算),其中ω=2.75 eV,a=1 nmTable 1.Radiation field of point electric dipole in homogeneous space.Results are obtained by finite element method with perfect matched layer(COMSOL Multiphysics software).Around the dipole source,two different levels of mesh size are used,for examplefiner mesh grid and general mesh grid.the first two columns arefor electric fields at the dipole location and the last two columns are the renormalized ones over a sphere with radius a=1 nm.The real(imaginary)part of the radiation field corresponds to the imaginary(real)part of the Green function.Here,ω =2.75 eV and the dipole strength is d=1[A.m]/jω.

接下來(lái),展示我們的算法在重整化散射格林函數(shù)求解中的應(yīng)用.量子點(diǎn)位于一個(gè)銀納米球外,模型如圖3所示.設(shè)背景為真空,其介電常數(shù)ε1=1,銀的介電常數(shù)由Drude模型給出[26,27],即:ε2(ω)=ε∞?ω(ω2? iγω),在這里ε∞=6,ωp=7.89 eV,γ=0.051 eV,設(shè)銀納米球半徑ra=20 nm,量子點(diǎn)距離銀納米球表面為l=2 nm,偶極矩方向沿著納米球的徑向,重整化體積仍然用半徑為a的球體表示.

圖3 納米球系統(tǒng)中重整化格林函數(shù)的仿真示意圖 量子點(diǎn)‘A’在銀納米球附近,考慮量子點(diǎn)的實(shí)際大小,用半徑為a的小球代替,量子點(diǎn)的中心距離球表面為l,設(shè)躍遷偶極矩的極化方向沿著銀納米球的徑向；銀納米球半徑設(shè)為ra=20 nm,ε1和ε2分別是空氣的和銀的介電常數(shù)Fig.3.Model of nano-sphere system used to calculate the renormalized Green function.A quantum dot‘A’is placed near a silver nano-sphere,where a represents the radius of the quantum dot and l is its center to the surface of silver nano-sphere.For simplicity,the transition dipole moment is set to be polarized along the radial direction of the nano-sphere and the radius of the nanosphere is ra=20 nm.The relative permittivity for air and silver are ε1=1 and ε2(ω)= ε∞ ? ω/(ω2 ? iγω)with ε∞ =6,ωp=7.89 eV,γ =0.051 eV,respectively.

計(jì)算結(jié)果如圖4(a)和圖4(b)所示.黑色實(shí)線描述的是解析的點(diǎn)電偶極子的散射場(chǎng)(解析散射格林函數(shù)見(jiàn)文獻(xiàn)[26,27]),紅色空心圓圈描述的是本文提出的有限元方法得到的點(diǎn)電偶極子的重整化散射場(chǎng).其求解步驟為:首先求點(diǎn)電偶極子的輻射場(chǎng)E(r,ω),利用(10)式將其重整化,得到Ereg(r0；ω),再利用(17)式得到重整化的格林函數(shù)Greg(r0,r0；ω),隨后利用(18)式(其中,

圖4 納米球系統(tǒng)中的散射場(chǎng) 黑色實(shí)線描述解析的點(diǎn)電偶極子的散射場(chǎng),可通過(guò)解析散射格林函數(shù)或者利用有結(jié)構(gòu)和無(wú)結(jié)構(gòu)情況下兩次有限元仿真的差值得到[26]；紅色空心圓圈代表本文的一次有限元仿真得到的點(diǎn)電偶極子的重整化散射場(chǎng),其中,a=1 nm；可以看到重整化散射場(chǎng)與量子點(diǎn)中心處的散射場(chǎng)符合得很好Fig.4.Scattering field of point electric dipole in nanosphere system. The black solid line is for the analytical result,which is obtained by analytical scattering Green function method or by the difference of two different finite element runs for every frequency point[26].The red circles represent the renormalized scattering fields by this work.Here,the radius of the quantum is a=1 nm.We see that the renormalized scattering fields agree well with the scattering fields at the center of the quantum dot.

不同于我們最近的類(lèi)似的工作[26],此方法只需要一次仿真,通過(guò)計(jì)算重整化的電場(chǎng),得到重整化的格林函數(shù),將其減去解析的均勻空間的重整化格林函數(shù),能得到重整化的散射格林函數(shù),而之前的工作需要兩次有限元仿真,利用有結(jié)構(gòu)和均勻空間兩種情況下偶極子的輻射場(chǎng)的差得到散射場(chǎng),進(jìn)而得到散射格林函數(shù).值得強(qiáng)調(diào)的是,盡管有限元方法能準(zhǔn)確計(jì)算偶極子的輻射場(chǎng),但是對(duì)每一個(gè)頻率點(diǎn),解線性方程組((25)式)所需的計(jì)算資源較多(例如:對(duì)于本例所示的納米球系統(tǒng),采用對(duì)稱(chēng)性后大約需要40 GB內(nèi)存,在我們的塔式工作站“Intel(R)Xeon(R)E5-2697 v3”上需30 min左右).通過(guò)以上對(duì)比,容易看到,相比于之前的方法,本方法不需計(jì)算均勻空間中的輻射場(chǎng),計(jì)算量約減小一半.

3 算法應(yīng)用——球-平板納米腔中的自發(fā)輻射和能級(jí)移動(dòng)

綜上可知,采用點(diǎn)電偶極子源,通過(guò)一次有限元仿真,能精確處理格林函數(shù)和散射格林函數(shù)的重整化問(wèn)題.本節(jié)利用該方法,研究球-平板表面等離激元納米腔中自發(fā)輻射和能級(jí)移動(dòng)問(wèn)題[23,44]自發(fā)輻射和能級(jí)移動(dòng)由(21)和(22)式表示.圖5是球-平板模型的示意圖,無(wú)限大的銀平板上方有一個(gè)銀納米球,其半徑為ra=20 nm,離銀板的距離為g,量子點(diǎn)在銀納米球和銀板的正中間,其偶極矩設(shè)為24D≈0.5 e.nm,設(shè)極化方向沿著垂直于銀板的方向.

圖5 球-平板納米腔系統(tǒng)的原理圖 ra=20 nm的銀納米球置于無(wú)限大的銀板上方,銀納米球與銀板之間的距離是g,量子點(diǎn)放在球和板的正中間,量子點(diǎn)的極化方向沿著垂直于銀平板表面的方向,量子點(diǎn)的半徑為aFig.5.Schematic diagram for the nano-cavity.A silver nano-sphere with radius ra=20 nm is above a silver planar.The gap between them is g.A quantum dot with radius a is placed at the center and is assumed to be polarized normal to the planar.

圖6 球-平板納米腔系統(tǒng)中的歸一化自發(fā)輻射率和能級(jí)移動(dòng) 考慮不同的量子點(diǎn)半徑a,其中a=0對(duì)應(yīng)利用散射格林函數(shù)的結(jié)果,a=0對(duì)應(yīng)利用重整化的散射格林函數(shù)的結(jié)果,(a)和(b)為g=4 nm的結(jié)果,(c)和(d)為g=2 nm的結(jié)果Fig.6.Spontaneous emission rate and the energy level shift of a quantum emitter in plasmonic gap nano-cavity.We consider different sizes of the quantum emitter.a=0 represent the results by the method in Ref.[26].Accordingly,results for a=0 is by our renormalized scattering Green function method.(a)and(b)Arefor large gap g=4 nm,(c)and(d)arefor a small gap g=2 nm.

當(dāng)球和平板間的間隙為g=4 nm時(shí),自發(fā)輻射和能級(jí)移動(dòng)分別見(jiàn)圖6(a)和圖6(b).我們考慮了三種不同的量子點(diǎn)半徑,其中黑色實(shí)線為量子點(diǎn)半徑a=0 nm時(shí)(即未采用重整化,即文獻(xiàn)[26]所示的通過(guò)兩次仿真得到散射格林函數(shù))的結(jié)果,紅色虛線(藍(lán)色空心圓圈)為量子點(diǎn)半徑a=0.5 nm(a=1.5nm)時(shí)的結(jié)果.我們發(fā)現(xiàn):當(dāng)量子點(diǎn)體積較小時(shí)(a=0.5 nm),其結(jié)果與利用嚴(yán)格散射場(chǎng)(a=0 nm時(shí)表示未重整化)得到的結(jié)果符合得很好,此結(jié)論對(duì)于更小的球-平板距離情況下也適用(g=2 nm情況見(jiàn)圖6(c)和圖6(d))；當(dāng)量子點(diǎn)非常大時(shí)(半徑a=1.5 nm),體積對(duì)自發(fā)輻射和能級(jí)移動(dòng)的影響仍然較小(對(duì)比圖6(a)和圖6(b)中藍(lán)色空心圓圈與黑色實(shí)線).這些結(jié)果表明,量子點(diǎn)的體積對(duì)自發(fā)輻射和能級(jí)移動(dòng)的影響較小.另外,類(lèi)似于文獻(xiàn)[26],我們也發(fā)現(xiàn)納米腔的最低階共振模式頻率(偶極模式,圖6(a)和圖6(c)的第一個(gè)峰所對(duì)應(yīng)的頻率)相對(duì)于單球情況發(fā)生了紅移,且空隙距離越小,紅移越明顯,這一現(xiàn)象可通過(guò)模式雜化理論解釋[45].更重要的是,我們發(fā)現(xiàn)了極大的自發(fā)輻射增強(qiáng)和能級(jí)移動(dòng),例如,當(dāng)g=2 nm時(shí),在偶極共振頻率附近,自發(fā)輻射率增強(qiáng)因子高達(dá)Γ/Γ0≈4.3×105,能級(jí)移動(dòng)可達(dá)?≈?320 meV,在高階共振頻率附近(ω=2.995 eV),Γ/Γ0≈ 2.02×106,?≈ 1000 meV此能級(jí)移動(dòng)可與原子的躍遷頻率相比擬,不難想象,由于?∝d2((22)式),進(jìn)一步增大偶極矩的大小,此能級(jí)移動(dòng)能超越原子的躍遷頻率,能使系統(tǒng)進(jìn)入超強(qiáng)耦合區(qū)域,系統(tǒng)中可能會(huì)出現(xiàn)本征能量為負(fù)值的束縛態(tài)[46].

4 結(jié) 論

本文提出了一種利用有限元計(jì)算重整化格林函數(shù)的方法.該方法通過(guò)計(jì)算點(diǎn)電偶極子的輻射場(chǎng),并將其重整化,得到重整化的格林函數(shù),減去均勻空間重整化格林函數(shù)的解析解,可求得重整化的散射格林函數(shù).應(yīng)用該方法到均勻空間中,數(shù)值解與解析解一致,證明該方法能處理均勻空間中格林函數(shù)的重整化問(wèn)題；在納米球系統(tǒng)中,我們發(fā)現(xiàn)重整化散射格林函數(shù)與解析散射格林函數(shù)也一致,證實(shí)該方法能處理散射格林函數(shù)的重整化問(wèn)題.相比于文獻(xiàn)[26]所示的方法,本方法約減小一半的計(jì)算量.應(yīng)用該方法到球-平板納米腔系統(tǒng)中,發(fā)現(xiàn)自發(fā)輻射和能級(jí)移動(dòng)與量子點(diǎn)的體積關(guān)系不大,因此,利用我們的方法,不需要精確的量子點(diǎn)的大小就能獲知系統(tǒng)中自發(fā)輻射和能級(jí)移動(dòng)的特性；另外,在此系統(tǒng)中,我們發(fā)現(xiàn)了極大的自發(fā)輻射率的改變和能級(jí)移動(dòng),例如,當(dāng)球離板很近(g=2 nm)時(shí),在高階共振頻率附近(ω=2.995 eV),自發(fā)輻射增強(qiáng)因子Γ/Γ0≈2.02×106,能級(jí)移動(dòng)?≈1000 meV.