基于COMSOL的光子晶体能带结构仿真计算

2019年第38卷第8期传感器与微系统（Transducer and Microsystem Technologies）

DOI：10．13873/J．1000—9787（2019）08—0111—03

基于COMSOL的光子晶体能带结构仿真计算*

付子义1，王晨旭1，长谷川弘治2

（1．河南理工大学电气工程与自动化学院，河南焦作454000；

2．室兰工业大学情报电子工学系，日本室兰0500071）

摘要：为更方便、快捷地求解光子晶体能带结构问题，基于有限元数值仿真软件COMSOL Multiphysics

系数偏微分方程模块，重新建立了数学模型，直接从亥姆霍兹方程出发，结合布洛赫态，推导出光子晶体偏

微分形式的本征方程，充分考虑了布洛赫波矢的传播情况，求解出相应的本征频率从而求解出能带，对

一维和二维光子晶体能带结构分别进行仿真计算。仿真计算结果与传统方法结果进行对比分析，验证了

此方案的可行性。

关键词：光子晶体；偏微分方程；本征方程；能带结构；COMSOL Multiphysics

中图分类号：TP391文献标识码：A文章编号：1000—9787（2019）08—0111—03

Simulation computation of photonic crystals energy

band structure based on COMSOL*

FU Ziyi1，WANG Chenxu1，HASEGAWA Koji2

（1．School of Electrical Engineering and Automation，Henan Polytechnic University，Jiaozuo454000，China；

2．Information and Electronic Engineering，Muroran Institute of Technology，Muroran0500071/Hokkaido，Japan）

Abstract：In order to solve the problem of photonic crystal band structure more conveniently and quickly，based

on the COMSOL Multiphysics coefficient partial differential equation（PDE）module of the finite element

numerical simulation software，the mathematical model is reestablished．The partial differential form of the photonic

crystal is derived from the Helmholtz equation and combined with the Bloch state．The eigenfunction，and the

propagation of Bloch wave vector is taken into full consideration，the corresponding eigenfrequency is solved so as

to solve the energy band．In this case，the energy band structure of one and two-dimensional photonic crystals is

simulated respectively．The simulation results are compared with the results of traditional methods，and the

feasibility of the scheme is verified．

Keywords：photonic crystals；partial differential equation（PDE）；eigenfunction；energy band structure；

COMSOL Multiphysics

0引言

当电磁波在光子晶体［1，2］中传播时由于布拉格衍射的影响，会受到调制而形成能带结构，即光子能带（photonic band），光子能带之间可能出现的带隙，即光子带隙（pho-tonic band gap，PBG），频率处于光子能带里的电磁波可以在光子晶体中几乎无损地传播，但是出于光子带隙的电磁波，却不能在光子晶体中传播。这一特性使得光子晶体成为一种理想的、可以按照人们的意愿来操控电磁波的工具。

基于光子晶体能带结构的计算，文献［3］将介电常数和电磁场以平面波的形式展开，将麦克斯韦方程组化成一个本征方程。文献［4］以差分原理为基础，从麦克斯韦旋度方程出发，将其转换为差分方程组，在一定体积内和一段时间上对连续电磁场的数据取样。文献［5］从麦克斯韦方程出发，经过变分原理得到微分方程的等效积分弱形式，再经过分片差值，分片求积得到代数方程组，最后对代数方程组进行求解计算。以上方法都是从麦克斯韦方程出发进行推导求解，考虑到光子晶体中电磁波处于时谐场，还需要进行时谐处理，结合布洛赫态得到一个矢量方程，还需要将梯度算子进行变换，转换成标量形式的本征方程［6］，推导过程繁琐，对光子晶体物理概念的理解帮助不大。并且没有对一维光子晶体中布洛赫波矢不沿坐标轴传播的情况和二维光子晶体中布洛赫波矢落于布里渊区域内的情况进行分析说明。

收稿日期：2018—06—22

*基金项目：国家自然科学基金资助项目（61501175）

111

传感器与微系统

第38卷

不同于以往从麦克斯韦方程出发，

本文直接由空间部分的亥姆霍兹方程出发，

分析电磁波的传播特性，将矢量形式的布洛赫态转化为标量形式的布洛赫态，推导出偏微分形式的本征方程简化了数学建模过程

［7］

，并推导分析了布

洛赫波矢传播的所有情况，

