王夢夢，胡 真

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

0 引言

光子晶體[1]是由不同折射率的介質在空間中周期性排列形成的人造材料.它以光的波長尺度為周期，具有獨特的光學性質，如存在光子帶隙[2]、負折射率[3,4]等，從而能夠被制造成各種光子晶體元件，在集成光路[5]、超微激光器[6]、光伏[7]等領域發(fā)揮了重要的作用.光子晶體通常有兩種制造方法，一種是在空氣中排列介質柱(介質柱光子晶體)[8]，另外一種是在薄板上挖取空氣洞(薄板光子晶體)[9,10].由于在空氣中排列介質柱的制造方法需要增加襯底，而在薄板上挖取空氣洞可以直接進行，所以實際應用更廣泛的是第二種薄板光子晶體[1].

近年來,人們發(fā)現在薄板光子晶體與另外一種介質的交界面上存在著表面模式[11]，即在邊界上存在一種可以傳播的電磁波.它可以被用于設計制造出性能更加優(yōu)異的光子晶體元件，如S S Xiao等[12]在2005年以及T W Lu等[13]在2009年分別設計出了具有更高Q因子的微腔結構；Z Y Zhang等[14]在2007年設計并制造出了一種消光比更大、系統(tǒng)Q因子更高的濾波器結構；H Kurt等[15]在2012年設計出了一種可以更大范圍調控帶寬和群組折射率的慢光結構；M N Erim等[16]在2013年設計出了一種在TE和TM模式下都具有高傳輸率的波導結構；Z Huang等[17]在2015年設計出了一種能夠把光更有效的輸入光子晶體波導中的耦合結構.

在分析與優(yōu)化設計薄板光子晶體表面模式的過程中，需要有效的數值方法.通常用于分析波導管與介質柱光子晶體表面模式的數值方法，都可以擴展后用于分析薄板光子晶體的表面模式，常用的方法有平面波展開法、有限差分法、有限元方法和時域有限差分法等.如H Chen等[18]在2006年利用平面波展開法分析了三角形周期排列薄板光子晶體的表面模式；M Che等[19]在2008年利用平面波傳輸矩陣方法分析了正方形和三角形周期排列介質柱光子晶體表面模式；Y C Hsu等[20]在2010年利用有限元方法分析了正方形周期排列介質柱光子晶體的表面模式；M N Erim等[16]在2013年采用平面波展開法和時域有限差分法分析了三角形周期排列薄板光子晶體的表面模式.

但是上述這些傳統(tǒng)的方法各有不足.時域有限差分法更適于計算脈沖的傳播等在時間域上更方便表示的問題，而表面模式分析描述的是波數與頻率之間的關系，因此更適合利用頻率域方法來進行計算.傳統(tǒng)的頻率域方法有：平面波展開法、有限元方法和有限差分法，但是這些方法對于非色散介質制作的薄板光子晶體，建立的特征值問題是線性的，而對于色散介質制作的薄板光子晶體，建立的特征值問題是非線性的，較難計算.并且上述時域、頻率域方法均需要在整個計算區(qū)域上進行離散，未知數個數比較多，得到的特征值問題所涉及的矩陣是比較大的.

近年來Y Huang等[21]提出了一種新的頻率域方法，叫做Dirichlet-to-Neumann (DtN)映射方法，這種方法通過利用單元晶格的DtN映射，并針對不同元件的特點，發(fā)展出一系列有效模擬各種光子晶體元件的數值方法.這些方法被廣泛用來分析光子元件，Huang等[21,22]2006年分析了光子晶體多層結構的傳播率與反射率問題，他們也在2007年分析了正方形周期排列介質柱光子晶體波導管；J Yuan等[23,24]2006年分析了正方形周期排列介質柱光子晶體的光子帶隙，他們也在2006年分析了三角形周期排列介質柱光子晶體的光子帶隙；Z Hu等[25]2015年分析了薄板光子晶體異質結構；王依翠等[26]2016年分析了正方形周期排列介質柱光子晶體的表面模式.單元晶格的DtN映射是一種算子，它將邊界上的波動場映射成波動場的法向導數.通過利用單元晶格的DtN映射，只需要在單元晶格的邊界上進行離散，不需要在單元晶格的內部進行離散，從而減少了未知數的個數，而且對于相同的單元晶格只需要構造一次DtN映射.光子晶體元件具有大量相同的單元晶格，所以特別適合利用DtN映射方法.

