李 靜， 胡 真

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

光子晶體［1］是一種新型的周期性人造材料，它的一個主要特征是具有光子帶隙，頻率落在帶隙中的光不能夠在光子晶體中傳播。光子晶體的這個特征已被用于設計制造各種光子晶體器件，例如波導管彎曲［2-3］、分支［4］、頻率濾波器［5］等。

在實際應用中，需要具有盡可能寬的帶隙光子晶體結構。二維光子晶體由某種材料按照一定周期形式排列構成，如有正方形排列、三角形排列和蜂巢狀排列等。其中蜂巢狀排列的光子晶體由于能夠用以設計出比較寬的完全帶隙結構［6-7］，即光子帶隙在兩種偏振模式下存在重疊，因此在理論和工業(yè)上都被廣泛研究［8-9］。在分析和設計蜂巢狀排列的光子晶體的帶隙結構中，有效的數(shù)值方法是必要的，而對帶隙結構的分析可以轉(zhuǎn)化為求解特征值問題。在標準的模型中［1］，特征值為ω2(ω 是角頻率)，Bloch 波矢是給定的參數(shù)。對于一個非色散介質(zhì)，這是一個線性特征值問題可以通過平面波展開法［10］、有限元法［11-12］、有限差分法［13-14］等進行求解;對于色散介質(zhì)，這是一個非線性特征值問題。當介質(zhì)是非色散的，可以利用柱面波展開法［15-16］等將問題轉(zhuǎn)化成非線性特征值問題;另外基于傳輸矩陣法［17］、散射矩陣法［18］等建立起特征值問題，特征值為Bloch 波矢的一個函數(shù)，頻率為給定的參數(shù)，這時即使介質(zhì)色散，得到的也是一個線性特征值問題。

Yuan 等［19-20］提出利用DtN 映射計算正方形排列和三角形排列的二維光子晶體的帶隙結構，通過建立合適的單元晶格DtN 映射，將帶隙結構問題轉(zhuǎn)化成特征值問題求解，特征值為Bloch 波矢的一個函數(shù)，頻率ω 為給定的參數(shù)。DtN 映射實質(zhì)是把單元晶格邊界上的波動場映射成波動場的法向?qū)?shù)，利用DtN 映射避免了在單元晶格內(nèi)部的離散，從而得到的是較小矩陣的特征值問題。

文中將利用DtN 映射計算蜂巢狀排列的光子晶體的帶隙結構，蜂巢狀光子晶體的單元晶格由3 個正六邊形組成。首先構造出這個特殊單元晶格的DtN 映射，然后將蜂巢狀光子晶體的帶隙結構問題轉(zhuǎn)化成特征值問題求解。特征值問題中的矩陣是比較小的，如果在單元晶格每條邊上離散N 個點，則矩陣的大小為(24N)× (24N)。

1 單元晶格上的DtN 映射

考慮標準的二維蜂巢狀排列的光子晶體結構，如圖1(a)所示。其中，白色圓柱代表空氣洞，背景介質(zhì)為硅板。對于一個折射率函數(shù)為n = n(x，y)的z方向上不變的二維結構，波沿xy 平面?zhèn)鞑サ那闆r下，控制方程是Helmholtz 方程:

其中，k0為自由空間波數(shù);n 為折射率函數(shù)。在E 偏振模式下，u 是電場的z 分量并且μ = 1;在H 偏振模式下，u 是磁場的z 分量并且μ = n2。

圖1 蜂巢結構光子晶體及其單元晶格Fig.1 Two-dimensional honeycomb photonic crystals and the unit lattice

文中方法是基于單元晶格上的DtN 映射，DtN映射是指把單元晶格邊界上的波動場映射成波動場的法向?qū)?shù)，蜂巢狀光子晶體的單元晶格Ω 由3個正六邊形Ω1，Ω2，Ω3組成，如圖1(b)所示。用Γ表示單元晶格Ω 的邊界AB…L(圖1(b)中的外圍)。DtN 映射滿足

其中，u 在Ω 中滿足Helmholtz 方程;v 為邊界Γ 的單位法向矢量。如果用uj1 ≤j ≤12 表示Γ 的12 條邊AB，BC，…，KL，式(2)可以寫成

