王雄輝，胡 真

(河海大學(xué) 理學(xué)院， 南京 211100)

光子晶體[1]是由介電常數(shù)不同的介質(zhì)在空間中按周期排列制成的新型人造周期性材料，它以光的波長(zhǎng)尺度為周期，其最重要的性質(zhì)是具有光子帶隙。頻率落在帶隙內(nèi)的光不能在光子晶體內(nèi)傳播，因此光子晶體具有獨(dú)特的控制光的能力。光子晶體波導(dǎo)管是通過(guò)在完整的光子晶體中引入線缺陷破壞其內(nèi)部周期性結(jié)構(gòu)形成的，頻率落在光子帶隙內(nèi)的光可以在線缺陷波導(dǎo)管中傳播，在遠(yuǎn)離線缺陷時(shí)會(huì)迅速衰減，因此線缺陷起到導(dǎo)波的效果。光子晶體波導(dǎo)管是連接集成光路各元件之間的橋梁，是光信號(hào)傳輸?shù)闹饕d體。利用光子晶體和光子晶體波導(dǎo)管可以制造各種光學(xué)元件，如高性能發(fā)光二極管[2]、超低功率激光器[3]、光子晶體諧振腔[4]、超高Q濾光器[5]等元件，這些元件具有集成度高、性能優(yōu)異等優(yōu)點(diǎn)，在大規(guī)模集成光路和全光網(wǎng)絡(luò)中具有廣闊的應(yīng)用前景。

近年來(lái)，由旋磁各向異性材料制成的旋磁光子晶體得到了廣泛關(guān)注，其中的旋磁各向異性材料是指材料在各個(gè)方向上具有不同的磁導(dǎo)率。旋磁光子晶體在外加磁場(chǎng)的作用下打破了時(shí)間反演對(duì)稱性，獲得了各向同性光子晶體所不具有的可控性，比如通過(guò)在旋磁光子晶體中引入線缺陷形成旋磁光子晶體波導(dǎo)管，波導(dǎo)管中的光波在外加磁場(chǎng)的作用下可以實(shí)現(xiàn)單向傳輸。利用旋磁光子晶體和旋磁光子晶體波導(dǎo)管可以制造出一系列可控光學(xué)元件。例如，Nazari[6]在2016年制造了一種可單向傳輸?shù)某o湊型光隔離器；Esmaieli等[7]在2011年制造了一種μ可切換的緊湊型對(duì)稱功率分配器；Zhang等[8]在2013年制造了一種具有沿特定端口定向傳輸?shù)乃亩丝谛殴庾泳w環(huán)形器；吳昌義[9]在2016年設(shè)計(jì)制造了具有高對(duì)比度、大帶寬、小尺寸的光子晶體磁控二選一光路選通開關(guān)。這些旋磁光子晶體元件具有高透射率、大帶寬、小尺寸、可控、高消光比和低插入損耗等優(yōu)勢(shì)，在光通信和光信號(hào)處理等領(lǐng)域中扮演著重要角色。

在模擬和分析旋磁光子晶體元件的過(guò)程中，有效的數(shù)值方法非常重要。目前已有的數(shù)值方法存在各種局限性。例如，時(shí)間域下模擬的時(shí)域有限差分方法[10]對(duì)傳輸和反射光譜這類在頻率域下更容易模擬的問題，使用時(shí)并不方便。標(biāo)準(zhǔn)頻域下模擬的數(shù)值方法包括有限元法[11]、平面波展開法[12]、傳輸矩陣法[13]等，其中有限元法和平面波展開法需要離散晶格內(nèi)部，建立邊值問題涉及的線性方程組規(guī)模大、未知數(shù)多且計(jì)算量大；傳輸矩陣法存在數(shù)值不穩(wěn)定的缺陷。近年來(lái)，Huang等[14]提出了能有效模擬光子晶體的Dirichlet-to-Neumann(DtN)映射方法，Yuan等[15-16]利用DtN映射方法計(jì)算二維光子晶體帶隙；Ye等[17]利用DtN映射方法計(jì)算二維三角形排列和蜂巢狀排列光子晶體的等頻線；Wang等[18]利用DtN映射方法計(jì)算二維蜂巢狀排列光子晶體波導(dǎo)；Hu等[19]利用DtN映射方法計(jì)算二維各向同性光子晶體元件的傳輸和反射光譜，分析了有損周期介質(zhì)結(jié)構(gòu)中連續(xù)譜附近具有復(fù)頻率的束縛態(tài)[20-21]。單元晶格的DtN映射將單元晶格邊界上的波動(dòng)場(chǎng)映射成邊界上的法向?qū)?shù)，避免在單元晶格內(nèi)部進(jìn)行離散，減少了計(jì)算量。對(duì)于相同的單元晶格，只需要構(gòu)造一次DtN映射。旋磁光子晶體元件有大量相同的單元晶格，故在計(jì)算旋磁光子晶體元件過(guò)程中有很大優(yōu)勢(shì)。

