蔣 丹，胡 真

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

光子晶體是一種具有周期性的新型人造材料［1］，其最重要的性質(zhì)是存在光子帶隙［2］。頻率落在光子帶隙的光波不能在光子晶體中傳播，從而使得光子晶體具有控制光的能力，被用來設計和制造各種光子晶體元件，在量子信息、超小型激光器、非線性光學等領域發(fā)揮了重要作用。

制造光子晶體的材料按性質(zhì)不同可分為各向同性和各向異性。各向同性材料的性質(zhì)不會隨外界因素變化而變化，由它們制造的光子晶體元件的性能是確定的，不具有可調(diào)控性。各向異性材料［3］不同于各向同性材料，其性質(zhì)受外界因素的影響會發(fā)生變化。例如，對各向異性材料施以電流或調(diào)整溫度等，會引起各向異性材料性質(zhì)的改變，從而引起由各向異性材料制造的光子晶體結(jié)構性能的改變。由此可以通過施加外界影響來實現(xiàn)對各向異性光子晶體元件的調(diào)控。

液晶［4］是一種常用的各向異性材料，它的性質(zhì)隨外界條件如電場、溫度［5］的變化而變化。近年來，涌現(xiàn)出了大量基于液晶制作的各向異性光子晶體元件，如光調(diào)制器、具有記憶效應的元件等。

在模擬和優(yōu)化設計各向異性光子晶體元件的過程中需要有效的數(shù)值方法。傳統(tǒng)的數(shù)值方法有各種局限性，比如:時域有限差分方法［6］需要很龐大的計算資源，且不容易算出高精度的結(jié)果;有限元方法［7］通常需要求解一個很大的方程組，不管是直接求解還是迭代求解都很繁瑣。

Xie等在文獻［8］中提出:可以利用各向異性光子晶體單元晶格的DtN映射來分析各向異性光子晶體的帶隙結(jié)構。DtN映射把單元晶格邊界上的波動場映射成自己的法向?qū)?shù)。在實際計算中，可以用一個矩陣來近似，這樣就避免了計算單元晶格內(nèi)部的波動場，大大減少了計算量。

本文將利用各向異性光子晶體單元晶格的DtN映射來計算周期排列的液晶柱的透射/反射光譜。這是一個邊界值問題，通過利用合適的邊界條件截取較小的計算區(qū)域，再利用DtN映射進一步減少未知數(shù)個數(shù)，最終得到一個較小的線性方程組并進行求解。這種方法可用于模擬更多的各向異性光子晶體元件。

1 控制方程

本文研究的是光波在周期排列液晶柱上的透射/反射問題。如圖1所示，整個結(jié)構是由液晶柱按照正方形晶格排列而成。其中，圓表示周期排列的液晶柱在xy平面上的橫截面，整個結(jié)構在z方向上沒有變化，光波從右側(cè)射入，從左側(cè)射出，共穿過m層液晶柱。

圖1 m層周期排列的液晶柱

線性電磁波(光波)在非導電、非磁性、無電流、無電荷的各向異性介質(zhì)中傳播，滿足麥克斯韋方程組，考慮時間諧波，對于各向異性的光子晶體，相對電容率可表示成矩陣假設z軸是各向異性介質(zhì)的一根主軸，相對電容率可寫成:

另外，由于圖1所示的整個結(jié)構在z方向上沒有發(fā)生變化，此時控制方程可表示為:

E偏振模式下的控制方程(2)和各項同性介質(zhì)的控制方程一致，在文獻［9］中已有研究。因此，本文只考慮H的偏振情形，控制方程為方程(3)。

2 線性方程組的建立

2.1 邊界條件

為了得到較小的計算區(qū)域，需要合適的邊界條件。圖1中實線段圍成區(qū)域表示結(jié)構中的一個周期，也是截取的計算區(qū)域。其中:左側(cè)邊界上y=0，在y＜0區(qū)域中材料的折射率為n1;右側(cè)邊界上y=mL，L為晶格常數(shù)，在y＞mL區(qū)域中材料的折射率為n2。

假設光波從y=mL的右側(cè)區(qū)域射入，入射波u(i)和 y軸的夾角是 θ0，可表示為u(i)(x，y)=ei［α0x－β0(y－mL)］，其中:α0=k0n2sin(θ0)，β0=k0n2·cos(θ0)。

反射波u(r)在y＞mL區(qū)域可表示為

傳輸波u(t)在y＜0區(qū)域可表示為

在如圖1所示的計算區(qū)域外部的4條邊界上(實線段)，邊界條件的建立過程與文獻［9］類似。

在左側(cè)邊界(y=0)上，v0，u為控制方程(3)的解，定義算子，使得±2，… ，從而得到邊界條件如下:

