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

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

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

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

1 引 言

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

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

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

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

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

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

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

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

其中,ε0為真空介電常數(shù),光子并矢格林函數(shù)G(r,r′,ω)滿足[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)可寫為

為求解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í)部對應(yīng)著系統(tǒng)能級,虛部對應(yīng)著衰減快慢.將電場的表達(dá)式((2)式)代入(1)式,利用(5)式,則能移算符矩陣元為:R±aa(ω) ≡ limη→0+?a|R(z=ω±iη)|a?= ?(r0,ω)±Γ(r0,ω)/2.經(jīng)過進(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ā)輻射率和能級移動.本文稱Γ(r；ω)和?(r；ω)為自發(fā)輻射率和能級移動.

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

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

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

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

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

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

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

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

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

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

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

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

則(13)式,即重整化的電場也可寫為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

依有限元理論,對于任意試探函數(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)做展開,如:

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

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

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

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

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

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

表1 真空中點(diǎn)電偶極子的輻射場: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ω.

接下來,展示我們的算法在重整化散射格林函數(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)電偶極子的散射場(解析散射格林函數(shù)見文獻(xiàn)[26,27]),紅色空心圓圈描述的是本文提出的有限元方法得到的點(diǎn)電偶極子的重整化散射場.其求解步驟為:首先求點(diǎn)電偶極子的輻射場E(r,ω),利用(10)式將其重整化,得到Ereg(r0；ω),再利用(17)式得到重整化的格林函數(shù)Greg(r0,r0；ω),隨后利用(18)式(其中,

圖4 納米球系統(tǒng)中的散射場 黑色實(shí)線描述解析的點(diǎn)電偶極子的散射場,可通過解析散射格林函數(shù)或者利用有結(jié)構(gòu)和無結(jié)構(gòu)情況下兩次有限元仿真的差值得到[26]；紅色空心圓圈代表本文的一次有限元仿真得到的點(diǎn)電偶極子的重整化散射場,其中,a=1 nm；可以看到重整化散射場與量子點(diǎn)中心處的散射場符合得很好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.

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

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

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

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

4 結(jié) 論

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