本研究擴(kuò)展了DtN映射方法，將其用于計(jì)算二維旋磁光子晶體元件的傳輸和反射光譜。方法步驟為：① 計(jì)算結(jié)構(gòu)中所有不同類型單元晶格的DtN映射，單元晶格的DtN映射可以用一個(gè)小矩陣近似表出；② 利用計(jì)算出的單元晶格的DtN映射構(gòu)造出超級(jí)晶格的DtN映射，在超級(jí)晶格的邊界上建立特征值問題，求解出波導(dǎo)管內(nèi)的Bloch模；③ 基于求解出的Bloch模建立邊界條件截?cái)喙庾泳w波導(dǎo)管，截取有限的計(jì)算域；④ 在截取的計(jì)算域中所有單元晶格的每條邊上建立一個(gè)線性方程，最終得到一個(gè)包含所有邊上波動(dòng)場(chǎng)的線性方程組。為了建立線性方程組，可分別在計(jì)算域的內(nèi)部和邊界上建立線性方程。由于計(jì)算域內(nèi)部的每條邊都位于2個(gè)相鄰單元晶格的交界面上，故利用2個(gè)單元晶格的DtN映射和這條邊上波動(dòng)場(chǎng)法向?qū)?shù)的連續(xù)性建立線性方程。對(duì)于計(jì)算域邊界上的邊，利用單元晶格的DtN映射、截?cái)嗖▽?dǎo)管的邊界條件和這條邊上波動(dòng)場(chǎng)法向?qū)?shù)的連續(xù)性建立線性方程。

1 控制方程

研究對(duì)象為由無(wú)限長(zhǎng)且平行的介質(zhì)柱按正方形周期排列組成的二維旋磁光子晶體元件。如圖1所示，二維旋磁光子晶體環(huán)形器即二維旋磁光子晶體元件中的一種，由Zhang等[8]提出。圖1中紅色部分表示旋磁光子晶體，藍(lán)色部分表示各向同性光子晶體，其中L是晶格常數(shù)，以此結(jié)構(gòu)為例進(jìn)行說(shuō)明。為了對(duì)結(jié)構(gòu)進(jìn)行數(shù)值分析，需要建立嚴(yán)格的邊界條件截?cái)嗖▽?dǎo)管，得到有限的計(jì)算域。圖1中黑色實(shí)線圍成的區(qū)域是該結(jié)構(gòu)的計(jì)算域，為建立邊界條件，需要在超級(jí)晶格的邊界上建立特征值問題求解出波導(dǎo)管內(nèi)的Bloch模，其中由綠色實(shí)線和紅色實(shí)線圍成的區(qū)域是超級(jí)晶格,代表波導(dǎo)管的1個(gè)周期。

對(duì)于旋磁光子晶體,磁導(dǎo)率μ是一個(gè)張量：

(1)

由于旋磁光子晶體和各向同性光子晶體在TM模式下的控制方程相同，不存在特殊的性質(zhì)，故僅在TE模式下進(jìn)行研究。對(duì)于一個(gè)結(jié)構(gòu)在z方向不變、波在xy平面?zhèn)鬏數(shù)亩S問題，它的控制方程為：

(2)

其中：u是電場(chǎng)的z分量；k0=w/c是自由空間波數(shù)，c是真空中的波數(shù)，w是角頻率；n(x,y)是折射率函數(shù)。

