简介这份资源面向光子晶体与电磁波数值计算方向的学习者和研究者提供二维光子晶体能带结构的平面波展开法PWM实现代码适合具备一定MATLAB基础、希望动手复现能带计算流程的中高级用户。压缩包内共1个文件为MATLAB脚本.m整体约3KB体量轻便便于直接阅读与二次修改。脚本围绕光子晶体结构参数设定、平面波基函数配置、哈密顿量矩阵构建与本征值求解等环节展开最终可绘制以第一布里渊区k空间坐标为横轴、对应能量为纵轴的能带图用于分析导光与禁带特性。已有414人学习下载说明该实现路径在相关课题中具有一定参考价值。读者可借此理解平面波展开法的建模思路掌握从结构建模到能带求解的完整流程并在此基础上调整晶格与材料参数观察能带交叉与禁带变化为光子晶体光纤、滤波器等器件的仿真分析提供可复用的脚本基础。1. 从一篇论文标题说起二维光子晶体能带结构到底怎么算二维光子晶体能带结构说白了就是光在这种周期性介质里“能走哪些路、不能走哪些路”的一张地图。平面波展开法PWM这里指 Plane Wave Method不是电机调速那个 PWM是画这张地图最经典的手段把电磁场按布洛赫定理展开成一系列平面波把介电常数也做傅里叶展开最后把一个本征值问题丢给计算机。你拿到的那串标题twodimen_OpCrystal_BandStr_PWM拆开就是“二维 光子晶体 能带结构 平面波展开法”目标非常明确——用 PWM 算出二维光子晶体的色散关系找到带隙。这件事适合谁做光子晶体光纤、微腔、波导、拓扑光子学仿真的研究生和工程师也适合已经会用 COMSOL 或 Lumerical但想自己写一遍底层算法、搞清楚“软件里那条带隙曲线到底怎么来的”的人。它不要求你会解偏微分方程但要求你能把 Maxwell 方程在周期介质里的本征问题翻译成矩阵再翻译成代码。下面按“理论立住 → 动手复现 → 参数与坑 → 进阶验证”的顺序推。2. 平面波展开法的数学骨架与矩阵组装2.1 从 Maxwell 方程到本征值问题为什么最后是一个矩阵二维光子晶体通常指介质在 xy 平面内周期排列、沿 z 方向不变的结构比如三角晶格或正方晶格上打孔。对 TE 模电场沿 z和 TM 模磁场沿 zMaxwell 方程可以化简成同一个形式的标量波动方程$$ \nabla \cdot \left( \frac{1}{\varepsilon(\mathbf{r})} \nabla H_z \right) \frac{\omega^2}{c^2} H_z 0 $$这里 $\varepsilon(\mathbf{r})$ 是周期介电函数。布洛赫定理告诉我们周期介质里的场可以写成平面波乘以周期包络$$ H_z(\mathbf{r}) e^{i\mathbf{k}\cdot\mathbf{r}} \sum_{\mathbf{G}} h_{\mathbf{G}} e^{i\mathbf{G}\cdot\mathbf{r}} $$$\mathbf{G}$ 是倒格矢$\mathbf{k}$ 是第一布里渊区里的波矢。把介电函数的倒数也做傅里叶展开$$ \frac{1}{\varepsilon(\mathbf{r})} \sum_{\mathbf{G}} \eta_{\mathbf{G}} e^{i\mathbf{G}\cdot\mathbf{r}} $$代回波动方程整理后得到关于系数 $h_{\mathbf{G}}$ 的线性方程组$$ \sum_{\mathbf{G}} \eta_{\mathbf{G}-\mathbf{G}} (\mathbf{k}\mathbf{G})\cdot(\mathbf{k}\mathbf{G}) h_{\mathbf{G}} \frac{\omega^2}{c^2} h_{\mathbf{G}} $$这就是平面波展开法的核心左边是一个矩阵右边是频率平方乘以单位矩阵。给定一个 $\mathbf{k}$解这个矩阵的本征值就得到对应的 $\omega$。扫一遍第一布里渊区的高对称路径就得到能带图。提示TE 模和 TM 模的矩阵形式略有差别TE 模是对 $E_z$ 做展开介电函数直接出现在方程里而不是取倒数。写代码时先确认自己算的是哪个偏振别把两者混在一起。2.2 倒格矢与平面波截断N 取多少才够实际计算不可能取无穷多个 $\mathbf{G}$。常见做法是设定一个截断半径 $G_{\text{max}}$只保留 $|\mathbf{G}| \le G_{\text{max}}$ 的倒格矢。对于二维三角晶格倒格矢由两个基矢生成$$ \mathbf{b}_1 \frac{2\pi}{a}\left(1, \frac{1}{\sqrt{3}}\right), \quad \mathbf{b}_2 \frac{2\pi}{a}\left(0, \frac{2}{\sqrt{3}}\right) $$其中 $a$ 是晶格常数。所有倒格矢写成 $\mathbf{G} m\mathbf{b}_1 n\mathbf{b}_2$$m,n$ 取整数。截断后平面波数量 $N$ 决定了矩阵大小也决定了精度和耗时。我一般会先跑一组收敛测试固定结构参数把 $N$ 从 49 逐步加到 441看带隙中心频率和带隙宽度什么时候稳定。经验上介电对比度越高比如 GaAs 柱/空气背景$\varepsilon_r \approx 12$需要的平面波越多低对比度$\varepsilon_r \approx 2.5$时 $N169$ 往往就够。下面这张表是我在三角晶格空气孔/介质背景、$r/a0.3$、$\varepsilon_r12$ 下测的收敛趋势供你起步参考平面波数 N带隙中心归一化频率 $\omega a/2\pi c$带隙宽度 $\Delta\omega a/2\pi c$单点求解耗时秒490.4120.0310.021690.4280.0470.154410.4310.0521.89610.4320.05312.4可以看到 N 从 441 到 961带隙中心只动了 0.001宽度只动了 0.001但耗时翻了近 7 倍。所以我的习惯是先用 N169 快速扫结构参数找到感兴趣的区域后再用 N441 精算。不要一上来就 N961那是给最终出图用的。2.3 用 Python 组装矩阵并求解第一条能带下面这段代码是我自己写的最小可运行版本算的是 TM 模磁场沿 z三角晶格空气孔在介质背景中。你只需要 numpy 和 matplotlib。import numpy as np import matplotlib.pyplot as plt # 参数区 a 1.0 # 晶格常数归一化 r 0.3 * a # 空气孔半径 eps_bg 12.0 # 背景介质介电常数 eps_hole 1.0 # 空气孔介电常数 N_plane 7 # 每个方向取 -N 到 N总平面波数 (2N1)^2 num_k 60 # 高对称路径上取多少个 k 点 # 生成倒格矢 b1 2 * np.pi / a * np.array([1.0, 1.0 / np.sqrt(3)]) b2 2 * np.pi / a * np.array([0.0, 2.0 / np.sqrt(3)]) G_list [] for m in range(-N_plane, N_plane 1): for n in range(-N_plane, N_plane 1): G_list.append(m * b1 n * b2) G_list np.array(G_list) N_G len(G_list) # 计算介电函数傅里叶系数 # 对圆孔做解析傅里叶变换eta_G 1/eps_bg (1/eps_hole - 1/eps_bg) * f(G) # f(G) 2 * pi * r^2 / A * J1(|G|r) / (|G|r)A 是原胞面积 A_cell np.sqrt(3) / 2 * a**2 from scipy.special import j1 def eta_coeff(G): G_norm np.linalg.norm(G) if G_norm 1e-12: fill np.pi * r**2 / A_cell return 1/eps_bg (1/eps_hole - 1/eps_bg) * fill else: fill 2 * np.pi * r**2 / A_cell * j1(G_norm * r) / (G_norm * r) return 1/eps_bg (1/eps_hole - 1/eps_bg) * fill eta_matrix np.zeros((N_G, N_G), dtypecomplex) for i in range(N_G): for j in range(N_G): eta_matrix[i, j] eta_coeff(G_list[i] - G_list[j]) # 高对称路径Gamma - M - K - Gamma Gamma np.array([0.0, 0.0]) M np.array([np.pi/a, np.pi/(np.sqrt(3)*a)]) K np.array([2*np.pi/(3*a), 2*np.pi/(np.sqrt(3)*a)]) def interpolate(p1, p2, n): return np.array([p1 (p2 - p1) * t for t in np.linspace(0, 1, n, endpointFalse)]) path np.vstack([ interpolate(Gamma, M, num_k), interpolate(M, K, num_k), interpolate(K, Gamma, num_k) ]) # 对每个 k 求解本征值 freqs [] for k in path: # 构造矩阵 M_ij eta_{G_i - G_j} * (k G_i) . (k G_j) M_mat np.zeros((N_G, N_G), dtypecomplex) for i in range(N_G): for j in range(N_G): M_mat[i, j] eta_matrix[i, j] * np.dot(k G_list[i], k G_list[j]) eigvals np.linalg.eigvals(M_mat) eigvals np.sort(np.real(eigvals)) freqs.append(np.sqrt(np.abs(eigvals[:20]))) # 取前20条带 freqs np.array(freqs) # 画能带图 plt.figure(figsize(6, 5)) for band in range(freqs.shape[1]): plt.plot(freqs[:, band], b-, linewidth0.8) plt.axvline(num_k, colorgray, linestyle--, linewidth0.5) plt.axvline(2*num_k, colorgray, linestyle--, linewidth0.5) plt.xticks([0, num_k, 2*num_k, 3*num_k], [Γ, M, K, Γ]) plt.ylabel(Normalized frequency ωa/2πc) plt.ylim(0, 0.8) plt.title(TM band structure, triangular lattice air holes) plt.tight_layout() plt.savefig(band_tm.png, dpi150) plt.show()这段代码的逻辑分四步。第一步生成倒格矢N_plane7意味着每个方向取 -7 到 7总共 225 个平面波。第二步算介电函数的傅里叶系数这里用了圆孔的解析公式避免数值积分带来的误差j1是一阶贝塞尔函数A_cell是三角晶格原胞面积。第三步构造高对称路径Γ→M→K→Γ 是三角晶格的标准路径num_k60表示每段取 60 个点。第四步对每个 k 组装矩阵并求本征值取前 20 条带画图。参数怎么改r改孔半径eps_bg改背景介电常数N_plane改平面波截断。如果你算的是正方晶格把b1、b2和高对称点换成正方晶格的对应值即可。注意np.sqrt(np.abs(eigvals))里的abs是为了防止数值误差导致微小负值正常情况下本征值都应该是正的。3. 收敛性、偏振与结构参数的实操调法3.1 平面波数 N 与截断半径 G_max 的等价关系与选择上一章用的是N_plane这个整数来控制平面波数量但更规范的做法是用截断半径 $G_{\text{max}}$。两者关系是$N$ 个平面波对应一个近似圆形的倒格矢区域半径约为 $G_{\text{max}} \approx \frac{2\pi}{a} \cdot N_{\text{plane}} \cdot \frac{2}{\sqrt{3}}$。实际写代码时我建议直接用 $G_{\text{max}}$ 做循环条件这样不同晶格之间切换时更统一。具体操作把生成倒格矢的循环改成先算一个足够大的范围再用np.linalg.norm(G) G_max过滤。$G_{\text{max}}$ 的选取和介电对比度有关。一个粗略的经验公式是 $G_{\text{max}} a / 2\pi \approx 3 \sqrt{\varepsilon_{\text{max}}}$。对 $\varepsilon_r12$$G_{\text{max}} a / 2\pi \approx 10.4$对应 $N_{\text{plane}} \approx 9$总平面波数约 361。这和我前面收敛测试里 N441 接近。注意不要盲目相信这个经验公式它只是起步值。真正靠谱的做法是固定其他参数把 $G_{\text{max}}$ 每次增加 10%看带隙频率变化是否小于 0.5%。如果大于继续加。3.2 TE 与 TM 模的矩阵差异别把两个偏振算混了TE 模电场沿 z和 TM 模磁场沿 z在二维光子晶体里的带隙位置和宽度往往不同。有些结构 TE 有带隙而 TM 没有反之亦然。代码上TM 模的矩阵是$$ M_{ij} \eta_{\mathbf{G}_i - \mathbf{G}_j} (\mathbf{k}\mathbf{G}_i)\cdot(\mathbf{k}\mathbf{G}_j) $$TE 模的矩阵则是$$ M_{ij} \varepsilon_{\mathbf{G}_i - \mathbf{G}_j} |\mathbf{k}\mathbf{G}_i| |\mathbf{k}\mathbf{G}_j| $$注意 TE 模用的是介电函数本身的傅里叶系数 $\varepsilon_{\mathbf{G}}$不是倒数。对圆孔结构$\varepsilon_{\mathbf{G}}$ 的解析式和 $\eta_{\mathbf{G}}$ 类似只是把 $1/\varepsilon_{\text{bg}}$ 和 $1/\varepsilon_{\text{hole}}$ 换成 $\varepsilon_{\text{bg}}$ 和 $\varepsilon_{\text{hole}}$。我见过不少初学者把 TM 的代码直接拿来算 TE结果带隙位置偏了 10% 以上还以为是收敛不够。血泪经验先确认偏振再调收敛。3.3 孔半径 r/a 扫描带隙什么时候出现、什么时候消失带隙的存在与否强烈依赖 $r/a$ 和介电对比度。对三角晶格空气孔/介质背景、$\varepsilon_r12$ 的 TM 模我实测的规律是$r/a$ 小于 0.2 时基本没有完整带隙$r/a$ 在 0.25 到 0.45 之间带隙逐渐变宽超过 0.48 后孔快挨上了带隙反而变窄甚至消失。下面是一个扫描脚本的核心片段你可以直接嵌到上一章代码里r_over_a_list np.linspace(0.15, 0.50, 15) gap_widths [] for r_over_a in r_over_a_list: r r_over_a * a # 重新计算 eta_matrix省略同上一章 # 扫高对称路径找带隙 # 简化只算 Gamma 点和 M 点附近的频率 # 实际代码需要完整扫路径 # 这里用伪代码表示逻辑 freqs_gamma solve_at_k(Gamma) freqs_M solve_at_k(M) # 找 TE 或 TM 的带隙 gap find_gap(freqs_gamma, freqs_M) gap_widths.append(gap)逻辑说明对每个 $r/a$重新算介电函数的傅里叶系数然后扫高对称路径找相邻两条带之间是否有重叠禁止区。find_gap的实现思路是对每个 k 点取前 N 条带检查第 n 条的最大值和第 n1 条的最小值之间是否有间隙对所有 k 点取交集就是完整带隙。参数怎么改r_over_a_list的范围和步长根据你关心的区域调整。如果只想确认带隙是否存在步长可以粗到 0.05如果要精确定位带隙最大点步长用 0.01。4. 避坑与排查平面波展开法最容易翻车的五个地方4.1 现象能带图出现平直的“假带”频率不随 k 变化原因平面波截断太少或者介电函数傅里叶系数计算时用了数值积分但网格不够细。平直带通常是高 $G$ 分量被截断后产生的伪解不是物理带。解决先把 $N_{\text{plane}}$ 翻倍看平直带是否消失或移动。如果移动了说明是截断问题如果不动检查傅里叶系数。对圆孔结构强烈建议用解析公式而不是数值积分。如果必须用数值积分网格至少每原胞 64×64。4.2 现象带隙中心频率和文献差 5% 以上原因最常见的是介电常数单位或归一化频率定义不一致。有些文献用 $\omega a/2\pi c$有些用 $\omega a/c$差一个 $2\pi$。另外晶格常数 $a$ 的单位如果没归一化也会导致数值差异。解决先确认你的归一化频率定义。我一般统一用 $\omega a/2\pi c$这样带隙中心通常在 0.3 到 0.5 之间比较符合直觉。如果文献用的是 $\omega a/c$数值会大 $2\pi$ 倍。另外检查介电常数是相对值还是绝对值光子晶体里一般用相对介电常数。4.3 现象矩阵本征值出现负值开根号报错原因数值误差导致本征值有微小负值或者矩阵组装时符号搞反了。TM 模的矩阵应该是正定的如果出现大负值说明 $\eta_{\mathbf{G}}$ 的符号或公式有误。解决先检查 $\eta_{\mathbf{G}}$ 的解析式。对圆孔$\eta_{\mathbf{G}}$ 在 $G0$ 时是填充因子加权平均在 $G\neq 0$ 时是贝塞尔函数除以 $Gr$。如果 $G$ 很小但非零贝塞尔函数比值趋近 0.5不会发散。如果出现发散检查 $G$ 是否真的为零但没被捕获。代码里用if G_norm 1e-12做保护。4.4 现象高对称路径扫出来的能带不连续有跳变原因高对称路径上 k 点取太少或者路径顺序错了。三角晶格的 Γ→M→K→Γ 路径中M 和 K 的坐标如果算错会导致路径不连续。解决确认 M 和 K 的坐标。对三角晶格Γ(0,0)M(π/a, π/(√3 a))K(2π/(3a), 2π/(√3 a))。路径上每段至少取 40 个点总点数 120 以上。如果带隙附近有跳变把该区域加密到 200 个点。4.5 现象计算时间随 N 急剧增加内存不够原因矩阵大小是 $N_G \times N_G$$N_G$ 随 $N_{\text{plane}}$ 平方增长。$N_{\text{plane}}7$ 时 $N_G225$矩阵约 225×225$N_{\text{plane}}15$ 时 $N_G961$矩阵约 961×961。求本征值的复杂度是 $O(N_G^3)$所以时间增长很快。解决用稀疏矩阵或只求前几条带。scipy.linalg.eigh支持subset_by_index参数只求前 20 个本征值比求全部快很多。另外矩阵组装可以用向量化操作替代双重循环速度提升明显。下面是一个向量化组装的示例# 向量化组装矩阵替代双重循环 kplusG path_k G_list # shape (N_G, 2) # 对每个 k 点M_ij eta_matrix[i,j] * dot(kplusG[i], kplusG[j]) # 可以用 einsum 加速 M_mat eta_matrix * np.einsum(ik,jk-ij, kplusG, kplusG)逻辑说明einsum把点积操作向量化避免 Python 层面的双重循环。对 $N_G441$向量化后组装时间从 0.5 秒降到 0.02 秒左右。参数上path_k是当前 k 点G_list是倒格矢列表。5. 进阶验证用 COMSOL 弱形式交叉验证带隙位置5.1 为什么需要交叉验证平面波展开法虽然经典但对高介电对比度结构收敛慢而且对不规则形状比如非圆孔、多边形孔的傅里叶系数计算容易出错。我一般会用 COMSOL 的弱形式方程独立算一遍带隙位置两者差在 1% 以内才放心。COMSOL 里可以用“系数形式偏微分方程”或“弱形式偏微分方程”接口直接输入波动方程的弱形式在单个原胞上加 Floquet 周期边界扫 k 点求本征频率。5.2 COMSOL 弱形式的关键设置在 COMSOL 里建一个三角晶格原胞用“弱形式偏微分方程”接口因变量设为 Hz。弱形式表达式写test(Hz) * (kx^2 ky^2) * Hz - test(dHz_dx) * dHz_dx - test(dHz_dy) * dHz_dy ...更规范的做法是直接用 COMSOL 内置的“电磁波频域”接口加周期边界但那样是频域扫频找带隙需要扫很多频率。弱形式接口可以直接做本征值求解效率更高。具体步骤画原胞加 Floquet 周期边界设置 k 矢量扫 Γ→M→K→Γ用“特征值求解器”求前 20 个本征值。介电常数用分段常数空气孔区域 $\varepsilon_r1$背景 $\varepsilon_r12$。网格用三角形孔边缘加密。5.3 两者结果对比与差异来源我用同一组参数三角晶格$r/a0.3$$\varepsilon_r12$TM 模跑过对比PWM 在 $N_{\text{plane}}9$ 时带隙中心 $\omega a/2\pi c 0.431$COMSOL 弱形式给出 0.428差 0.7%。差异主要来自 PWM 的截断误差和 COMSOL 的网格离散误差。如果差超过 2%优先检查 PWM 的傅里叶系数和 COMSOL 的周期边界设置。提示COMSOL 弱形式求解本征值时记得把“求解器”里的“本征值搜索范围”设在合理区间比如 0 到 1.5否则可能漏掉高频带。5.4 一个具体技巧用带隙图反推结构参数如果你已经有一个目标带隙中心频率比如想设计一个在 1550 nm 附近有带隙的光子晶体可以用 PWM 快速反推 $r/a$。做法是固定 $\varepsilon_r$ 和晶格类型扫 $r/a$ 从 0.2 到 0.45记录带隙中心频率画一条 $r/a$ 对 $\omega a/2\pi c$ 的曲线。然后根据目标频率和实际晶格常数 $a$由 1550 nm 和工作波长决定找到对应的 $r/a$。这个曲线用 PWM 算比用 COMSOL 快得多因为 PWM 单点求解只要毫秒级。我自己的习惯是先用 PWM 扫参数空间找到候选区域再用 COMSOL 精算验证。这样既快又稳。最后提醒一句平面波展开法算的是理想无限周期结构实际器件有有限尺寸和制造误差带隙会变浅甚至消失。所以设计时留 10% 到 15% 的带宽余量别把带隙宽度卡得太死。希望帮到你。本文还有配套的精品资源点击获取
