基于Dirichlet-to-Neumann映射計算波導中缺陷激發(fā)的連續(xù)譜束縛態(tài)*

陸明鑫， 胡真

(河海大學 理學院，江蘇 南京211100)

1 引言

在正能量的連續(xù)域內(nèi)可以存在孤立解，即局域模式的能量落在連續(xù)體內(nèi)，這種嵌模不與擴展波耦合，具有無限的壽命，就像零寬度的諧振[1]，量子力學中將這種狀態(tài)稱為連續(xù)域束縛態(tài)(Bound states in the continuum)，簡稱BIC.最近，人們認識到BIC本質(zhì)上是一種波動現(xiàn)象，在光學上，由于其無限高質(zhì)量因子(Q因子)和對系統(tǒng)參數(shù)變化的魯棒性，使其在生物傳感[2-4]和開關(guān)[5]以及非線性光學[6]等方面上有很多應(yīng)用.目前，BIC在光學領(lǐng)域中的研究已經(jīng)取得了很大進展[7-8].

在眾多存在BIC的材料中，光子晶體[9]因其設(shè)計簡單和易于制作而備受關(guān)注.當在光子晶體波導中引入兩個對稱的點缺陷時，由于缺陷之間的相互作用，可以用兩能級模型[10]的近似結(jié)構(gòu)來描述它們的耦合機制.通過調(diào)節(jié)缺陷的結(jié)構(gòu)參數(shù)，使強烈的共振干擾將共振態(tài)與連續(xù)體的耦合完全抑制，這時便激發(fā)出了BIC.

分析與計算光子晶體波導中的BIC需要高效的數(shù)值方法.求解光子晶體波導中由缺陷激發(fā)的BIC，等價于求解缺陷的實共振頻率，而這是一個腔膜問題，因此可以作為麥克斯韋方程的特征值問題來求解.傳統(tǒng)的數(shù)值方法，如有限差分法[11-12]、有限元法和平面波展開法[13]等，求解時都會產(chǎn)生大矩陣特征值問題，因此求解比較困難.此外，對于色散介質(zhì)，矩陣項依賴于頻率，特征值問題會變得非線性，將更加難以求解.

2006年Huang等[14]提出了一種特殊的頻域方法——Dirichlet-to-Neumann (DtN)映射方法，可以用來精確模擬二維光子晶體器件.DtN映射方法通過將邊界上的波場映射為邊界上對應(yīng)的法向?qū)?shù)，以避免對計算區(qū)域內(nèi)部的計算，將原來的特征值問題變成一個小矩陣線性特征值問題，其中特征值是Bloch波矢量，ω是一個給定參數(shù)，這種方法降低了計算的工作量和難度.目前，DtN映射方法已經(jīng)被廣泛應(yīng)用于分析各類光子晶體結(jié)構(gòu)[15-18].

本文將DtN映射方法應(yīng)用于二維光子晶體波導中由缺陷激發(fā)的BIC的計算與分析，其主要思想便是利用嚴格的邊界條件截斷波導管，從而得到一個有限的計算區(qū)域，再通過單元晶格的DtN映射進一步將計算域縮小到在單元晶格上的所有邊上，最終在邊界上建立起求解小矩陣的特征值問題.基于該方法，便能確定激發(fā)波導中BIC的點缺陷的半徑和介電常數(shù)范圍.

2 DtN映射方法

在二維光子晶體波導中引入兩個對稱的點缺陷為諧振腔，可以激發(fā)出BIC，以點缺陷平行于波導排列的結(jié)構(gòu)(圖1(a))為例，進行DtN映射方法構(gòu)建.該結(jié)構(gòu)是一個由空氣中的介質(zhì)棒按正方形晶格排列的光子晶體，其晶格常數(shù)為L，介質(zhì)棒半徑r，介電常數(shù)為ε，背景空氣折射率為n.其中黑色填充圓表示波導中具有可變介電常數(shù)εd的附加局部缺陷棒，可以引起左右波導中傳播模式的散射，從而激發(fā)出BIC.

(a)完整二維光子晶體波導 (b)截斷區(qū)域

對于在z方向上不變，且波在xy方向上傳播的二維問題，E-偏振下，控制方程為