對(duì)于各向同性光子晶體，磁導(dǎo)率μ為常數(shù)。此時(shí)，TE模式下的Helmholtz控制方程為

(3)

2 單元晶格的DtN映射

圖1所示的二維結(jié)構(gòu)中包含3種不同類型的正方形單元晶格，分別是旋磁各向異性單元晶格、各向同性單元晶格和缺陷單元晶格。為了對(duì)該結(jié)構(gòu)進(jìn)行有效數(shù)值計(jì)算，首先需要構(gòu)造出這3種不同類型單元晶格的DtN映射。單元晶格Ω的DtN映射Λ本質(zhì)上是將單元晶格邊界Γ上的波動(dòng)場(chǎng)u映射成邊界上波動(dòng)場(chǎng)u的法向?qū)?shù)，即

(4)

其中：u是控制方程(2)和(3)的解；Γ是Ω的邊界；υ表示邊界Γ上的單位法向量。

旋磁光子晶體和各向同性光子晶體的區(qū)別在于構(gòu)成材料磁導(dǎo)率的不同，各向同性光子晶體的磁導(dǎo)率是一個(gè)常數(shù)，而旋磁光子晶體的磁導(dǎo)率是一個(gè)張量。旋磁光子晶體在外加磁場(chǎng)的作用下具有各向同性光子晶體不具有的可控性。各向同性單元晶格DtN映射構(gòu)造過(guò)程可以視為旋磁各向異性單元晶格DtN映射構(gòu)造過(guò)程的一個(gè)特例。因此，接下來(lái)以旋磁各向異性單元晶格(如圖2所示)為例來(lái)說(shuō)明單元晶格DtN映射的構(gòu)造過(guò)程，其中介質(zhì)柱半徑為a。

圖2 旋磁各向異性單元晶格示意圖

單元晶格的DtN映射構(gòu)造思想是將控制方程的通解寫成一組特解的線性組合。對(duì)于只包含1個(gè)介質(zhì)柱的旋磁各向異性單元晶格Ω，控制方程(2)和(3)的通解u(x,y)可近似表示成特解的線性組合形式：

(5)

其中：點(diǎn)(x,y)對(duì)應(yīng)的極坐標(biāo)是(r,θ)；φj(r)eijθ是控制方程(2)和(3)的特解；φj(r)用Bessel函數(shù)Jj和Yj表示：

(6)

Jj(k0n3a)=AjJj(k0n2a)+BjYj(k0n2a)

(7)

(8)

Λ=BA-1

(9)

假設(shè)在單元晶格Ω的每條邊上離散N個(gè)點(diǎn)，則DtN映射Λ可以近似為一個(gè)4N×4N的矩陣。構(gòu)造完只包含1個(gè)介質(zhì)柱的旋磁各向異性單元晶格Ω的DtN映射后，簡(jiǎn)要介紹只包含1個(gè)介質(zhì)柱的各向同性單元晶格和不包含介質(zhì)柱的缺陷單元格晶格的DtN映射構(gòu)造過(guò)程。對(duì)于只包含1個(gè)介質(zhì)柱的各向同性單元晶格，其DtN映射構(gòu)造過(guò)程是旋磁各向異性單元晶格DtN映射構(gòu)造過(guò)程的一個(gè)特例，此時(shí)，P=-1,Q=0,n3=n1，重復(fù)式(6)—(9)完成構(gòu)造。對(duì)于不包含介質(zhì)柱的缺陷單元晶格，P=-1,Q=0，其DtN映射構(gòu)造過(guò)程和各向同性單元晶格的DtN映射構(gòu)造過(guò)程類似，唯一的區(qū)別在于特解的表達(dá)式φj(r)=Jj(k0nr)，此時(shí)n=n2=1。

3 超級(jí)晶格的DtN映射及特征值問題

構(gòu)造完單元晶格的DtN映射后，需要構(gòu)造出超級(jí)晶格的DtN映射來(lái)建立特征值問題，求解出波導(dǎo)管內(nèi)的Bloch模。超級(jí)晶格的DtN映射M滿足：

(10)