本文擴展了DtN映射方法，將其用于分析三角形周期排列薄板光子晶體的表面模式.在計算過程中，首先對標準的正六邊形單元晶格和非標準的五邊形單元晶格分別建立DtN映射，然后構造出超級晶格的DtN映射，最后在超級晶格的邊界上建立較小矩陣的特征值問題來進行求解.

與平面波展開法、有限元方法和有限差分法等傳統(tǒng)頻率域方法相比，本文方法主要具有以下兩個方面的優(yōu)點.首先，本文建立的特征值問題所涉及的矩陣比較小，從而能夠加快求解速度.這是因為上述傳統(tǒng)頻率域方法都需要在整個計算區(qū)域上進行離散，而本文的方法只需要在計算區(qū)域的邊界上離散，從而大大減少了未知數的個數.對于薄板光子晶體元件的優(yōu)化設計等需要大量重復運算的問題來說，本文的算法具有更大的優(yōu)勢.

其次，本文的方法無論對色散介質還是非色散介質制作的薄板光子晶體，建立的特征值問題都是線性的，容易求解.而上述傳統(tǒng)頻率域方法對于色散介質制作的薄板光子晶體，建立的特征值問題是非線性的，難以求解.這是因為上述傳統(tǒng)頻率域方法構造的特征值問題是以波數β為給定參數，以ω2為特征值，ω是角頻率，而本文方法構造的特征值問題是以角頻率ω為給定參數，以波數β的函數eiβL為特征值.

1 單元晶格的DtN映射

本文研究的二維薄板光子晶體結構如圖1(a)所示，灰色是薄板介質，白色的圓表示空氣洞，整個結構在垂直于紙面的z方向上沒有變化.在實際制造中薄板介質是有限厚度的三維結構，但是三維結構的模擬計算非常困難，而二維結構的模擬計算也能夠在很大程度上反應出實際三維元件的性能，所以往往利用有效折射率方法[27]將三維問題轉化為二維問題進行模擬分析.在圖1(a)所示的結構中，薄板介質的有效折射率為n2，空氣洞的半徑為a，L是晶格常數，它的邊界面位于薄板光子晶體與空氣之間，邊界面與第一排空氣洞邊緣之間的距離為d.圖1(b)所示為選取的超級晶格，即圖1(a)所示的結構的一個周期.

(a)三角形周期排列薄板光子晶體

(b)超級晶格圖1 三角形周期排列薄板光子晶體及其超級晶格

圖1(b)所示的超級晶格由標準的正六邊形單元晶格和非標準五邊形單元晶格組成，這些單元晶格如圖2所示.

(a)標準無缺陷正六邊形單元晶格 (b)標準有缺陷正六邊形單元晶格

(c)非標準無缺陷五邊形單元晶格 (d)非標準有缺陷五邊形單元晶格圖2 超級晶格(圖1(b))中包含的四種單元晶格

標準正六邊形單元晶格的DtN映射構造過程在Y Huang等[21,22]和J Yuan等[23,24]的論文中有詳細的介紹，非標準五邊形單元晶格的DtN映射構造過程與它們類似.下面簡要描述如何構造這些單元晶格的DtN映射.

在本文的問題中，波沿xy平面?zhèn)鞑ィ贖-偏振模式下，控制方程是Helmholtz方程：

(1)

式(1)中：u是磁場的z分量，n=n(x,y)為折射率函數，k0=ω/c為自由空間波數，ω為角頻率，c是真空下的光速.圖2(a)、2(c)分別表示無缺陷的正六邊形單元晶格及五邊形單元晶格，它們包含空氣洞，空氣洞的半徑為a,單元晶格中的點到圓心的距離為r，當ra時n=n2，圖2(b)、2(d)分別表示有缺陷的正六邊形單元晶格和五邊形單元晶格，它們不包含空氣洞，內外介質相同，n=n1=n2=1.

在圖2所示的單元晶格Ω(i)(i=1,2,3,4)中，Helmholtz方程(1)的通解u可以表示成特解的線性組合，表達形式如下：

(2)

式(2)中：φj(r)eijθ是Helmholtz方程在Ω(i)中的一個特解，特解可以選用柱面波解，cj是未知系數，φj與貝塞爾函數Jj、Yj有關，可以表示為

上式中：Aj,Bj是系數，可以由邊界面r=a上的連續(xù)性條件求出;(r,θ)是(x,y)在極坐標系下的坐標.

單元晶格Ω(i)的DtN映射Λ(i)實際上是一種算子，將單元晶格Ω(i)邊界Γ(i)上的波動場u映射成 的法向導數.即