在右側(cè)邊界(y=mL)上，u=u(i)+u(r)。定義算子 s0，使得:s0eiαjx= βjeiαjx，j=0，± 1，± 2，…，從而得到邊界條件:

在上側(cè)邊界(x=L)和下側(cè)邊界(x=0)處，利用結(jié)構的周期性，可以得到邊界條件:

其中ρ=eiα0L。

2.2 各向異性光子晶體單元晶格的Dirichlet-to-Neumann映射

圖1計算區(qū)域中的每一個小正方形表示一個單元晶格，半徑是a的液晶柱位于正方形中間，正方形邊長為L。L也是2個相鄰液晶柱中心點之間的距離。

各向異性光子晶體單元晶格的DtN算子的構造過程在文獻［8］中已有提及。主要構造過程為:首先對電容率矩陣的子矩陣作對角化，然后旋轉(zhuǎn)、拉伸坐標系，方程(5)可轉(zhuǎn)化為標準的Helmholtz方程:

這樣，方程(3)的通解為

其中:

這里:n0是液晶柱以外的材料的折射率;是待定系數(shù)。利用Hz和電場的切向分量在邊界面r=a處連續(xù)，再利用適當?shù)腇ourier展開，可以得到。令u=Hz，Γ為單元晶格Ω

DtN映射Λ將單元晶格邊界上的波動場映射成其法向?qū)?shù)，一旦在單元晶格Ω的邊界Γ的每條邊上離散N個點，DtN映射Λ就可用4N×4N階矩陣近似。

2.3 線性方程組的建立

在圖1所示的計算區(qū)域中，利用各向異性光子晶體單元晶格的DtN算子Λ避免了對單元晶格內(nèi)部的計算，只需計算單元晶格邊上的波動場。單元晶格的邊分為2類:一類位于計算區(qū)域的內(nèi)部(虛線段);另一類位于計算區(qū)域的邊界上(實線段)。對計算區(qū)域內(nèi)部單元晶格的邊，可用DtN映射建立起關于波動場的方程;對邊界上單元晶格的邊，可用DtN映射和邊界條件建立起關于波動場的方程。

以v1所在的邊為例來說明如何對內(nèi)部邊(虛線段)建立方程。由于這條邊屬于單元晶格Ω1，Ω1的 DtN 映射可表示為:Λ［v0，u01，v1，u11］T=把Λ分為4×4分塊矩陣，得到v1法向?qū)?shù)的表達式:

其中u01，v0，u11，v1分別表示 Ω1的上、左、下、右4條邊上的波動場。另外，v1所在的邊還屬于單元晶格Ω2，從而又能利用Ω2的DtN映射得到v1法向?qū)?shù)的一個表達式，而且由于Ω1和Ω2為相同的單元晶格，它們的DtN映射是完全一樣的，有的邊界，從u|Γ和?vu|Γ中消去未知系數(shù)cm可以得到DtN映射Λ:

其中u02，v1，u12，v2分別表示 Ω2的上、左、下、右 4條邊上的波動場。由式(12)、(13)可得到方程:

類似地，在計算域內(nèi)部單元晶格的每條邊上(虛線段)都可以建立類似式(14)的方程。

以v0所在的邊為例來說明如何對邊界上單元晶格的邊(實線段)建立方程，在左側(cè)邊界(y=0)上，由于這條邊位于單元晶格Ω1內(nèi)，利用Ω1的DtN算子Λ，v0的法向?qū)?shù)可表示為

聯(lián)立式(6)表示的邊界條件得到方程:

類似地，對邊界上單元晶格的每條邊上都可建立類似式(16)的方程。最終，可以得到一個線性方程組:

若單元晶格每條邊離散N個點，則系數(shù)矩陣A的大小為(3m+1)N×(3m+1)N。

從Ax=b中解出v0和vm(v0是y=0處的波動場，vm是y=mL處的波動場)，利用式(4)、(5)分別求出R0和T0，從而得到透射/反射率。

3 數(shù)值算例

當z軸為主軸，各向異性液晶的電容率矩陣可寫成矩陣(1)，此時液晶是單軸介質(zhì)。矩陣(1)有3個特征值，其中2個相等特征值與平凡折射率no有關，另一個不等的特征值與非凡折射率ne有關，液晶的光軸和x軸的夾角用φ表示，則矩陣(1)中的元素可表示為:

為了驗證算法的有效性，分別考慮真空中呈周期排列的液晶柱的透射/反射光譜和液晶填充的多孔板的透射/反射光譜2種情況。

首先，模擬光波在真空中呈周期排列的液晶柱中的傳播狀況。結(jié)構如圖1所示，圓形表示周期排列的液晶柱，圓外是真空。真空中折射率n0=1;y＜0處的折射率n1=1;y＞mL處的折射率n2=1。光波從右側(cè)射入。參數(shù)的選擇和文獻［8］一致，層數(shù) m=16，半徑 a=0.3L，ne=2.223，no=1.590，單元晶格每條邊上離散的點數(shù)取N=15，光軸與x軸的夾角為，算出透射/反射光譜如