超級(jí)晶格的DtN映射M的構(gòu)造思想是利用單元晶格的DtN映射和超級(jí)晶格內(nèi)部所有公共邊上波動(dòng)場(chǎng)法向?qū)?shù)連續(xù)性，在公共邊上建立線性方程組消去公共邊上的波動(dòng)場(chǎng)，從而得到超級(jí)晶格的DtN映射M。

圖3所示的是圖1中由綠色實(shí)線和紅色實(shí)線圍成的超級(jí)晶格，其中x軸在垂直方向，該超級(jí)晶格由m-1個(gè)各向同性單元晶格和1個(gè)缺陷單元晶格組成，圖3中綠色實(shí)線Γ0和Γ1分別為超級(jí)晶格的下邊界和上邊界。在超級(jí)晶格的左、右兩條紅色邊界(S0P0和SmPm)上采用簡(jiǎn)單的零Dirichlet邊界條件，即波動(dòng)場(chǎng)u=0。這是因?yàn)轭l率是在帶隙范圍內(nèi)選取的，該頻率的光波不能在光子晶體中傳播，所以在遠(yuǎn)離波導(dǎo)管處波動(dòng)場(chǎng)迅速衰減到0。

圖3 由m個(gè)單元晶格組成的超級(jí)晶格(m是奇數(shù))示意圖

由于超級(jí)晶格內(nèi)部的每條公共邊都位于相鄰2個(gè)單元晶格的交界面上，故利用2個(gè)單元晶格的DtN映射以及這條邊上波動(dòng)場(chǎng)法向?qū)?shù)的連續(xù)性建立方程，由此得出m-1條內(nèi)部邊上的線性方程，最終得到1個(gè)線性方程組：

(11)

其中：C0是(m-1)N×(m-1)N矩陣；C1是(m-1)N×2mN矩陣；V1是公共邊上的波動(dòng)場(chǎng)；u0和u1分別是超級(jí)晶格下邊界Γ0和上邊界Γ1的波動(dòng)場(chǎng)。

另一方面，可以利用每個(gè)正方形單元晶格的DtN映射表示超級(jí)晶格的下邊界Γ0和上邊界Γ1波動(dòng)場(chǎng)的法向?qū)?shù)：

(12)

聯(lián)立式(11)和式(12)消去V1，得到超級(jí)晶格的DtN映射M表達(dá)式為：

(13)

其中：M是2mN×2mN矩陣；D0是2mN×2mN矩陣；D1是2mN×(m-1)N矩陣。

構(gòu)造完超級(jí)晶格DtN映射之后，需要在超級(jí)晶格2條邊界上建立線性特征值問題求解Bloch模。1個(gè)Bloch模就是Helmholtz控制方程(3)的1個(gè)特解，且滿足：

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

(14)

其中：β是傳播常數(shù)；φ(x,y)是在x軸上以晶格常數(shù)L為周期的周期函數(shù)。利用擬周期條件u1=μu0,?υu(píng)1=μ?υu(píng)0將超級(jí)晶格的DtN映射M分成2×2塊的分塊矩陣，Bloch模由特征值方程(15)解出。

(15)

其中：μ=eiβL是特征值，且該特征值是Bloch波矢函數(shù)；I是恒等算子。

4 建立線性方程組

圖1中的波導(dǎo)管在x和y方向是無(wú)限延伸的，它們起著端口的作用，使光可以進(jìn)入或離開該結(jié)構(gòu)。在實(shí)際數(shù)值模擬中，需要在xy平面截?cái)嗖▽?dǎo)管并使用合適的邊界條件，因此接下來(lái)考慮在光子帶隙范圍內(nèi)給定頻率下的邊值問題。下面介紹如何在切斷1根波導(dǎo)管的邊界上建立邊界條件。

為簡(jiǎn)化敘述，假設(shè)計(jì)算域由0≤x≤mL,0≤y≤mL給出，其中L是晶格常數(shù)，入射波從計(jì)算域左側(cè)的波導(dǎo)管輸入。下面以計(jì)算域的左側(cè)邊界(x=0,0≤y≤mL)為例來(lái)說(shuō)明如何建立截?cái)喙庾泳w波導(dǎo)管的邊界條件。在波導(dǎo)管的1個(gè)周期(即0≤x≤L)，y方向上包含m個(gè)晶格常數(shù)。利用求解出的Bloch模，將波導(dǎo)管中波動(dòng)場(chǎng)展開為