(1)

其中，k0=w/c是自由空間波數(shù)；w是角頻率；c是真空中的光速；n(x,y)是折射率函數(shù)；u是在E-偏振下方向z上的電場分量(唯一非零分量).

考慮上述結(jié)構(gòu)，缺陷晶格放置在了光子晶體波導附近.因此，為了計算缺陷模式，對附近波導建立精確的入射和出射的邊界條件尤其重要.本文通過使用嚴格的邊界條件來截斷PhC波導，將計算域截斷為圖1(b)所示矩形S，假設(shè)計算域由0

(2)

(3)

其中L+和L-是兩個作用于y上的算子，它們依賴于給定的頻率，并且當對y進行離散時，便可以用矩陣逼近.

(a)單元晶格 (b)超級晶格

確定了截斷PhC波導的邊界條件后，將計算域劃分為正方形晶格并計算DtN映射[14].對于如圖2(a)所示，若每邊離散N個點，DtN映射Λ可近似為一個4N×4N的矩陣，即

(4)

一般來說，邊界上每條邊上的點在選取上都會避開頂點進行均勻選擇.

在建立了單元晶格的DtN映射后，對于計算域中內(nèi)部的邊，可以利用兩個相鄰單元晶格的DtN映射來匹配?vu；在計算域邊界的邊上，則使用一個單元的DtN映射和嚴格的邊界條件即(2)，(3)來匹配?vu.最終在單元晶格的邊上，得到了一個這樣的線性方程組

A(ω)U=0，

(5)

其中，U是計算域S中所有邊上u的列向量；A是依賴于ω的矩陣.以圖1為例，若計算域是p行q列的單元晶格組成的，則總邊數(shù)m=2pq+p+q，那么U是一個mN×1的列向量，A是一個mN×mN的方陣.

結(jié)構(gòu)的共振頻率ω需要通過求解方程(5)的非零解來得到，當A是一個奇異矩陣時，可以使用非線性方程

λ1(A(ω,r,ε))=0

(6)

來迭代求解根ω.其中，λ1是A(ω)的特征值的絕對值中最小的特征值.當固定半徑r時，對于公式(6)，需要求解出對應(yīng)于實頻率ω的一組解(ω,ε).這是因為BIC是有著完美的局域性且沒有能量泄露，而復頻率ω表示為ω=ω0-iγ，其中虛部則表示能量泄露，故判斷BIC的條件是γ=0.

對于固定半徑r，以缺陷晶格的介電常數(shù)ε為變量，當找到對應(yīng)于實頻率模式的共振模，便找到了BIC.因此，當以半徑為變量，對范圍內(nèi)的每個半徑值，只要能找到對應(yīng)于實頻率ω的解，便有了一組能夠激發(fā)BIC的點缺陷參數(shù)對(r,ε)，如此就確定了能夠激發(fā)波導中BIC的點缺陷的半徑和介電常數(shù)范圍.

3 邊界條件

以圖1(b)為例，討論x=3L右側(cè)的半無限光子晶體波導的邊界條件，對于波導的一個周期xj

(7)

其中，uj和?xuj是xj處的電場值及其法向?qū)?shù)；M是一個2mN×2mN的矩陣.

基于算子M，可以計算波導的Bloch模，對于周期性的結(jié)構(gòu)，波導中的波場是Bloch模的疊加，一個Bloch模就是方程(1)的一個特解，即

φ(x,y)=Φ(x,y)eiβx，

(8)

其中Φ在x上關(guān)于L具有周期性，β是常數(shù).因此對于μ=eiβL來說，有

φj+1=μφj,?xφj+1=μ?xφj，

(9)

因此Bloch模可以通過特征值方程

(10)

求解；其中，特征值μ=eiβL，φ和?xφ為xj處的值，I是單位矩陣.通過利用Bloch模，可以將光子晶體中微腔的波場展開為

(11)

當頻率ω為復頻率時，條件|μ|≤1并不總是有效的，因為此時原本在實頻率下隨x的增加而衰減的模也會隨著x的增加而增長，故需要進行修正.由于BIC的共振頻率的虛部通常是非常小的，對于這樣的復頻率ω，它的Bloch模的特征值與它對應(yīng)實部的頻率的Bloch模的特征值非常接近，因此修正的方法是計算并篩選復頻率ω實部所對應(yīng)的Bloch 模來作為復頻率所對應(yīng)的Bloch模[18].