(3)

式(3)中：u是Helmholtz方程在Ω(i)中的解；v為單位法向量.

在計算過程中，需要對每一個不同的單元晶格Ω(i)(i=1,2,3,4)分別構造出它的DtN映射Λ(1).對于正六邊形單元晶格Ω(2)(圖2(a))、(圖2(b))，如果每條邊上離散N個點，那么總共可以將方程(1)的通解u用s=6N個特解的線性組合來表示.對于五邊形單元晶格Ω(3)(圖2(c))、Ω(4)(圖2(b))來說，由于邊的長度并不是完全相同的，所以在選取離散點的時候，盡可能使得每條邊上離散點的間距相近.假設在AB、AE兩條邊上分別離散N個點，BC、DE兩條邊上分別離散N1個點，CD邊上離散N2個點，那么總共可以把方程(1)的通解u用s=2N+2N1+N2個特解的線性組合來表示，若s為偶數，則有：

(4)

式(4)中：(xk,yk)是邊界Γ(i)上的離散點，k=1,2,…,s.

對式(4)兩邊同時求法向導數，得到

(5)

將式(4)和式(5)寫成如下的矩陣形式

(6)

式(6)中：u=[u(x1,y1),u(x2,y2),…,u(xs,ys)]T，C=[c1,c2,…,cs]T是未知數向量，在式(6)的兩個方程組中消去C，可以得到DtN映射的近似矩陣形式

Λ(i)=BA-1

(7)

式(7)中：Λ(i)(i=1,2,3,4)是s×s的矩陣.由于相同單元晶格的DtN映射只需要構造一次，且超級晶格(如圖1(b)所示)只包含四種不同的單元晶格，所以總共只需構造四次單元晶格的DtN映射.

2 超級晶格的DtN映射及特征值問題

在構造完成圖2所示四種單元晶格的DtN映射之后，開始構造超級晶格的DtN映射.為了描述的更清楚，將圖1(b)所示的超級晶格單獨提出并放大，如圖3所示.

圖3 包含m個單元晶格的超級晶格

圖3所示的超級晶格由m個單元晶格構成，邊QmPm表示薄板光子晶體與空氣的邊界面，邊界面的左邊由有缺陷的正六邊形和五邊形單元晶格構成，右邊由無缺陷的正六邊形和五邊形單元晶格構成，整個結構是左右對稱的.

超級晶格的DtN映射把超級晶格上側邊界和下側邊界上的波動場映射成邊界上波動場的法向導數，即滿足：

(8)

式(8)中：u=[u0,1,u0,2,…,u0,2m-2]T，u1=[u1,1,u1,2,…,u1,2m-2]T,u0,1,u0,2,…,u0,2m-1，表示邊P1P2,P2P3,…,P2m-2,P2m-1上的波動場，u1,1,u1,2,…,u1,2m-1,表示邊Q1Q2,Q2Q3,…,Q2m-2,Q2m-1上的波動場.由于表面模式在遠離邊界面QmPm的位置上是衰減的，所以在m足夠大的時候，即在超級晶格左右邊界Q1O1、O1P1、Q2m-1O2、O2P2m-1遠離QmPm時，可以在其上取零邊界條件.

構造超級晶格DtN映射M的基本思想是，利用單元晶格的DtN映射在超級晶格內部的公共邊上建立起線性方程組，然后消去公共邊上的波動場，從而得到超級晶格的DtN映射.記公共邊Q3P2,Q4P5,…,Q2m-3P2m-2上的波動場為w1,w2,…,wm-1，下面以公共邊QmPm為例來說明如何建立線性方程，這條邊上的波動場為wm/2.

如圖3所示，公共邊QmPm位于無缺陷五邊形單元晶格Ω(3)(圖2(c))和有缺陷五邊形單元晶格Ω(4)(圖2(d))的交界處，在上一節(jié)中已經建立了這兩個單元晶格的DtN映射，分別為Λ(3)和Λ(4).

由于公共邊QmPm屬于單元晶格Ω(3)，且單元晶格Ω(3)有五條邊，所以將DtN映射Λ(3)根據每條邊的離散點數分成5×5的分塊形式，如下式所示：

(9)

取公式(9)中第三行的展開式，得到公共邊QmPm上波動場wm/2法向導數的表達式：

(10)

同樣由于公共邊QmPm也屬于單元晶格Ω(4)，將單元晶格Ω(4)的DtN映射Λ(4)也分成5×5的分塊，如下式所示:

(11)

取公式(11)中第一行的展開式，得到公共邊QmPm上波動場wm/2法向導數的另一種表達式：

(12)

公共邊QmPm是位于兩種介質的交界面處，左側介質的折射率為n1，右側介質的折射率為n2，連續(xù)性條件為：

(13)

結合式(10)、式(12)和式(13)，化簡整理得到線性方程：

(14)

對圖3中超級晶格內部的每一條公共邊都可以按照上述過程建立線性方程，從而得到線性方程組：

(15)

式(15)中：w=[w1,w2,…,wm-1]T.

對于位于超級晶格上側邊界上的邊P1P2,P2P3,…,P2m-2P2m-1和位于下側邊界上的邊Q1Q2,Q2Q3,…,Q2m-2Q2m-1來說，由于它們都屬于某一種單元晶格，所以可以利用該單元晶格的DtN映射，將它們波動場的法向導數用相應單元晶格邊界上的波動場的線性組合來表示，形如：

(16)

利用式(15)消去公共邊上的波動場w，再將其代入式(16)得到:

(17)

即超級晶格的DtN映射M為：

(18)

式(18)中：M是(4mN-4N1)×(4mN-4N1)維的矩陣.

得到超級晶格的DtN映射之后，接下來就可以在超級晶格的上下邊界上建立起特征值問題來求解表面模式.一個表面模式就是Helmholtz方程(1)的形如

u(x,y)=Φ(x,y)eiβx

(19)

的特解.其中β是布洛赫波數，即表面模式的傳播常數，Φ是x方向上的周期函數，在y軸正負方向上衰減到零.

利用擬周期邊界條件

其中λ=eiβL，將式(8)化成如下所示的特征值問題：

(20)

式(20)中：λ是特征值,I是恒等算子，Mij(i,j=1,2)是超級晶格DtN映射M的2×2分塊矩陣，維度為(2mN-2N1)×(2mN-2N1).

3 數值算例

本文計算三角形周期排列薄板光子晶體的表面模式，結構圖如圖1(a)所示.H Chen等[18]在2006年用平面波展開法研究了這個結構.在這個結構中，它的邊界面位于薄板光子晶體與空氣層之間，邊界面與第一排空氣洞邊緣之間的距離d=0.7L，空氣洞的半徑為a=0.3L，L為晶格常數，薄板材料的有效折射率為n2=3.5.在標準正六邊形單元晶格的每條邊選取N=5個離散點，在非標準五邊形單元晶格AB、AE邊上選取N=5個離散點，BC、DE邊上選取N1=11個離散點，CD邊上選取N2=8個離散點，選取超級晶格的層數m=14，得到精確到小數點后第四位的計算結果，如圖4(a)所示.

(a)色散曲線

(b)圖4(a)中點P1、P2對應表面模式的磁場分布圖圖4 d=0.7L時圖1(a)所示三角形周期排列薄板光子晶體與空氣層之間邊界面上表面模式的計算結果

接下來，本文將利用收斂性測試說明在計算時如何選擇超級晶格層數m及離散點數N，使得計算結果精確到小數點后第四位.仍然以圖4(a)中P1為例，取標準化頻率ωL/(2πc)=0.2505，首先固定超級晶格層數m=14，改變標準正六邊形單元晶格每條邊上的離散點數N，五邊形單元晶格離散點數相應改變，計算結果如表1所示；然后固定離散點數N=5，改變超級晶格層數m，得到的計算結果如表2所示.如果只要求計算結果收斂到小數點后第四位，從表1、表2可以看出βL/(2π)=0.359 8是滿足條件的具有四位有效數字的計算結果，為了達到這一標準并且盡可能使得計算所涉及的未知數個數比較少，從表1中可以看出選取m=14,N=5是比較合適的.

表1 固定m改變N

表2 固定N改變m

在第二個算例中，改變邊界面與第一排空氣洞邊緣之間的距離，使得d=0.1L，空氣洞的半徑及背景材料有效折射率不變，同樣選擇超級晶格層數m=14，標準單元晶格每條邊離散點數N=5，得到精確到小數點后第四位的計算結果，如圖5(a)所示.本文結果仍然與H Chen等[18]文章中的結果一致.在色散曲線圖(圖5(a))中選擇點P3(0.469 2,0.207)，畫出它對應的表面模式的磁場分布圖，如圖5(b)所示.

(a)色散曲線

