基于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的系数偏微分模块的内置算法直接对方程进行数值计算,该模块完全基于数学方程求解,具有极大的灵活性,对于一些新提出的数学模型有较好的借鉴意义。这种方法为研究更复杂光子晶体模型打造光学类比平台,提供了一种更为便捷的数值仿真手段。