圖2(a)、(b)所示。

為了模擬外界因素的變化對光子晶體元件的影響，其他參數(shù)保持不變，通過調(diào)整光軸與x軸的夾角 φ，取，算出透射/反射光譜，如圖2(c)、(d)所示。觀察對比透射光譜圖2(a)和(c)，發(fā)現(xiàn)光的傳播情況在0.4附近有明顯變化，表現(xiàn)在某些頻率的光波由不能在各向異性光子晶體中傳播變成可以傳播。比如頻率為0.391 80的光波，當時，透射率為0.006 6;當時，透射率為0.990 1。而頻率為0.392 75的光波，當時，透射率為0.005 7;當時，透射率為0.991 5。

第2個數(shù)值算例是模擬光波在液晶填充的多孔板中的傳播狀況。結(jié)構如圖1所示，圓形表示用液晶填充的空氣洞，圓外是硅板。硅板的折射率n0=3.4;y＜0處的折射率n1=1;y＞mL處的折射率n2=1。光波從右側(cè)射入。層數(shù)m=16，半徑a=0.4L，L為晶格常數(shù)，ne=1.707 2，no=1.529 2，，在單元晶格的每條邊上離散19個點，得到透射/反射光譜如圖3(a)、(b)所示。

為了模擬外界因素的變化對光子晶體元件的影響，通過調(diào)整光軸與x軸的夾角φ，取，計算出透射/反射光譜，如圖3(c)、(d)所示，觀察對比透射光譜圖3(a)和(c)，發(fā)現(xiàn)光的傳播情況在0.24附近透射率有較大變化，表現(xiàn)在某些頻率的波由可以在各向異性光子晶體中傳播變成不能傳播。比如頻率為0.239 80的光波，當時的透射率為0.992 7時的透射率僅為0.000 2;而頻率為0.242 45的光波，當時的透射率為0.990 7，而當時的透射率僅為0.003 1。也就是說，施加外界因素影響可以調(diào)控特定頻率光波的傳播情況，使得光波在某些條件下可以傳播，在某些條件下不能傳播。

圖2 φ=和φ=時真空中呈周期排列液晶柱的透射/反射光譜

圖3 φ=和φ=時液晶填充的多孔板的透射/反射光譜

4 結(jié)束語

本文主要研究基于各向異性光子晶體單元晶格的Dirichlet-to-Neumann(DtN)映射，發(fā)展了一種能有效分析周期排列液晶柱的透射/反射光譜的數(shù)值方法。通過利用DtN算子避免了晶格內(nèi)部的計算，這樣只需計算晶格邊界上的波動場，使得未知數(shù)的數(shù)量大大減少。因此，這是一種快速高效的數(shù)值算法，可推廣用于模擬各種各向異性光子晶體元件，并進一步用于其優(yōu)化設計。

［1］JOANNOPOULOS J D，MEADE R D，WINN J N.Photonic Crystals:Modeling the Flow of Light［M］.Princeton，NJ::Princeton Uiversity Press，1995.

［2］龍濤，劉啟能.各向異性矩形光子晶體禁帶結(jié)構及量子效應［J］.激光與紅外，2011，41(2):197－201.

［3］慕振峰.各向異性介質(zhì)波導及含左手材料光子晶體的傳輸特性［D］.北京:北京工業(yè)大學，2013.

［4］FRANCESCA S，SHANE M E，ROBERTO C.Nematic Liquid Crystals Embedded in Cubic Microlattices:Memory Effects and Bistable Pixels［J］.Adv Funct Mater，2013，23(32):3990－3994.

［5］KATSUMI Y，YUKI S，YOSHIAKI K.Temperature tuning of the stop band in transmission spectra of liquid-crystal infiltrated synthetic opal as tunable photonic crystal［J］.Applied PhysicsLetters，1999，75(7):932－936.

［6］SINGH G，TAN E L，CHEN Z N.Efficient Complex Envelope ADI-FDTD Method for the Analysis of Anisotropic Photonic Crystals［J］.IEEE Photonics Technology Letters，2011，24(12):801－803.

［7］CHERBI L，BELLALIA A，UTHMAN M，et al.Modeling of two rings-photonic crystal fibers with the scalar-finite element method［J］.Journal of optoele ctronics and advanced materials，2013(15):1385－1391.

［8］XIE H，LU Y Y.Modeling two-dimensional anisotropic photonic crystals Dirichlet-to-Neuman-n maps［J］.Journal of the optical society of America a-optics image science and vision，2009，26(7):1607－1614.

［9］HUANG Y X，LU Y Y.Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps［J］.Journal of lightwave technology，2006，24(9):3448－3453.