如果在單元晶格Ω 每條邊離散N 個點，那么DtN 映射Λ 可以用(12N)× (12N)的矩陣近似。單元晶格Ω 由3 個正六邊形Ω1，Ω2，Ω3組成，為了計算單元晶格Ω 的DtN 映射Λ，首先需要計算3 個正六邊形Ω1，Ω2，Ω3的DtN 映射Λ(j)(j =1，2，3)。這3 個正六邊形的DtN 映射的構造過程和文獻［19-21］相同。若Ω1，Ω2和Ω3的邊界分別為Γ1(ABCDOL)，Γ2(KLOHIJ)和Γ3(ODEFGH)，那么Ω1，Ω2，Ω3的DtN 映射Λ(1)，Λ(2)和Λ(3)滿足

為了構造Λ(j)(j = 1，2，3)的近似矩陣，首先寫出Helmholtz 方程式(1)的通解，用極坐標(r，θ)表示:

式中

其中，n1為圓柱內(nèi)部材料的折射率;n2為背景材料的折射率;系數(shù)Aj，Bj可以利用圓柱邊界上的連續(xù)性條件求出。如果在邊界Γj(j =1，2，3)的每條邊上離散N 個點，通解式(5)可以近似寫成

有了Λ(j)(j = 1，2，3)，就可以計算蜂巢狀光子晶體單元晶格DtN 映射Λ。用u13，u14和u15表示邊DO，LO 和HO，則DtN 映射式(4)可以寫成:

將Λ(j)(j = 1，2，3)分成6 ×6 塊，每一塊是一個N ×N 的矩陣，表示為(1 ≤j ≤3，1 ≤s，t ≤6)，并且記

利用式(6)，可以寫出單元晶格Q 邊界上每一條邊上的解的法向?qū)?shù)j = 1，…，12，則有

為了得到單元晶格Ω 的DtN 映射Λ，就要消除式(8)中的q，即u13，u14和u15，將S 分塊成矩陣S11(12N ×12N)和S12(12N ×13N)，則有

其中，q = ［u13，u14，u15］T可以通過式(6)建立方程組消除。如對于u13，由于u13為DO 邊上的解，DO 邊同時屬于正六邊形Ω1和六邊形Ω3，利用Λ(1)和Λ(3)分別寫出?vu13的表達式，則有

同樣，關于u14和u15可以得到另外兩個方程，從而得到方程組:

因此蜂巢狀光子晶體的單元晶格Ω 的DtN 映射Λ 為

其中，Λ 為(12N)× (12N)的矩陣。

2 特征值問題

為了分析二維蜂巢狀光子晶體的帶隙結構，考慮方程式(1)的Bloch 模式解:

其中，(α，β)為Bloch 波矢;Ψ 為周期函數(shù)，與折射率函數(shù)n 遵循同樣的周期。接下來利用DtN 映射Λ建立特征值問題。這里特征值問題的建立過程和文獻［19-20］類似。由二維蜂巢狀光子晶體結構的周期性，利用Bloch 模式解式(13)得到

其中ρα= eiαL/2，。如果把Λ 分塊成12 ×12 的矩陣，每一塊矩陣為N ×N 的，則可以得到如下方程:

將式(14)～式(16)帶入式(17)，消除ui(7 ≤i ≤12)，?vuj(1 ≤j ≤12)，化簡得到

其中，M1，1，M1，2，…，M1，6，…，M6，6與ρα和ρβ相關。

蜂巢狀排列的介質(zhì)柱的布里淵區(qū)如圖2 所示。

對蜂巢狀排列的光子晶體，只需要尋找不可約布里淵區(qū)域邊界上的Bloch 模式解。圖2 中不可約布里淵區(qū)域是三角形陰影部分ΓMK，其中Γ(α = 0，β =0)，M(α = 0，β = 2π/)，K(α = 2π/3L，β =2π/)。從Γ 到M 有ρα= 1，ρβ為特征值;從M 到K 有ρβ= －1，ρα為特征值。對于K 到Γ 的情況，由于光子晶體的對稱性，通過旋轉(zhuǎn)K 到K'(α =4π/4L，β = 0)計算K 到Γ 的情況，此時ρβ= 1，ρα為特征值。

圖2 蜂巢狀排列的介質(zhì)柱的布里淵區(qū)Fig.2 Brillouin zone of a honeycomb lattice of dielectric columns

最后無論在三角形區(qū)域ΓMK 的哪一條邊上求解特征值，都可以寫成

其中，U = (u1，u2，…，u6)T，特 征 值λ 和 矩 陣M0，…，M4定義如下:

1)從Γ 到M，λ = ρβ，0 ＜β ≤2π/，并且

