基于Dirichlet-to-Neumann映射计算旋磁光子晶体的带隙结构*

杨华娇， 胡真

(河海大学 理学院，江苏 南京 211100)

光子晶体(Photonic Crystals)是一种新型的周期性人造材料，周期为光波长量级[1].光子晶体最重要的性质之一是具有光子带隙(Photonic Band Gap)，频率落在光子带隙范围内的光波将无法在光子晶体内部进行传播，从而光子晶体具有了控制光的传播的能力，被制成各种光子晶体元件应用于各个领域，如无损耗光波导[2]、高辐射亮度超微激光器[3]、高光电转化效率薄膜光伏电池[4]等.

在过去的二十年中，对于二维光子晶体的研究主要集中在各向同性材料制成的光子晶体.近年来，人们发现基于旋磁各向异性材料制作的二维旋磁光子晶体具有更多特殊的性质.旋磁各向异性材料是指该材料的相对磁导率会随着方向的变化而改变.利用二维旋磁光子晶体，Wang等设计并制作出了只能单向传播的波导管[5-6]；He等设计了一种单向交叉波导分路器[7]；Chen等设计了一种无色散的慢光波导管[8].

分析设计二维旋磁光子晶体最基本的问题是计算旋磁光子晶体的带隙结构.计算带隙结构需要有效的数值方法，传统的数值方法有：平面波展开法[9]、有限元法[10]、多重散射法[11]等等.但是这些方法各有不足：平面波展开法和有限元方法对于非色散介质而言，建立的是线性特征值问题，但对于色散介质建立的却是非线性特征值问题；多重散射法需要格点求和技术，计算过程复杂.

近年来Lu等[15]提出了一种能有效分析各种光子晶体结构的新的数值方法，叫作Dirichlet-to-Neumann(DtN)映射方法.单元晶格的DtN映射将单元晶格边界上的波动场映射成边界上的法向导数，通过DtN映射避免了在晶格内部的离散，从而达到简化运算的效果.DtN映射方法已经被广泛用于分析二维各向同性光子晶体结构，

本文将DtN映射方法应用于计算二维旋磁光子晶体的带隙结构.旋磁光子晶体的相对磁导率μ为一个矩阵，控制方程及不同介质边界面上的连续性条件均比较复杂，因此需要重新在旋磁光子晶体的单元晶格上构造出相应的DtN映射，再基于该DtN映射建立特征值问题来求解旋磁光子晶体的带隙结构.利用扩展的DtN映射建立的特征值问题以波数β的函数为特征值，以频率ω为给定的参数.通过旋磁光子晶体单元晶格的DtN映射，得到较小矩阵的特征值问题，且无论是色散介质还是非色散介质，该方法建立的特征值问题均为线性的.

1 基本方程

旋磁光子晶体的相对磁导率μ为3×3维的对称矩阵

(1)

本文仅考虑二维问题，且只研究TE模式，利用Maxwell方程组，可得TE模式下的控制方程为

(2)

其中u为电场在z方向上的分量，是二维梯度算子，

2 二维旋磁光子晶体单元晶格的DtN映射

为了建立特征值问题来求解二维旋磁光子晶体的带隙结构，需要首先构造出其单元晶格的DtN映射.单元晶格的DtN映射Λ将单元晶格边界Γ上的波动场u映射成u在Γ上的法向导数，即

Λu|Γ=∂vu|Γ

(3)

其中v表示边界Γ上的单位法向量.旋磁光子晶体是由旋磁各向异性材料在空间中周期排列而构成的人造材料，常用的周期排列方式有正方形、三角形和蜂巢状三种(如图1所示)；图中线段框出的部分为对应结构的单元晶格，a为晶格常数，其中正方形和三角形周期排列的磁光子晶体的单元晶格都是只包含一个圆柱的晶格，他们的DtN映射可以直接构造；蜂巢状周期排列的旋磁光子晶体的单元晶格由三个小正六边形晶格组成(两个包含圆柱，一个不包含圆柱)，其DtN映射可以利用这三个小晶格的DtN映射构造得到.

(a)正方形周期排列介质柱型 (b)三角形周期排列薄板型 (c)蜂巢状周期排列介质柱型

图1 旋磁光子晶体

Fig.1 Gyromagnetic photonic crystals

2.1 只包含一个圆柱的晶格的DtN映射

控制方程(4)在只含有一个圆柱的结构上的通解u(x,y)可以用特解φj(r)eijθ的线性组合来表示，即