(b)圖5(a)中點P3對應表面模式的磁場分布圖圖5 d=0.1L時圖1(a)所示三角形周期排列薄板 光子晶體與空氣層之間邊界面上表面模式的計算結果

應用本文發(fā)展的數值方法，在薄板光子晶體與PMC材料的邊界面上還發(fā)現了一種新的表面模式.結構圖如圖6(a)所示，紅色直線代表邊界面，在邊界面的上側是PMC材料，邊界面的下側是三角形周期排列薄板光子晶體.邊界面與第一排空氣洞邊緣之間的距離d=0.6L，空氣洞的半徑及背景材料的有效折射率與前兩個算例一致.此時超級晶格只包含標準無缺陷的正六邊形單元晶格Ω(1)(如圖2(a)所示)和非標準無缺陷五邊形單元晶格Ω(3)(如圖2(c)所示)這兩種單元晶格，所以在實際計算時只需要求出這兩種單元晶格的DtN映射Λ(1)與Λ(3)即可.

同樣選擇超級晶格層數m=14，標準單元晶格每條邊的離散點N=5,得到精確到小數點后第四位的計算結果，色散曲線圖如圖6(b)所示.在色散曲線上選取點P4(0,348 1,0.273 5)、P5(0.356 1,0.235 5)，分別畫出它們對應的表面模式的磁場分布圖，如圖6(c)所示，與圖4(b)相似，圖中顏色的亮度表示|Hz(x,y)|/Hz_max值的大小.在圖像中可以看出顏色較亮區(qū)域集中在邊界面上，在垂直于邊界面的薄板光子晶體內部顏色逐漸變暗.由于磁場無法穿透PMC材料，邊界面上側的磁場值全部為零，所以沒有畫出邊界面上側部分的磁場分布圖.從圖6(c)可以看出,表面模式的能量集中在邊界面上，在薄板光子晶體內部快速衰減，只能沿著薄板光子晶體的邊界面進行傳播.

(a)三角形周期排列薄板光子晶體與PMC材料的邊界面

(b)色散曲線

(c)圖6(b)中點P4、P5對應表面模式的磁場分布圖圖6 三角形周期排列薄板光子晶體與PMC材料邊界面上表面模式的計算結果

在分析計算圖6(a)所示的三角形周期排列薄板光子晶體與PMC材料的邊界面結構的過程中，可以發(fā)現將直線型邊界面變成鋸齒狀邊界面，仍然存在著表面模式.新的結構圖如圖7(a)所示，鋸齒狀邊界面沿著正六邊形單元晶格的邊界面分布，邊界面上側為PMC材料，下側為三角形周期排列薄板光子晶體.此時超級晶格中只包含一種標準正六邊形單元晶格Ω(1)(圖2(a))，所以只需要計算Λ(1)即可.空氣洞的半徑和薄板材料的有效折射率仍然與前三個算例一致，同樣選擇超級晶格層數m=14，標準正六邊形單元晶格每條邊的離散點N=5,得到精確到小數點后第四位的計算結果，色散曲線圖如圖7(b)所示.取圖7(b)中色散曲線上的一點P6=0.459,0.202 2，畫出它對應的表面模式的磁場分布圖，如圖7(c)所示.

(a)三角形周期排列薄板光子晶體與PMC材料的鋸齒狀邊界面

(b)色散曲線

(c)圖7(b)中點P6對應表面模式的磁場分布圖圖7 三角形周期排列薄板光子晶體與PMC材料的鋸齒狀邊界面上表面模式的計算結果

4 結論

本文擴展了DtN映射方法，將其用于計算三角形周期排列薄板光子晶體的表面模式.以標準與非標準單元晶格的DtN映射為基礎，構造出超級晶格的DtN映射，最后在超級晶格的兩條邊界上建立矩陣較小的線性特征值問題來進行求解.應用這種方法在薄板光子晶體與PMC材料的交界面上還發(fā)現了一種新的表面模式.

在本文研究的薄板光子晶體結構中，挖取的空氣洞截面都是標準圓形，因此可以在單元晶格中用Helmholtz方程的柱面波特解的線性組合來表示通解，從而構造出單元晶格DtN映射.如果空氣洞的截面不是標準圓形，比如橢圓、三角形等等，這種構造DtN映射的方法將不再適用.接下來將進一步研究在空氣洞截面不是標準圓形的薄板光子晶體結構中如何建立有效的數值方法分析計算表面模式，此時需要利用邊界積分方程方法來構造單元晶格的DtN映射.