其中 C1= Λ9，3－ Λ4，3－ Λ4，10+ Λ9，10，

C2= Λ9，4－ Λ4，4－ Λ4，9+ Λ9，9，

C3= Λ10，3－ Λ3，3－ Λ3，10+ Λ10，10，

C4= Λ10，4－ Λ3，4－ Λ3，9+ Λ10，9。

2)從M 到K，λ = ρα(0 ＜α ≤2π/3L)，并且M4，M3，M2，M1和M0同樣是由DtN 映射Λ 的分塊矩陣拼出來的大矩陣。

3)從K'到Γ，λ = ρα(0 ＜α ≤4π/3L)，且矩陣M0，M2，M4和情況2 的矩陣M0，M2，M4相同，矩陣M1，M3和情況2 的矩陣M1，M3符號相反。

將式(19)式改寫為下面的線性特征值問題:

其中，I 為單位矩陣;V = ［λ3UT，λ2UT，λ1UT，UT］T。對于單元晶格Ω 每個邊離散N 個點，矩陣Mj(0 ≤j ≤4)都是(6N)× (6N)的矩陣。因此線性特征值問題式(20)的矩陣是(24N)× (24N)的。

3 數(shù)值算例

通過計算兩種不同的二維蜂巢狀排列的光子晶體的帶隙結構，驗證文中建立算法的有效性。

首先考慮標準的二維蜂巢狀光子晶體的帶隙結構(見圖1(a))。在此結構中，背景硅板上空氣洞的大小一致，半徑為a = 0.25L，L 是晶格常數(shù)。文中分別計算兩種偏振模式下的帶隙結構。首先是H 偏振模式(u 表示磁場的z 分量)的情況，控制方程為式(1)，其中μ = n2。背景介質(zhì)硅板的有效折射率為n2= 2.91。此時由于標準情況下空氣洞的大小相同，所以正六邊形Ω2和Ω3完全一樣，因此它們的DtN 映射也是一樣的，即Λ(3)= Λ(2)。在單元晶格Ω每條邊離散7 個點(N = 7)的情況下，特征值(20)的矩陣是(24N)× (24N)的。計算結果如圖3(a)所示。其中，縱坐標是標準化頻率ωL/(2πc)，ω 為角頻率，c 為真空光速;橫坐標為不可約布里淵區(qū)域的邊界。虛線部分是光子帶隙，頻率范圍為［0.523 3，0.544 8］。

然后計算E 偏振模式(u 表示電場的z 分量)的帶隙結構，控制方程為式(1)，其中μ = 1。背景介質(zhì)硅板的有效折射率為n2= 2.52，單元晶格Ω 每條邊同樣離散7 個點。計算結果如圖3(b)所示。

圖3 兩種偏振模式下的光子能帶結構Fig.3 Band structure of two-dimensional honeycomb photonic crystals is shown in Fig.1.1(a)

圖3 (b)中，虛線部分為光子帶隙，頻率范圍為［0.282 4，0.315 5］和［0.500 4，0.541 0］。蜂巢狀光子晶體在E 偏振模式和H 偏振模式下的光子帶隙有重疊(頻率范圍為［0.523 3，0.541 0］)，因此蜂巢狀光子晶體在頻率范圍［0.523 3，0.541 0］里存在完全帶隙。文中計算結果和文獻［7］的結果是吻合的。

再考慮不標準的二維蜂巢狀光子晶體的帶隙結構，這時光子晶體中空氣洞的大小是不一樣的，如圖4 所示，其中較大的空氣洞的半徑為r1=0.38L，較小的空氣洞的半徑為r2= 0.14L。此時構成單元晶格Ω 的3 個正六邊形中，Ω2為包含較小空氣洞的正六邊形，Ω3為包含較大空氣洞的正六邊形，它們的DtN 映射Λ(2)與Λ(3)并不相同，需要分別計算。在單元晶格Ω 每條邊上同樣離散7 個點，計算結果如圖5 所示。圖5(a)為H 偏振模式下的帶隙結構，有效折射率n2= 2.91，光子帶隙的頻率范圍為［0.311 8，0.419 2］;圖5(b)為E 偏振模式下的帶隙結構，有效折射率n2= 0.52，光子帶隙在頻率范圍［0.311 6，0.344 5］內(nèi)。因此蜂巢狀光子晶體在頻率范圍［0.311 8，0.344 5］里存在完全帶隙。文中計算結果和文獻［7］的結果完全一致的。

