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

蔣 丹，胡 真

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

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

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

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

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

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

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

1 控制方程

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

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

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

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

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

2 線性方程組的建立

2.1 邊界條件

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

假設(shè)光波從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所示的計(jì)算區(qū)域外部的4條邊界上(實(shí)線段)，邊界條件的建立過(guò)程與文獻(xiàn)［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é)構(gòu)的周期性，可以得到邊界條件:

其中ρ=eiα0L。

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

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

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

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

其中:

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

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

2.3 線性方程組的建立

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

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

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

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

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

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

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

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

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

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

3 數(shù)值算例

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

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

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

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

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

第2個(gè)數(shù)值算例是模擬光波在液晶填充的多孔板中的傳播狀況。結(jié)構(gòu)如圖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個(gè)點(diǎn)，得到透射/反射光譜如圖3(a)、(b)所示。

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

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

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

4 結(jié)束語(yǔ)

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

［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é)構(gòu)及量子效應(yīng)［J］.激光與紅外，2011，41(2):197－201.

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