然后基于有限元仿真软件COMSOL 的系数偏微分模块求解本征值。1

光子晶体能带结构的数学推导

在光子晶体中，研究的对象是电磁波，电磁波的行为可

以由亥姆霍兹（Helmholtz ）方程［8］

准确描述

［Δ

2

+ω2ε（r ）μ］E （r ）=0（1）

式中

ω为电磁波的的频率，由于光子晶体的介电常数在

空间周期性分布，

ε（r ）为周期函数，称为相对介电函数，μ为磁导率，

E （r ）为电场强度。1．1一维光子晶体

本文的研究对象一维光子晶体是一种多层膜结构，这

个系统在xy 平面上是匀质的，

在Z 方向上具有周期性，并且由介电常数ε1，

ε2交替组成，空间周期为a 。在本文中，将一维光子晶体分为两种情况来进行仿真计算，一种是沿坐标轴传播，

一种是不在坐标轴上传播，分别如图1所示

。图1二种一维光子晶体仿真方式

光子晶体晶格的周期性通过在光子晶体晶胞边界充分使用布洛赫态来保证。

由周期函数E （z +a ）=E （z ）

（2）

可以得到布洛赫态E （r ）=^y

E （z ）e －ik ·r 。在这里k ·r =k z ^

z ·r =k z z （3）可以得出标量形式的布洛赫态

E （r ）=^y

E （z ）e －ik z z （4）将得到的布洛赫态代入到亥姆霍兹方程中［d 2d z

2－2j k z d d z －k 2z +ω2

ε（z ）μ］E （z ）=0（5）作为本文理论模型的验证，采用文献［6］所用参数，以便对比。在文献［

6］中，频率由无量纲量（ωa /2πc ）所表示，这里将上式进行归一化处理

a 2［d 2d z

2－2j k z d d z －k 2z +ω2ε（z ）μ］E （z ）=0（6）

接下来讨论不在坐标轴上传播的情况如图1（b ）所示，同理，得出布洛赫态

E （r ）=^x

E （z ）e －ik y y e －ik z z （7）

将布洛赫态代入亥姆霍兹方程中

［d 2dz

2－2j k z d d z －k 2z －k 2y +ω2

ε（z ）μ］E （z ）=0（8）

将上式进行归一化处理得到

［d 2

d （z a ）（z a ）－2j k z a d

d （z

a

）

－（k z a ）2－（k y a ）2+

（

ωa c

）2

εr μr ］E （z ）=0（9）

1．2二维光子晶体

本文研究对象为简单立方晶格结构的光子晶体，每

一个光子晶体原胞只含有一个散射体［9］

。这里将无限长

电介质圆柱体放置在空气电介质背景材料中周期性排列组成二维光子晶体。圆柱体横截面半径为r ，

晶格常数为a ，r =0．2a 。圆柱体轴线与z 轴平行，二维周期性结构分布在x －y 坐标平面内。

二维光子晶体的布洛赫态E （r ）=E （x ，y ）e －ik x x e －ik y y

（10）

将布洛赫态代入亥姆霍兹方程中可得

［ 2

x 2－2j k x x －k 2x + 2

y

2－2j k y

y －k 2y +ω2ε（x ，y ）］E （x ，y ）=0（11）

对上式进行归一化处理

［ 2

（x a ）（x a ）－2j k x a （x a ）－（k x a ）2+

2

（y a ）（y a

）－

2j k y a

（y a

）－（k y a ）2+（ωa c ）2

εr μr ］E （x ，y ）=0（12）

二维正方晶格的基矢为a →1=a ^x ，a →

2=a ^y ，则倒易空间的

基矢为b →1=（2π/a ）^x

，b →2=（2π/a ）^y 。充分利用光子晶体的对称性，消除冗余，将波矢限制在第一布里渊区，如图2（a ）所示。

考虑Q 点在第一布里渊区域内

k x =|k |cos θ^x =π2a ^x ，k y =|k |sin θ^y

=槡3π6a ^y （13）

考虑波矢沿着第一布里渊区的边界移动由于带隙往往产生于第一布里渊区域［10］

的边界处，因

此，只需要计算仿真出边界处的值就可以表现出全部需要

研究的能带情况

k x =

12b 1=πa ^x ，k y =12

b 1=

πa ^

y

（14）

图2

二维晶体二种仿真情况

211