一种用于三维双各向异性光子晶体仿真的高效数值方法

《Computer Physics Communications》：A Novel Efficient Numerical Method for 3D Bi-Anisotropic Photonic Crystal Simulation

　　摘要 光子晶体因其能够精确控制光线的传播而广泛应用于光学通信和量子信息处理等领域。本文提出了一种基于斜Yee网格的、保持厄米正定性的插值策略，用于三维双各向异性光子晶体的计算。该方法将交错网格中各种电磁场分量的空间位置统一起来，从而显著减小离散化矩阵的规模。借鉴我们之前的工作，我们采用显式奇异值分解来压缩离散化特征值问题的零空间，进而为三维双各向异性光子晶体提供了一种高效的特征值求解方法。通过多个数值实验展示了所提方法的效率以及在高病态参数条件下的局域现象。本文相关的代码已发布在GitHub上（https://github.com/FAME-GPU/3D-Bi-Anisotropic-PhC）。

引言 光子晶体是一种具有周期性介电结构的人工微结构材料[1]、[2]。由于它们能够形成内部光子带隙，从而实现对电磁波传播的精确操控，因此受到了广泛关注[2]、[3]、[4]。正因为如此，光子晶体已被广泛应用于各种前沿领域，如高效光波导、新型激光器以及先进的传感装置[5]、[6]。除了传统的各向异性介质外，将双各向异性介质引入光子晶体为更复杂的灯光操控打开了新的大门[7]。近年来，双各向异性、手性以及磁光光子结构也受到了极大关注。这些结构的本构关系中的磁电耦合项能够同时调节电响应和磁响应，产生在各向同性材料中难以实现的新颖光学效应[8]。受复杂的磁电耦合效应影响，三维双各向异性光子晶体不仅表现出强烈的光局域现象[9]、[10]以及稳定的拓扑波导特性[4]、[11]，还与光学旋转、圆二色性、手性传感、非对称传播、拓扑边缘态以及超表面和超材料中的相位和偏振调控密切相关[12]、[13]。因此，准确计算三维双各向异性光子晶体的能带结构和本征模式，是分析和设计这类先进光子器件的基础前提。

在三维双各向异性光子晶体中，光的传播可以用频域下的无源麦克斯韦方程组来描述：?×E(r)=??ωB(r)，?·B(r)=0，?×H(r)=?ωD(r)，?·D(r)=0，其中E、H、D和B分别代表电场、磁场、电位移矢量和磁感应强度。虚数单位记为?=?1，ω为角频率。由材料的固有属性决定的本构关系描述了光子晶体内的磁电耦合效应，其表达式为D(r)=ε(r)E(r)+ξ(r)H(r)，B(r)=ζ(r)E(r)+μ(r)H(r)。材料参数的分布函数定义为ε(r)={εin,r∈Ω1,εout,r∈Ω0}，μ(r)={μin,r∈Ω1,μout,r∈Ω0}，ξ(r)={?γξin,r∈Ω1,0,r∈Ω0}，ζ(r)={??γζin,r∈Ω1,0,r∈Ω0}，其中Ω1和Ω0分别表示单元格内散射体与背景介质所占的子区域（如图1a所示）。矩阵εin、μin属于C3×3阶，代表散射体的相对介电常数和相对磁导率，而εout和μout则对应背景介质的相应参数，通常设为单位矩阵I3。变量ξin和ζin表示散射体的电磁耦合参数，需满足ξ=ζ*∈C3×3的条件，γ∈R则为手性参数[14]、[15]。此外，根据布洛赫定理[16]，电场E、磁场H、电位移矢量D和磁感应强度B都需满足准周期边界条件：F(x+a?)=e?2πk·a?F(x)，其中?=1,2,3，F∈{E,?H,?D,?B}，{a?}?=13∈R3表示光子晶体的晶格平移向量，2πk∈R3则表示第一布里渊区内的布洛赫波矢[3]。

构建并求解频域特征值问题(1)是计算光子学中获取特征频率ω(k)及相应本征模式的核心内容。现有的离散化方法各有优缺点。平面波展开法[17]、[18]概念上较为简单，但对于不连续介质收敛速度较慢，且在存在磁电耦合的双各向异性介质中计算成本较高。有限元法[19]、[20]虽然能够灵活处理复杂几何结构，但需要使用矢量边元素来消除虚假模式，这会导致自由度大幅增加，而双各向异性还会进一步加剧这一问题。多散射理论[21]、[22]在特定几何结构下效率很高，但难以适用于任意形状的散射体。