圖4 不標準的二維蜂巢狀光子晶體的帶隙結構Fig.4 Two-dimensional honeycomb photonic crystals composed of two kinds of air holes

圖5 兩種偏正模式下的光子能帶結構Fig.5 Band structure of two-dimensional honeycomb photonic crystals is shown in Fig.4

4 結 語

主要介紹利用DtN 映射建立起特征值問題計算二維蜂巢狀排列的光子晶體的帶隙結構。由于DtN映射避免了在光子晶體的單元晶格內(nèi)部的離散，所以得到的是較小矩陣的特征值問題，可以快速求解，因此這種方法有利于分析和設計具有更好的帶隙結構的蜂巢狀光子晶體。

［1］Joannopoulos J D，Meade R D，Winn J N.Photonic Crystals:Molding the Flow of Light［M］.Princeton，NJ:Princeton University Press，1995.

［2］Smajic J，Hafner C，Erni D.Design and optimization of an achromatic photonic crystal bend［J］.Optics Express，2003，11(12):1378-1384.

［3］Ikuno H，Naka Y.Finite-difference time domain method applied to photonic crystals［J］.Electromagnetic Theory and Applications for Photonic Crystals，2006(103):401-443.

［4］Koshiba M，Tsuji Y，Hikari M.Time-domain beam propagation method and its application to photonic crystal circuits［J］.Journal of Lightwave Technology，2000，18(1):102.

［5］Ogusu K，Takayama K. Transmission characteristics of photonic crystal waveguides with stubs and their application to optical filters［J］.Optics letters，2007，32(15):2185-2187.

［6］Kalra Y，Sinha R K.Modelling and design of complete photonic band gaps in two-dimensional photonic crystals［J］. Pramana，2008，70(1):153-161.

［7］WEN F，David S，Checoury X，et al.Two-dimensional photonic crystals with large complete photonic band gaps in both TE and TM polarizations［J］.Optics Express，2008，16(16):12278-12289.

［8］Sondergaard T，Broeng J，Bjarklev A，et al. Suppression of spontaneous emission for a two-dimensional honeycomb photonic bandgap structure estimated using a new effective-index model［J］.Quantum Electronics，1998，34(12):2308-2313.

［9］YE J Y，Mizeikis V，XU Y，et al.Fabrication and optical characteristics of silicon-based two-dimensional photonic crystals with honeycomb lattice［J］.Optics Communications，2002，211(1):205-213.

［10］Johnson S，Joannopoulos J. Block-iterative frequency-domain methods for Maxwell＇s equations in a planewave basis［J］. Optics Express，2001，8(3):173-190.

［11］Axmann W，Kuchment P.An efficient finite element method for computing spectra of photonic and acoustic band-gap materials:I.scalar case［J］.Journal of Computational Physics，1999，150(2):468-481.

［12］Dobson D C.An efficient method for band structure calculations in 2D photonic crystals［J］.Journal of Computational Physics，1999，149(2):363-376.

［13］GUO S，WU F，Albin S，et al.Photonic band gap analysis using finite-difference frequency-domain method［J］.Optics Express，2004，12(8):1741-1746.

［14］YU C P，CHANG H C.Applications of the finite difference mode solution method to photonic crystal structures［J］.Optical and Quantum Electronics，2004，36(1):145-163.

［15］Nicorovici N A，McPhedran R C，Botten L C.Photonic band gaps for arrays of perfectly conducting cylinders［J］.Physical Review E，1995，52(1):1135.

［16］Ohtaka K，Ueta T，Amemiya K.Calculation of photonic bands using vector cylindrical waves and reflectivity of light for an array of dielectric rods［J］.Physical Review B，1998，57(4):2550.

［17］Pendry J B.Calculating photonic band structure［J］.Journal of Physics:Condensed Matter，1996，8(9):1085-1108.

［18］Botten L C，Nicorovici N A，McPhedran R C，et al. Photonic band structure calculations using scattering matrices［J］. Physical Review E，2001，64(4):046603.

［19］YUAN J，LU Y Y.Photonic bandgap calculations with Dirichlet-to-Neumann maps［J］.JOSA A，2006，23(12):3217-3222.

［20］YUAN J，LU Y Y. Computing photonic band structures by Dirichlet-to-Neumann maps:the triangular lattice［J］. Optics Communications，2007，273(1):114-120.

［21］HUANG Y，LU Y Y. Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps［J］. Journal of Lightwave Technology，2006，24(9):3448.