(4)

其中，cj为未知系数，(r,θ)是点(x,y)在极坐标下的坐标，φj(r)与Bessel函数Jj、Yj有关，表示为

(5)

Jj(q)Aj+Yj(q)Bj=Jj(s)

(6)

其中，s=k0n1a,q=k0n2a.

(7)

为了构造晶格的DtN映射，还需要将通解u(x,y)的法向导数用特解法向导数的线性组合表示出来：

(8)

由(4)式和(8)式消去未知系数cj,从而得到只包含一个圆柱的晶格的DtN映射Λ.例如，假设在边界Γ上共离散p个点(p为偶数)，离散点的坐标为(xk,yk),(rk,θk)为该点的极坐标，其中k=1,2,…,p,则由(4)式得到

(9)

再由(8)式得到

(10)

将(9)和(10)改写为矩阵形式

(11)

Λ=BAV-1

若在正边形晶格的每条边上离散N个点，则p=4N,Λ近似为(4N)×(4N)的矩阵；若在六边形晶格的每条边上离散N个点，则p=6N,Λ近似为(6N)×(6N)的矩阵.

2.2 不包含圆柱的晶格的DtN映射

对不包含圆柱的晶格而言，其DtN映射的构造过程和包含圆柱的晶格的DtN映射构造过程类似，唯一的区别在于方程(4)的特解更加简单，为

φj(r)=Jj(k0nr)

2.3 蜂巢状单元晶格的DtN映射

蜂巢状单元晶格由两个只包含一个圆柱的小晶格和一个不包含圆柱的小晶格组成.为了构造出蜂巢状单元晶格的DtN映射，需要先建立Ω2、Ω3(如图2(c)所示)这两个小晶格的DtN映射Λ(2)、Λ(3),再利用2.2节中介绍的方法构造不包含圆柱的晶格Ω1的映射Λ(1).利用这三个小晶格的DtN映射可以在公共边上建立如下方程：

A1V=A2U

(12)

其中，U为蜂巢状单元晶边界上的波动场，V为三条公共边上的波动场.

蜂巢状单元晶格边界上波动场的法向导数可以由对应小晶格的DtN映射得到

∂vU=B1U+B2V

(13)

通过方程(12)和(13)消去公共边上的波动场V,从而得到蜂巢状单元晶格的DtN映射

3 特征值问题

构造了二维旋磁光子晶体单元晶格的DtN映射之后，便可以建立求解旋磁光子晶体带隙结构的特征值问题.这一建立过程与求解各向同性光子晶体带隙结构的建立过程类似，其中求解正方形、三角形和蜂巢状周期排列的各向同性光子晶体带隙结构的特征值问题分别在Yuan[16]、Yuan[17]和李静[18]的文章中有详细介绍.下文简要介绍这一过程.

分析旋磁光子晶体的带隙结构，需要求解方程(4)在相应结构中的Bloch模式解：

u(x,y)=ei(αx+βy)Ψ(x,y)

其中，(α,β)为Bloch波矢，Ψ为周期函数.

将在第2节中得到的单元晶格的DtN映射Λ按照单元晶格的边数m进行分块，使其每一个分块均为N×N的矩阵，且记单元晶格每条边上的波动场分别为u1,u2,…,um(正方形m=4；六边形m=6,蜂巢状m=12)，从而可将(3)式改写为

(14)

再分别利用正方形、三角形和蜂巢状周期排列的旋磁光子晶体的拟周期边界条件和(14)式，便可以得到求解他们带隙结构的特征值问题.由于这三种结构建立特征值问题的过程类似，只以正方形周期排列的旋磁光子晶体结构为例说明此过程.

图2 正方形周期排列旋磁光子晶体的布里渊区Fig.2 Brillouin zone of a square photonic crystals

正方形单元晶格的拟周期边界条件为：

u3=ρβu1,∂vu3=ρβ∂vu1；u4=ραu2,∂vu4=ρα∂vu2

(15)

其中，ρα=eiαL,ρβ=eiβL.

将(13)式代入(12)式，便可得到求解正方形周期排列旋磁光子晶体带隙结构的特征值问题：

其中，I为单位矩阵，V=[λUT,UT]T,U=[u1,u2]T,特征值λ和矩阵M0、M1、M2的取值与正方形周期排列旋磁光子晶体的布里渊区有关.