(16)

已知Bloch模是成對(duì)出現(xiàn)的，即對(duì)于傳播常數(shù)為β的Bloch模，也存在另一個(gè)-β的傳播常數(shù)。在|μ|≤1中選擇合適的β，當(dāng)|μ|=1時(shí),選擇實(shí)的β使得Bloch模在x方向上有正的平均功率通量，平均功率通量可通過(guò)Poynting向量計(jì)算；當(dāng)|μ|<1時(shí)，對(duì)應(yīng)隨著x增大指數(shù)衰減的Bloch模，確保Bloch模式展開式(16)只包含向x=+∞傳播的模和隨著x增加而衰減的模。波導(dǎo)管內(nèi)的波動(dòng)場(chǎng)分為向前傳播的波動(dòng)場(chǎng)和向后傳播的波動(dòng)場(chǎng)，即u=u++u-。其中

(17)

(18)

同樣地，得到x=0處u-的邊界條件為：

(19)

根據(jù)u=u++u-和式(18)，式(19)消去u-得到x=0處邊界條件為：

(20)

其中u+是入射波。對(duì)于計(jì)算域的另外3條邊界，由于這3條邊截?cái)嗟牟▽?dǎo)管內(nèi)只有輸出波，故在構(gòu)造其邊界條件時(shí)，不需考慮入射波情況，3條邊的邊界條件為：

(21)

(22)

(23)

對(duì)于圖1所示的計(jì)算域，需要寫出計(jì)算域中每個(gè)單元晶格邊界上波動(dòng)場(chǎng)u的線性方程。為建立線性方程組，分別在計(jì)算域的內(nèi)部和邊界上建立線性方程。對(duì)于計(jì)算域內(nèi)部的每條邊，由于這條邊位于2個(gè)相鄰單元晶格的交界面上，這時(shí)利用2個(gè)單元晶格的DtN映射以及這條邊上波動(dòng)場(chǎng)法向?qū)?shù)的連續(xù)性建立方程。對(duì)于計(jì)算域邊界上的邊，利用單元晶格的DtN映射、截?cái)嗖▽?dǎo)管的邊界條件和這條邊上波動(dòng)場(chǎng)法向?qū)?shù)的連續(xù)性建立另一個(gè)方程。最終得到包含所有邊上波動(dòng)場(chǎng)的線性方程組：

AX=b

(24)

其中:X為截?cái)嘤騼?nèi)所有邊上波動(dòng)場(chǎng)的列向量。若單元晶格每條邊離散N個(gè)點(diǎn)，則系數(shù)矩陣A的大小為2(m+1)mN×2(m+1)mN。內(nèi)部邊的每個(gè)方程只和2個(gè)相鄰單元晶格其余的少量邊有關(guān)，因此A是一個(gè)稀疏矩陣。最后求出傳輸和反射光譜。

5 數(shù)值算例

以四端口旋磁光子晶體環(huán)形器和波導(dǎo)交叉結(jié)構(gòu)為例來(lái)驗(yàn)證算法有效性。所有結(jié)構(gòu)均由旋磁光子晶體和各向同性光子晶體組成，且這些光子晶體均由相應(yīng)材料制成的介質(zhì)柱在空氣中按照正方形點(diǎn)陣排列構(gòu)成。選擇從左側(cè)波導(dǎo)管輸入的傳輸模為入射波，并分別計(jì)算左側(cè)端口的反射率和右側(cè)、上側(cè)、下側(cè)端口的傳輸率。

以表1和表2為例來(lái)說(shuō)明超級(jí)晶格層數(shù)m和單元晶格每條邊上離散點(diǎn)數(shù)N的選取過(guò)程。選取圖4中的點(diǎn)A1對(duì)應(yīng)的頻率wL/(2πc)=0.365 8來(lái)進(jìn)行說(shuō)明：固定超級(jí)晶格層數(shù)m，改變單元晶格每條邊上離散點(diǎn)數(shù)N；固定單元晶格每條邊離散點(diǎn)數(shù)N，改變超級(jí)晶格層數(shù)m。

