基于平面波展开法的二维光子晶体能带仿真系统

资 源 简 介

该项目利用平面波展开法（Plane Wave Expansion Method）深入分析二维光子晶体在不同极化模式下的光学特性。其核心任务是求解受周期性介电常数调制的麦克斯韦方程组，通过将空间变化的介电函数和电磁场展开为傅里叶级数，将其转化为矩阵本征值问题进行数值求解。系统支持正方形晶格、三角形晶格等多种常见的二维结构，能够分别计算并绘制TE（横电波）和TM（横磁波）模式的能带结构图。该工具允许用户通过调整介质柱的半径、晶格常数、介质与背景的折射率对比度等关键参数，精确观测光子禁带（Photonic Ba

详 情 说 明

二维光子晶体能带特性仿真系统

项目介绍

本项目是一个基于平面波展开法（Plane Wave Expansion, PWE）的数值仿真工具，专门用于计算和分析二维光子晶体的能带结构。通过求解受周期性介电常数调制的麦克斯韦方程组，该系统能够精确模拟光子在周期性微结构中的传播特性。软件能够直观地展示在不同极化模式（TE和TM）下，由于布拉格散射产生的光子禁带（Photonic Bandgap），为光子晶体器件如波导、谐振腔和传感器提供理论设计依据。

功能特性

• 多晶格支持：内置正方形（Square）和三角形（Triangular）两种主流二维晶格结构的几何参数化模型。
• 全极化分析：支持TE（横电波，电场位于xy平面）和TM（横磁波，电场沿z轴）两种模式的独立计算。
• 参数化控制：用户可灵活调整晶格常数、介质柱半径、介质与背景的折射率比值。
• 自动化路径采样：自动生成不可约布里渊区高对称点路径（如正方形的Γ-X-M-Γ或三角形的Γ-M-K-Γ）。
• 禁带自动检索：系统能自动分析计算结果，识别全方位光子禁带的位置、宽度及其相对带宽。
• 可视化输出：同步输出倒格矢空间采样路径图、TM/TE能带结构图，并以阴影标示禁带区域。

系统要求

• 软件环境：MATLAB R2016a 或更高版本。
• 硬件要求：建议内存8GB以上，以支持高阶平面波展开矩阵的运算。
• 工具箱：仅需MATLAB核心组件，无需额外工具箱。

实现逻辑与功能模块说明

1. 参数初始化与几何建模

系统首先定义物理常数与结构模型。通过设置晶格常数 $a$ 和介质柱半径 $r$，计算结构填充率 $f$。在正方形晶格中，$f = pi r^2 / a^2$；在三角形晶格中，$f = 2pi r^2 / (sqrt{3} a^2)$。同时设定介电常数对比度（$epsilon_{rod}$ 与 $epsilon_{back}$），作为后续傅里叶变换的基础。

2. 倒格矢空间构建

根据选择的晶格类型，系统生成基矢 $b_1$ 和 $b_2$。通过在倒空间各维度取 $-N_{max}$ 到 $N_{max}$ 的阶数，构建双索引平面波集合（共 $(2N_{max}+1)^2$ 个倒格矢）。这些倒格矢用于将周期性变化的介电函数和空间电磁场展开为傅里叶级数。

3. 介电常数矩阵计算（Ho方法）

系统采用改进的收敛算法（Ho方法），即先计算介电常数在倒空间的傅里叶分量矩阵 $varepsilon(G-G')$，再对其求逆获得 $eta$ 矩阵。

• 对于圆形介质柱，其傅里叶变换通过一阶贝塞尔函数 $J_1$ 实现。
• 当 $G=G'$ 时，取平均介质常数。
• 当 $G neq G'$ 时，利用公式涉及矢量差的模长进行填充。

4. 采样路径生成

系统定义了不可约布里渊区的特征路径。

• 正方形晶格：路径为 $Gamma(0,0) to X(pi/a, 0) to M(pi/a, pi/a) to Gamma$。
• 三角形晶格：路径为 $Gamma(0,0) to M(pi/a, pi/sqrt{3}a) to K(4pi/3a, 0) to Gamma$。

每一段路径根据预设的采样点数进行线性插值，生成连续的 $k$ 矢量序列。

5. 本征方程求解

针对采样路径上的每个 $k$ 点，系统构造本征值矩阵：

• TM模式：求解算子为 $eta cdot text{diag}(|k+G|^2)$，直接针对电场 $E_z$ 建立矩阵。
• TE模式：求解算子涉及 $eta(G-G') cdot [(k+G) cdot (k+G')]$。该部分代码采用了矢量化加速技术，避免多重循环，显著提升了大型矩阵的运算速度。

通过调用本征值求解器获取前若干阶本征频率，并利用 $a/(2pi)$ 进行归一化处理。

6. 带隙分析与结果输出

主程序调用专用算法遍历计算得出的所有能带：

• 禁带判定：检查相邻能带之间是否存在频率重叠。若第 $n+1$ 带的最小值大于第 $n$ 带的最大值，则判定存在禁带。
• 参数统计：计算并打印每个禁带的起始频率、终止频率、绝对宽度以及相对带宽（百分比）。

关键算法细节分析

二维傅里叶变换的解析实现

系统放弃了离散快速傅里叶变换(FFT)，转而使用解析解。对于圆形介质柱，通过布里渊区内的积分导出 $2f J_1(x)/x$ 形式的系数，这保证了倒空间矩阵在高频部分的精确度和稳定性。

矩阵对角化加速

在TE模式的矩阵构建中，代码使用了外积运算 (Kx_vec * Kx_vec' + Ky_vec * Ky_vec')。这种矩阵化的表达方式充分利用了MATLAB对线性代数运算的底层优化，在处理大量平面波展开阶数（如 $N_{max}=7$ 产生的 225x225 矩阵）时，比嵌套循环快数倍。

自动化的带隙检索逻辑

find_gaps 函数实现了全局搜索，它不仅仅检查高对称点，而是基于整个 $k$ 轴采样数据。通过计算 min(bands(:, n+1)) - max(bands(:, n))，确保所发现的禁带是在整个 Brillouin Zone 路径上都不存在状态的“真禁带”。

可视化增强

仿真结果通过三个并排子图展示：左侧显示倒空间倒格矢分布与实际扫描路径；中间和右侧分别用蓝线和红线绘制TM和TE能带，并使用透明色块（FaceAlpha）覆盖禁带区域，使用户能一眼识别出光子受限的频率区间。