任瑩蓉， 胡真

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

光子晶體[1]是一種以光波的波長尺度為周期的新型人造周期性材料,其最顯著的特征是具有光子帶隙.頻率落在光子帶隙內的光波不能在光子晶體中傳播,這使得光子晶體具有了控制光的能力,被廣泛應用于高輻射超小型激光器[2]、高光電轉化率光伏電池[3]和量子信息元件[4]等諸多領域.

光子晶體表面模式[1]是指在光子晶體的邊界面上所形成的表面波,它能夠在邊界面上傳播,在垂直于邊界面的方向上指數衰減.對于傳統的各向同性光子晶體,表面模式總是成對出現,即在同一頻率下,既存在正向傳播的表面模式,也存在反向傳播的表面模式,故光波在邊界面上既能正向傳播也能反向傳播.旋磁光子晶體由在各個方向上磁導率不同的旋磁各向異性材料制成.人們發(fā)現在外加磁場的作用下,旋磁光子晶體因其破壞了時間反演對稱性,從而能夠存在單向表面模式[5-6],即在一個頻率下,只存在沿著一個方向傳播的表面模式,故光波在邊界面上只能沿著一個方向傳播,當邊界面位于旋磁光子晶體與各向同性光子晶體或兩種不同旋磁光子晶體的交界處時,即使結構存在缺陷,光波也不能反向散射,傳輸率仍然能夠達到100%,整個結構是魯棒的.基于這些優(yōu)點,Chen等人設計了無損耗單向波導[6]; Zhang等人設計了高精度環(huán)形諧振器[7]; Ren等人設計了無反射三端口通道降頻濾波器[8].

分析與計算光子晶體表面模式需要有效的數值方法.現有的數值方法包括有限元法[7-8]、平面波展開法[9]、傳輸矩陣法[10]、散射算子方法[11]等；有限元法和平面波展開法均需要離散計算區(qū)域內部,因此涉及的矩陣不容易求解,且對于色散介質,建立的是復雜的非線性特征值問題;傳輸矩陣法存在數值不穩(wěn)定性；散射算子方法的計算則比較困難.

近年來,Lu等[12]提出了一種新的頻率域數值方法,即Dirichlet-to-Neumann(DtN)映射方法.DtN映射本質上是將單元晶格邊界上的波動場映射成波動場的法向導數,從而能夠避免對計算區(qū)域內部進行離散,極大地減少了計算量.DtN映射方法已經被用于分析各類光子晶體結構[13-18],但從未被用于旋磁光子晶體單向表面模式的分析與計算.

本文將DtN映射方法用于計算旋磁光子晶體的單向表面模式，首先挑選出合適的超級晶格,即結構的一個周期,然后構造出超級晶格的DtN映射,最后在超級晶格的邊界上建立起線性特征值問題來進行求解.

1 控制方程

為方便討論僅考慮z方向上不變的二維結構[12-18]，圖1所示的結構為一種由旋磁光子晶體(右半部分)和各向同性光子晶體(左半部分)組成的二維結構,單向表面模式存在于這兩種光子晶體的交界面處.

圖1 由旋磁光子晶體(藍)和各向同性光子晶體(綠)組成的二維結構

(1)

其中u是電場的z分量,n(x,y)是折射率函數,k0是自由空間波數.

對于各向同性光子晶體,磁導率μ=μ0是一個常數,E-極化下的控制方程為Helmholtz方程

(2)

2 單元晶格的DtN映射

圖1的二維結構中,旋磁光子晶體和各向同性光子晶體均由正方形周期排列的介質柱組成,其中各向同性光子晶體旋轉了45°后放置在旋磁光子晶體左側.圖1中綠色實線和紅色實線圍成的區(qū)域即為挑選出的超級晶格,它包含三種不同的單元晶格,如圖2所示.

為了構造超級晶格的DtN映射,首先需要構造出單元晶格的DtN映射.單元晶格的DtN映射Λ將邊界Γ上的波動場u映射成邊界上波動場u的法向導數,即:

(3)

其中ν是邊界Γ上的單位法向量.

(a) 旋磁各向異性單元晶格(正方形) (b) 各向同性單元晶格(五邊形) (c) 各向同性單元晶格(正方形)

對于只包含一個半徑為a的介質柱的單元晶格,方程(1)和(2)的通解u可由特解的線性組合表出:

(4)

(5)

(6)

如果在單元晶格的邊界Γ上離散p個點,則可以用p個特解的線性組合來近似方程(1)和(2)的通解u,然后利用(4)中的通解u及其法向導數?νu的表達式消去未知系數cj.最終,可以用一個p×p矩陣來近似DtN映射Λ.

3 超級晶格的DtN映射和特征值問題

在得到三種單元晶格的DtN映射后,可以構造出超級晶格的DtN映射M:

(7)

構造的基本思想[16-17]是利用超級晶格內部的公共邊上波動場法向導數的連續(xù)性,在公共邊上建立起線性方程組,從而消去公共邊上的波動場,最終得到超級晶格的DtN映射M.

為了便于敘述超級晶格的DtN映射M的構造過程,需要引入一些記號:如圖3所示,將2m+1個單元晶格(從右到左)依次記為Ω(i)(i=1,2,...,2m+1),用u0,i(i=1,2,...,2m+1)表示底邊界Q0Q1,…,Q2mQ2m+1的波動場u,用u1,i(i=1,2,...,2m+1)表示上邊界P0P1,…,P2mP2m+1的波動場u,用wi(i=1,2,...,2m)表示公共邊(虛線)的波動場u,則u0=[u0,1,...,u0,2m+1]T,u1=[u1,1,...,u1,2m+1]T,w=[w1,...,w2m]T.

