基于PWE法的二维光子晶体能带计算与禁带分析工具

资 源 简 介

本项目基于电磁场理论中的平面波展开法（Plane Wave Expansion Method, PWE），实现对二维光子晶体光子能带结构的精确计算与禁带分析。程序首先根据用户设定的几何参数（如晶格常数、介质柱或空气孔半径）和物理参数（背景介电常数、散射体介电常数）构建周期性介电结构模型，支持正方晶格和三角晶格等多种布拉菲晶格形式。核心算法通过将麦克斯韦方程组转化为倒易空间中的代数特征值问题，在第一布里渊区的不可约边界上进行波矢量扫描，求解出对应的本征频率。功能涵盖了TE模式（电场平行于横截面）和TM模式（磁场平行于横截面）的独立计算，能够绘制出清晰的ω-k色散关系曲线图（能带图）。系统具备自动检测光子禁带的功能，能够识别并输出带隙的起始频率、终止频率及绝对带宽，并计算相对带宽（Gap-to-Midgap Ratio）以评估禁带性能。该工具不仅适用于理论教学演示周期性结构对光波的调制作用，还能辅助科研人员设计具有特定频率禁带的光学器件，如光子晶体波导、反射镜及高Q值微腔，体现了通过禁带调控实现宽频带光控制的应用价值。

详 情 说 明

二维光子晶体能带结构与禁带特性仿真分析系统

项目简介

本项目是一个基于MATLAB开发的仿真与分析工具，旨在利用平面波展开法（Plane Wave Expansion Method, PWE）对二维光子晶体的光子能带结构进行精确计算。系统能够模拟电磁波在周期性介电结构中的传播特性，求解麦克斯韦方程组的特征值问题，绘制出TE（横电）和TM（横磁）模态的色散关系曲线，并具备自动识别光子禁带（Photonic Bandgap）及其带宽特性的功能。该工具适用于光子晶体波导设计、微腔研究以及光物理教学演示。

核心功能与特性

• 多晶格与结构支持：支持正方晶格（Square）和三角晶格（Triangular）两种布拉菲晶格形式，可配置介质柱（Rod）或空气孔（Hole）型散射体结构。
• 双模态独立求解：分别构建TE和TM偏振模态的本征方程，并行计算对应的本征频率。

* TM模态：基于标量波动方程求解，适用于电场平行于介质柱轴向的场景。 * TE模态：基于矢量波动方程求解，处理磁场平行于介质柱轴向的场景，采用优化的Ho法（逆介电常数矩阵）以提高收敛精度。

• 自动化禁带分析：系统自动扫描能带数据，识别带隙的起始频率和终止频率，计算绝对带宽以及相对带宽（Gap-to-Midgap Ratio），并能检测TE与TM模态的重叠区域（完全光子禁带）。
• 倒易空间扫描：沿第一布里渊区的不可约边界（例如 Gamma-M-K-Gamma）进行高分辨率波矢量扫描。
• 参数化仿真配置：用户可灵活调整平面波展开阶数（精度控制）、晶格常数、填充率（半径）以及背景/散射体材料的介质常数。

系统要求

• 运行环境：MATLAB R2018b 或更高版本。
• 工具箱依赖：需包含基础数学计算模块（用于矩阵运算和特征值求解）。

主程序实现逻辑与算法分析

主程序 main 是整个系统的核心入口，其执行流程严格遵循PWE算法的物理逻辑，具体实现细节如下：

1. 仿真参数初始化

程序首先定义全局配置结构体 cfg，涵盖以下关键参数：

• 几何参数：设定晶格类型（如Triangular）、散射体类型（如Hole）、晶格常数 a 及散射体半径 r。代码中默认配置展示了一个典型的产生完全带隙的三角晶格空气孔结构（r=0.48a）。
• 材料参数：分别定义背景介电常数 eps_bg 和散射体介电常数 eps_inc。例如模拟硅基底（epsilon=11.56）上的空气孔。
• 算法控制：设置平面波展开阶数 PWC，该参数决定了基底函数的数量（矩阵维数为 $(2P+1)^2$）；设定布里渊区路径的分辨率 k_res 和求解的能带数量 num_bands。

2. 预处理与代数模型构建

在进入循环计算前，程序主要完成倒易空间和材料矩阵的构建：

• 几何计算：根据晶格配置，计算正格矢量、倒格矢量以及不可约布里渊区的关键对称点（k_nodes）。
• 平面波网格生成：生成倒格矢量 $G$ 的网格集合，用于展开电磁场。
• 介电常数傅里叶变换：

* 构建介电常数的卷积矩阵 eps_matrix，用于TM模计算。 * 计算逆介电常数矩阵 eta_matrix（即 $epsilon^{-1}$ 的傅里叶展开矩阵），这是TE模态计算的关键，用于处理电位移矢量 $D$ 的连续性条件，能显著提高级数展开的收敛速度。

3. 能带结构求解（核心算法）

程序沿生成的波矢量路径 k_path 逐点进行本征值求解：

#### TM 模态计算 TM模的波动方程被转化为广义特征值问题： $$ A cdot x = lambda cdot B cdot x $$

• 左端项 (LHS)：构建对角矩阵，对角元为 $|k+G|^2$，代表自由空间的动能项。
• 右端项 (RHS)：直接使用介电常数矩阵 eps_matrix。
• 求解策略：利用 MATLAB 的 eigs 函数求解广义特征值问题，提取最小模值的 $N$ 个特征值，并转化为归一化频率 $omega a / 2pi c$。

#### TE 模态计算 TE模涉及矢量场运算，方程形式更为复杂：

• 矢量点积：计算 $(k+G_i) cdot (k+G_j)$ 的点积矩阵，体现了磁场矢量的方向性。
• 矩阵构造：将点积矩阵与逆介电常数矩阵 eta_matrix 进行元素级对应的矩阵乘法（或矩阵卷积，视具体实现维度而定，代码中体现为构建系数矩阵）。
• 求解策略：构建标准特征值问题矩阵，同样利用 eigs 求解最小实部特征值。

4. 禁带特性分析

能带计算完成后，程序调用分析模块执行以下操作：

• 逐个分析TM和TE的能带数据，寻找能带索引 $n$ 的最大值与 $n+1$ 的最小值之间的非重叠区域。
• 计算禁带宽度和中心频率，得出相对带宽比率。
• 通过求交集算法（intersect_gaps），自动判定是否存在使得光在任意偏振态下均无法传播的“完全光子禁带”。

5. 结果可视化

利用图形窗口输出直观的仿真结果：

• 双子图布局：左侧显示TM模态，右侧显示TE模态。
• 色散曲线：以不同颜色（TM红/TE蓝）绘制 $omega-k$ 关系图。
• 禁带标记：在带隙频率范围内绘制填充色块（如浅粉色），直观展示光被截止的频率区间。
• 坐标轴：横坐标标注高对称点（Gamma, M, K），纵坐标为归一化频率。

使用方法

1. 确保MATLAB当前路径包含 main.m 及相关辅助函数文件。
2. 打开 main.m，根据需求在“1. 仿真参数设置”区域修改晶格类型、半径或材料参数。
3. 直接运行脚本。
4. 观察控制台输出的进度信息、禁带数据统计，以及弹出的能带结构图。