在ΓX、XM、MΓ边上特征值λ与矩阵M0、M1、M2的取值分别如下所示：

(1)从Γ到X,λ=eiαa(0<α≤π/a)且

(2)从X到M,λ=eiβa(0<β≤π/a)且

(3)从M到Γ,λ=eiαa=eiβa(0<α=β≤π/a)且

建立求解三角形和蜂巢状周期排列的旋磁光子晶体带隙结构的特征值问题的过程与此类似，差别仅在于拟周期边界条件以及相对应的布里渊区略有不同.

4 数值算例

分别以正方形、三角形和蜂巢状周期排列的二维旋磁光子晶体为例验证方法的有效性，并且基于该方法优化设计出了带隙更宽的蜂巢状周期排列的旋磁光子晶体.

4.1 正方形周期排列的介质柱型旋磁光子晶体

首先计算二维正方形周期排列的介质柱型旋磁光子晶体的带隙结构.Wang等用有限元法也计算过这一旋磁光子晶体的带隙结构，旋磁各向异性材料采用的是YIG材料[5].YIG的相对介电常数为εr=15,介质柱半径为r=0.11a,a为晶格常数.当外加磁场的作用下引发YIG强烈的各向异性，此时μr=14,μk=-12.4.计算结果如图3所示，阴影部分表示光子带隙，标准化频率的带隙范围为[0.324 2,0.448]、[0.527,0.576]和[0.611,0.647 5],结果与Wang等一致.

图3 计算结果 图4 计算结果

4.2 三角形周期排列的薄板型旋磁光子晶体

接下来，计算二维三角形周期排列的薄板型旋磁光子晶体的带隙结构.Lu等用有限元法和平面波展开法计算过这一旋磁光子晶体的带隙结构，采用的旋磁材料也为YIG材料，材料参数与4.1中相同，且空气洞半径为r=0.44a[14].同样，在外加磁场的作用下光子晶体磁化.计算结果如图4所示，阴影部分表示光子带隙，标准化频率的带隙范围为[0.341 5,0.361 1],与Lu等一致.

4.3 蜂巢状周期排列的介质柱型旋磁光子晶体

最后，计算二维蜂巢状周期排列的介质柱型旋磁光子晶体的带隙结构.Ao等用多重散射法计算过这一旋磁光子晶体的带隙结构，采用的旋磁各向异性材料为铁氧体，相对介电常数εm=15,介质柱半径r=0.2a(其中a=10 mm)[13].铁氧体光子晶体被外加的恒定磁场沿+z方向饱和磁化后便会呈现出各向异性的性质，其中H0=500 Oe,H0是外加磁场，4πMs=1 750 G(4πMs是饱和磁化强度

ω0=γH0是进动频率(旋磁比γ=2.8×106rad·s-1·G-1).计算结果如图5所示，第三条能带和第四条能带间的标准化的带隙范围为[0.252 9,0.272 5]，与Ao等一致.

由于蜂巢状旋磁光子晶体的单元晶格由三个小正六边形晶格构成，所以有更多的结构参数可供调整，特别适合进行优化设计.下面尝试对图1(c)所示的蜂巢状旋磁光子晶体结构进行优化，仅通过调整小晶格Ω2、Ω3中圆柱的半径大小，来扩大其带隙范围.当r2=r3=0.27 a时，得到了较宽的带隙结构，计算结果如图6所示，第二条能带和第三条能带间的带隙范围扩大到了[0.221 3,0.257 5].

图5 计算结果 图6计算结果

5 结 语

扩展了DtN映射方法，建立了二维旋磁光子晶体单元晶格的DtN映射，并将其用于计算正方形、三角形及蜂巢状周期排列的旋磁光子晶体的带隙结构，方法避免了在单元晶格内部的离散，只需要在边界上进行离散，从而得到矩阵较小的特征值问题，且计算快速并易于调整参数，非常有利于优化设计出具有更宽带隙的二维旋磁光子晶体.目前考虑的旋磁光子晶体都是由圆柱周期排列构成的，柱体的截面均为圆形，所以可以使用柱面波展开法来构造DtN映射，接下来将进一步研究非圆柱构成的旋磁光子晶体的带隙结构，如柱体的截面为三角形、椭圆等其他形状，这时需要重新设法来构造单元晶格的DtN映射.此外，还可以研究三维旋磁光子晶体的带隙结构，比如在有限厚度的板上周期性的挖取空气洞形成的结构等.