基于COMSOL的光子晶體能帶結構仿真計算*

付子義, 王晨旭, 長谷川弘治

(1.河南理工大學 電氣工程與自動化學院,河南 焦作 454000； 2.室蘭工業大學 情報電子工學系,日本 室蘭 0500071)

0 引 言

當電磁波在光子晶體[1,2]中傳播時由于布拉格衍射的影響,會受到調制而形成能帶結構,即光子能帶(photonic band),光子能帶之間可能出現的帶隙,即光子帶隙(photonic band gap,PBG),頻率處于光子能帶里的電磁波可以在光子晶體中幾乎無損地傳播,但是出于光子帶隙的電磁波,卻不能在光子晶體中傳播。這一特性使得光子晶體成為一種理想的、可以按照人們的意愿來操控電磁波的工具。

基于光子晶體能帶結構的計算,文獻[3]將介電常數和電磁場以平面波的形式展開,將麥克斯韋方程組化成一個本征方程。文獻[4]以差分原理為基礎,從麥克斯韋旋度方程出發,將其轉換為差分方程組,在一定體積內和一段時間上對連續電磁場的數據取樣。文獻[5]從麥克斯韋方程出發,經過變分原理得到微分方程的等效積分弱形式,再經過分片差值,分片求積得到代數方程組,最后對代數方程組進行求解計算。以上方法都是從麥克斯韋方程出發進行推導求解,考慮到光子晶體中電磁波處于時諧場,還需要進行時諧處理,結合布洛赫態得到一個矢量方程,還需要將梯度算子進行變換,轉換成標量形式的本征方程[6],推導過程繁瑣,對光子晶體物理概念的理解幫助不大。并且沒有對一維光子晶體中布洛赫波矢不沿坐標軸傳播的情況和二維光子晶體中布洛赫波矢落于布里淵區域內的情況進行分析說明。

不同于以往從麥克斯韋方程出發,本文直接由空間部分的亥姆霍茲方程出發,分析電磁波的傳播特性,將矢量形式的布洛赫態轉化為標量形式的布洛赫態,推導出偏微分形式的本征方程簡化了數學建模過程[7],并推導分析了布洛赫波矢傳播的所有情況,然后基于有限元仿真軟件COMSOL的系數偏微分模塊求解本征值。

1 光子晶體能帶結構的數學推導

在光子晶體中,研究的對象是電磁波,電磁波的行為可以由亥姆霍茲(Helmholtz)方程[8]準確描述

式中ω為電磁波的的頻率,由于光子晶體的介電常數在空間周期性分布,ε(r)為周期函數,稱為相對介電函數,μ為磁導率,E(r)為電場強度。

1.1 一維光子晶體

本文的研究對象一維光子晶體是一種多層膜結構,這個系統在xy平面上是勻質的,在Z方向上具有周期性,并且由介電常數ε1,ε2交替組成,空間周期為a。在本文中,將一維光子晶體分為兩種情況來進行仿真計算,一種是沿坐標軸傳播,一種是不在坐標軸上傳播,分別如圖1所示。

圖1 二種一維光子晶體仿真方式

光子晶體晶格的周期性通過在光子晶體晶胞邊界充分使用布洛赫態來保證。

由周期函數

E(z+a)=E(z)

(2)

可以得出標量形式的布洛赫態

將得到的布洛赫態代入到亥姆霍茲方程中

作為本文理論模型的驗證,采用文獻[6]所用參數,以便對比。在文獻[6]中,頻率由無量綱量(ωa/2πc)所表示,這里將上式進行歸一化處理

接下來討論不在坐標軸上傳播的情況如圖1(b)所示,同理,得出布洛赫態

將布洛赫態代入亥姆霍茲方程中

將上式進行歸一化處理得到

本文研究對象為簡單立方晶格結構的光子晶體,每一個光子晶體原胞只含有一個散射體[9]。這里將無限長電介質圓柱體放置在空氣電介質背景材料中周期性排列組成二維光子晶體。圓柱體橫截面半徑為r,晶格常數為a,r=0.2a。圓柱體軸線與z軸平行,二維周期性結構分布在x-y坐標平面內。

二維光子晶體的布洛赫態

E(r)=E(x,y)e-ikxxe-ikyy