在此背景下，利用Yee网格进行空间离散化的有限差分法[23]得到了极为广泛的应用，因为它能够自然保留电磁场的旋度关系，实现方式简单，且与FFT加速技术兼容性良好。借助有限差分法的优势，我们此前曾用它来离散化麦克斯韦方程组，并对离散化系统进行了严格的矩阵分析[24]、[25]，包括采用显式奇异值分解和零空间自由变换等方法[26]，从而为三维光子晶体实现了高效的特征值求解器[27]。此外，我们还提出了一个内在的斜坐标系[28]，以便与准周期边界条件(3)完美匹配。

在充分利用有限差分法的优势以及我们之前的研究成果基础上，我们继续致力于开发针对三维双各向异性光子晶体的高效特征值求解器。Yee方案虽然严格满足了麦克斯韦方程组的无散度条件，但它会导致场采样点在空间上出现错位，因此无法直接用于离散化各向异性介质的本构关系(2)。[29]提出了一种双Lebedev方案，通过将标准Yee网格与三个经过位移处理的Yee网格及其对偶网格叠加在一起，来解决空间错位问题。尽管这种方法有效，但却使自由度增加了八倍，从而导致离散矩阵的规模大幅增大。[30]

为了解决离散化双各向异性本构关系时存在的固有难题，以及双Lebedev方案导致的矩阵维度急剧膨胀问题，本文提出了一种全新的高效方法。本研究的主要贡献如下：• 该离散化方案是在斜坐标系下构建的，能够保持电磁参数的正定性，从而保留连续系统的固有物理特性。• 相比于[29]中的双Lebedev方案，我们采用了插值策略来解决双各向异性本构关系中电磁场与介质参数的空间错位问题，使得离散化系统的维度仅为前者的1/8。• 基于我们之前的研究，我们提出了一种统一且高效的三维双各向异性光子晶体特征值求解框架，可高效处理正定和不定电磁参数系统。

在本文中，AT、A′和A*分别表示矩阵A的转置、共轭和共轭转置；若矩阵A是非奇异的，则A?T:=(A?1)T=(AT)?1；In为n维单位矩阵；O和0分别表示适当维数的零矩阵和零列向量。A⊕B表示矩阵A和B的直和；A?B表示矩阵A和B的克朗克积；vec(T)表示张量T的向量化；diag(a)用于构建以向量a的元素为对角线元素的対角矩阵；δji为克朗克德尔塔函数，当i=j时δji=1，否则为0。

章节片段 斜坐标系下麦克斯韦方程组的离散模型 由于光子晶体具有内在的周期性，根据布洛赫定理[16]，在晶格中传播的电磁波会满足准周期边界条件(3)。这一性质使得全空间的电磁场问题可以简化为单个单元格内的边值问题。典型的单元格几何结构及其计算域如图1a所示。为了保留光子晶体的内在周期性，采用离散化麦克斯韦方程组更为合理。

零空间自由的广义特征值问题 在离散化了旋度算子和本构关系之后，我们将(14)代入(7a)，从而得到麦克斯韦特征值问题(1)的离散化广义特征值问题(GEP)：[O??C*?CO][eh]=ω[NintXintZintMint]?1[eh]。对于任意的波矢k∈R3，GEP(16)本质上至少拥有2n维的零空间，因为算子C的零空间维度至少为n。GEP(16)中庞大的零空间在物理上代表了静态电磁模式。在光子带隙中……

数值实验 除非另有说明，以下所有实验均采用晶格常数a=1，网格尺寸n1=n2=n3=60，共计864,000个自由度。在所有实验中，各模型的背景介质均为空气，其相对介电常数和磁导率分别为ε=μ=I3。本文的所有实验均在配备有Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz处理器和256 GB内存的服务器上进行。相关算法是在MATLAB R2024b环境中实现的。

结论 总的来说，我们提出了一种用于三维双各向异性光子晶体特征值问题的高效算法。该算法采用一种新的、能够保持厄米正定性的插值策略，将交错分布的电磁场分量统一到相同的位置，从而使离散化矩阵的维度降至双Lebedev方案下的1/8。通过我们之前提出的显式奇异值分解方法，离散化旋度算子所对应的巨大零空间也被完全压缩。

未引用参考文献 [34]

CRediT作者贡献说明 Ke Gui：撰写——初稿、验证、软件实现、方法设计。Tiexiang Li：撰写——审稿与编辑、监督、项目管理、资金筹措。Wen-Wei Lin：监督、项目管理、资金筹措。Xing-Long Lyu：撰写——初稿、方法设计、实验研究、概念构思。

利益冲突声明 作者声明不存在任何可能影响本文研究结果的已知财务利益或个人关系。

Ke Gui | Tiexiang Li | Wen-Wei Lin | Xing-Long Lyu 东南大学数学学院及邵忠尧中心，中国南京211189