基于方程(11)的Bloch模展開，有

(12)

其中，φk=Φk(xj,y)，μk=eiβkL，因此可以定義算子T使得

Tφk=μkφk,k=1,2,….

由于T是線性的，則有uj+1=Tuj，且由(7)中給出的DtN算子M可以得到

?xuj=L+uj，其中L+=M11+M12T.

(13)

通過同樣的方法，也可以得到左邊界上的L-算子.綜上，得到了終止PhC波導的嚴格邊界條件.

4 數(shù)值算例

本節(jié)計算了幾種結(jié)構(gòu)算例[20]中能夠激發(fā)BIC的點缺陷的結(jié)構(gòu)參數(shù)范圍.

第一個算例如圖1所示，計算區(qū)域為3×15個正方形單元晶格.利用本文方法，以點缺陷半徑rd為變量，得出當rd在[0.056,0.5]范圍內(nèi)變化時，都可以找到一個介電常數(shù)εd，使得結(jié)構(gòu)中出現(xiàn)BIC，如圖3所示.取參考文獻中的半徑rd=0.18L，得到εd=3.035 8，wL/2πc=0.374 3，這與文獻[19]中的結(jié)果一致，該點電場圖如圖3中右側(cè)所示.

圖3 圖1所示結(jié)構(gòu)εd與rd的關(guān)系圖和rd=0.18L時電場圖

第二個算例如圖4(a)所示，計算區(qū)域為5×15個正方形單元晶格.利用本文方法，以點缺陷半徑rd為變量，得出當rd在[0.044,0.5]和[0.127,0.5]范圍內(nèi)變化時，都可以找到一個介電常數(shù)εd，使得結(jié)構(gòu)中分別出現(xiàn)了奇BIC和偶BIC，如圖4(b)和4(c)所示.取參考文獻中的半徑rd=0.18L，分別在εd=1.990 9，wL/2πc=0.374 1激發(fā)出了奇BIC，以及εd=5.615 2，wL/2πc=0.327 3時激發(fā)出了偶BIC，結(jié)果與參考文獻[19]中的結(jié)果一致，其中這兩個BIC的電場圖如圖4(b)和4(c)中插圖所示.

(a) 缺陷平行波導結(jié)構(gòu)圖 (b) 奇BIC的εd與rd關(guān)系圖 (c) 偶BIC的εd與rd關(guān)系圖

算例三是一個缺陷相對于坐標x→-x，y→-y反對稱排列的結(jié)構(gòu)，如圖5(a)所示.這里將計算區(qū)域截斷為3×15個正方形單元晶格.利用本文方法，以點缺陷半徑rd為變量，得出當rd在[0.051,0.5]范圍內(nèi)變化時，都可以找到一個介電常數(shù)εd，使得結(jié)構(gòu)中出現(xiàn)BIC，如圖5(a)所示.取參考文獻中的半徑rd=0.18L，得到εd=2.517 9，wL/2πc=0.374 2，這與參考文獻[19]中的結(jié)果一致，其中這點電場圖如圖5(b)中插圖所示.

(a)缺陷反對稱排列結(jié)構(gòu)圖 (b) εd與rd的關(guān)系圖

5 結(jié)語

將DtN映射方法應(yīng)用于二維光子晶體波導中BIC的計算與分析,確定了在光子晶體波導中能夠激發(fā)波導中BIC的點缺陷的結(jié)構(gòu)參數(shù).本文的方法使得不需要在整個平面上進行計算，而只需要在截斷區(qū)域內(nèi)計算，而即便在較小的計算域中，也只需要在單元晶格的邊界上進行離散，而不需要考慮晶格內(nèi)部，從而最終在單元晶格邊界上建立起小規(guī)模的特征值問題，降低了計算量和難度.本文只是對理想的二維結(jié)構(gòu)進行計算分析，下一步考慮將這種方法擴展到三維的結(jié)構(gòu)上去，研究三維波導管中的連續(xù)譜束縛態(tài).