(10)

將布洛赫態代入亥姆霍茲方程中可得

E(x,y)=0

(11)

對上式進行歸一化處理

考慮Q點在第一布里淵區域內

考慮波矢沿著第一布里淵區的邊界移動

由于帶隙往往產生于第一布里淵區域[10]的邊界處,因此,只需要計算仿真出邊界處的值就可以表現出全部需要研究的能帶情況

圖2 二維晶體二種仿真情況

2 計算仿真與分析

利用上一部分建立的數學模型,結合COMSOL系數偏微分(partial differential equation, PDE)模塊,可以開展對一維和二維光子晶體能帶結構的仿真計算。使用者可以根據實際問題在數學方程上進行修改,具有極大的靈活性,因此某些特殊的不便于用現成專題模塊描述的問題可通過數學模塊獲得合理的解決。

2.1 一維光子晶體能帶結構仿真

在一維光子晶體中,考慮到波矢完全沿著Z軸方向傳播。在這種情況下,k‖=0,因此只重點考慮波矢分量Kz即可。根據Comsol計算出的本征值,繪制出在坐標軸上傳播的兩種情況如圖3(a),(b)所示。當波矢量不沿著介電材料呈周期性變化的方向傳播時,在此考慮波矢沿y軸傳播,根據表3繪制出不在軸上傳播的情況如圖3(c)所示。在軸上傳播和不在軸上傳播的重點不同在于不存在帶隙。這是因為不在軸上傳播破壞了周期性結構。另一個不同涉及到能帶的簡并。由于兩種模式僅由晶體所具有的旋轉對稱性不同,所以它們一定會簡并。然而對于以任意波矢K方向傳播的模式,這個對稱性被破壞,因為這兩個帶之間不存在旋轉對稱性關系了,所以它們通常具有不同的頻率。

圖3 一維的光子晶體能帶結構

圖3中，ε0為真空介電常數。縱坐標頻率ω進行了歸一化處理由無量綱量(ωa/2πc)表示,其中a為空間周期,c為真空中電磁波的傳播速度。橫坐標波矢k同樣進行歸一化處理,由無量綱量(ka/2π)表示。

在一維光子晶體帶隙求解中,基于COMSOL系數偏微分模塊進行仿真計算時,將單元數分為100和1 000兩種情況,由以上結果可以看出,根據本文推導出的數學模型所計算的結果與文獻所用傳統方法所得結果基本一致,單元數的增加對仿真結果的影響可以忽略不計。

2.2 二維光子晶體能帶結構仿真計算

僅沿著第一布里淵區域邊緣繪制K‖的原因是給定帶狀(決定帶隙)的最大值和最小值總是出現在布里淵區域的邊緣。

圖4 二維光子晶體能帶結構

圖4中，ε0表示真空介電常數。縱坐標頻率ω進行了歸一化處理由無量綱量(ωa/2πc)表示,其中a為空間周期,c為真空中電磁波的傳播速度。橫坐標波矢k同樣進行歸一化處理,由無量綱量(ka/2π)表示。

在二維光子晶體結構中,布洛赫波矢被限制在第一布里淵區,并沿著第一布里淵區邊界運動(Γ-X-M),在此,將第一布里淵區分為三部分來進行仿真計算(Γ-X),(X-M),(Γ-M),由以上結果可以看出,根據本文推導出的數學模型所計算的結果與文獻所用傳統方法所得結果基本一致。

3 結 論

結合COMSOL Multiphysics系數偏微分模塊,重新建立了一套數學模型,通過與傳統方法所得結果的對比,驗證了此數學模型的正確性和可信性。并且對布洛赫波矢在光子晶體中的不同情況進行了推導說明。不同于以往從麥克斯韋方程出發,需要經過繁瑣的數學推導與變換,本文直接由空間部分的亥姆霍茲方程出發,推導出偏微分形式的本征方程,通過COMSOL的系數偏微分模塊的內置算法直接對方程進行數值計算,該模塊完全基于數學方程求解,具有極大的靈活性,對于一些新提出的數學模型有較好的借鑒意義。這種方法為研究更復雜光子晶體模型打造光學類比平臺,提供了一種更為便捷的數值仿真手段。