圖3 由2m+1個單元晶格組成的超級晶格(m是偶數)

在超級晶格左右兩邊紅色邊界(邊SP2m+1,SQ2m+1,P0Q0)處采用零邊界條件,即波動場u=0.這是因為問題的計算是在光子帶隙內進行,頻率落在光子帶隙內的光波不能在光子晶體內傳播,故當左右紅色邊界遠離兩種光子晶體的交界處時,波動場快速衰減到0.

以圖4所示的單元晶格Ω(m)和Ω(m+1)之間的公共邊PmQm為例在超級晶格內的公共邊上建立線性方程組.

圖4 單元晶格Ω(m)與Ω(m+1)的公共邊

對單元晶格Ω(m)和Ω(m+1)的DtN映射Λ(m)和Λ(m+1)分別進行4×4和5×5分塊,則PmQm上波動場wm的法向導數?νwm既可以表示成Ω(m)的每條邊上波動場wm-1,u0,m,wm,u1,m的線性組合,也可以表示成Ω(m+1)的每條邊上波動場wm,u0,m+1,wm+1,u1,m+2,u1,m+1的線性組合,再根據?νwm的連續(xù)性能夠得到除去PmQm以外的其余7條邊的線性方程.對于超級晶格內其它公共邊,可以用同樣的方法構造線性方程.最終,所有公共邊上的波動場w可以由u0,u1線性表出:

(8)

其中P0是m(N+N2)×m(N+N2)矩陣，P1是m(N+N2)×2(m(N+N2)+N1)矩陣。

另一方面,超級晶格邊界上波動場的法向導數?νu0,?νu1有如下表達式:

(9)

其中Q0是2(m(N+N2)+N1)×2(m(N+N2)+N1)矩陣,Q1是2(m(N+N2)+N1)×m(N+N2)矩陣.將(8)式代入(9)式中可求得超級晶格的DtN映射M為:

(10)

且M是2(m(N+N2)+N1)×2(m(N+N2)+N1)矩陣.

對于圖1所示結構,其表面模式是控制方程的一種形如下式的特解:

u(x,y)=e-iβxΨ(x,y)

(11)

其中β是實的傳播常數,Ψ(x,y)是在x方向上以L為周期的周期函數,且當y→±時,u衰減到0.

將(7)式中的矩陣M進行2×2分塊,再利用擬周期性邊界條件u1=e-iβLu0,?νu1=e-iβL?νu0,最終得到線性特征值問題

(12)

通過求解得到特征值e-iβL，從而得到計算結構的單向表面模式.

4 數值算例

以存在著單向表面模式的兩種不同的邊界面結構為例來驗證方法的有效性,第一種邊界面位于旋磁光子晶體和各向同性光子晶體的交界處,第二種邊界面位于兩種不同的旋磁光子晶體的交界處.

(a) 色散曲線 (b) 點A1對應的電場分布圖

圖5(a)中的紅色實線上的點A1(0.377 2,0.560 4)表示在歸一化頻率ωL/(2πc)=0.560 4時,存在傳播常數βL/(2π)=0.377 2的單向表面模式.點A1的電場分布如圖5(b)所示,從圖5(b)中可以看出能量集中在兩種光子晶體的交界處,在垂直于交界的方向上快速衰減到0.

圖6 由兩種不同的蜂巢狀旋磁光子晶體組成的結構

第二個算例結構如圖6所示,由Chen等人[6]提出.該結構由兩種不同的蜂巢狀旋磁光子晶體組成.蜂巢狀旋磁光子晶體(晶格常數L=15 mm)均由旋磁各向異性材料制成的介質柱在空氣中按蜂巢狀周期排列組成,介質柱的折射率均為(15.26),右半部分介質柱的半徑分別為a1=3 mm,a2=2.5 mm,左半部分介質柱的半徑分別為a3=3.6 mm,a4=1.9 mm.當外加磁場H0=500 Oe時,μr=0.78μ0,μk=0.93μ0.取m=11,N=7[16-17],計算結果如圖7(a)所示,其中單向表面模式存在于兩種光子晶體的交界處,其頻率范圍為[0.261 5,0.322 6].計算結果同文獻[6]中的結果一致.圖7(b)給出了點B1(0.386 9,0.288 4)對應的單向表面模式的電場分布圖.

這兩個算例表明,利用DtN映射方法能夠得到與有限元方法[5-6]完全一致的計算結果,但DtN映射方法的優(yōu)點在于它只需要離散計算區(qū)域的邊界,避免了對計算區(qū)域內部進行離散,從而極大地減少了未知數的個數,使得求解特征值問題時涉及的矩陣規(guī)模非常小,提高了計算速度.

(a) 色散曲線 (b) 點B1對應的電場分布圖

5 結語

將DtN映射方法用于計算旋磁光子晶體單向表面模式.只需要離散超級晶格的邊界,在邊界上建立起一個以Bloch波矢的函數為特征值的線性特征值問題,該特征值問題涉及的矩陣規(guī)模小,求解速度快,且無論是對于色散介質還是非色散介質,該特征值問題都是線性的,因此DtN方法非常適合用于旋磁光子晶體的分析與計算.