圖4 算例1結(jié)構(gòu)傳輸和反射光譜

表1 固定超級(jí)晶格總數(shù)m改變單元晶格邊界上的點(diǎn)的個(gè)數(shù)N

表2 固定單元晶格邊界上的點(diǎn)的個(gè)數(shù)N改變超級(jí)晶格總數(shù)m

從表1和表2可以看出，當(dāng)m=15,N=11時(shí)，計(jì)算結(jié)果具有4位有效數(shù)字。

頻率為wL/(2πc)=0.357 1、wL/(2πc)=0.361 6和wL/(2πc)=0.372 8的電場(chǎng)圖見圖5。 3個(gè)頻率分別對(duì)應(yīng)能量從端口A、B和C輸出的最大透射率。從圖5(a)可以看出，當(dāng)頻率為wL/(2πc)=0.357 1時(shí)，大部分能量從A端口定向輸出。由于單向波導(dǎo)的存在，能量從左側(cè)D端口輸入后，在第1個(gè)交叉口處只能向其中1個(gè)方向轉(zhuǎn)向，此頻率下在A端口定向輸出。圖5(b)給出了頻率為wL/(2πc)=0.361 6時(shí)的電場(chǎng)圖，能量主要從B端口輸出。同時(shí)給出了頻率為wL/(2πc)=0.372 8時(shí)的電場(chǎng)圖，如圖5(c)所示。可以看出，大部分能量從C端口定向輸出。

算例2的結(jié)構(gòu)如圖6(a)所示，為波導(dǎo)交叉結(jié)構(gòu)[8]。該結(jié)構(gòu)中有4根YIG棒位于十字交叉口處，YIG棒半徑a2=0.19L，其余半徑a1=0.2L。為了得到4位有效數(shù)字，選取包含15×15個(gè)單元晶格的計(jì)算域。該計(jì)算域包含480條邊，取N=11，得到如圖6(b)所示的傳輸和反射譜，計(jì)算結(jié)果與文獻(xiàn)[8]中結(jié)果一致。從圖6(b)可以看出，當(dāng)頻率為wL/(2πc)=0.366 0時(shí)，A端口的透射率高達(dá)0.97，該頻率下的電場(chǎng)圖如圖6(c)所示?？梢钥闯觯?dāng)能量從左側(cè)D端口輸入時(shí)，大部分能量從A端口輸出。

圖5 不同頻率電場(chǎng)圖

圖6 算例2結(jié)構(gòu)和相關(guān)圖譜

6 結(jié)論

擴(kuò)展了DtN映射的方法，將其用于計(jì)算二維旋磁光子晶體元件的傳輸和反射光譜。DtN映射方法只需在計(jì)算域邊界上進(jìn)行離散，減少了未知數(shù)，最終得到包含計(jì)算域中所有邊上波動(dòng)場(chǎng)的線性方程組。該線性方程組的系數(shù)矩陣是稀疏的，系數(shù)矩陣規(guī)模小，計(jì)算速度快，有利于計(jì)算二維旋磁光子晶體元件的傳輸和反射光譜。研究對(duì)象均為介質(zhì)柱截面為圓形的二維旋磁光子晶體結(jié)構(gòu)，方便進(jìn)一步研究介質(zhì)柱截面非圓形(如橢圓形、三角形和六邊形等)的旋磁光子晶體結(jié)構(gòu)。對(duì)于介質(zhì)柱截面為圓形的光子晶體結(jié)構(gòu)，采用柱面波展開法來(lái)構(gòu)造單元晶格的DtN映射，但對(duì)于介質(zhì)柱截面非圓形的光子晶體結(jié)構(gòu)，該方法不再適用，此時(shí)單元晶格的DtN映射可通過(guò)邊界積分方程方法來(lái)重新構(gòu)造。由于主要研究z方向上不變的二維結(jié)構(gòu)，因此下一步將考慮用垂直方向模展開法來(lái)計(jì)算三維結(jié)構(gòu)。