基于DtN映射計算周期排列液晶柱的透射/反射光譜

蔣 丹，胡 真

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

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

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

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

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

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

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

1 控制方程

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

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

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

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

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

2 線性方程組的建立

2.1 邊界條件

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

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

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

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

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

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

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

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

其中ρ=eiα0L。

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

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

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

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

其中:

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

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

2.3 線性方程組的建立

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

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

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

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

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

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

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

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

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

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

3 數值算例

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

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

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

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

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

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

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

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

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

4 結束語

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

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

［3］慕振峰.各向異性介質波導及含左手材料光子晶體的傳輸特性［D］.北京:北京工業大學，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.