简介这份资源面向光子晶体与电磁波数值计算方向的学习者和研究者提供二维光子晶体能带结构的平面波展开法PWM实现代码帮助理解周期性介质结构中光子允许与禁带频率的形成机制适用于光学器件设计入门与课程实践。压缩包为zip格式共1个文件即一份MATLAB脚本.m整体约3KB体量轻便便于直接阅读与二次修改。脚本预计涵盖晶体结构参数设定、平面波基函数配置、哈密顿量矩阵构建与本征值求解并输出第一布里渊区能带图可用于分析能带交叉与禁带特性。目前已有414人学习下载说明该主题在光子晶体研究中具有一定关注度。读者可借助这份代码快速搭建二维光子晶体能带计算流程理解平面波展开法的核心步骤并在此基础上调整折射率、晶格常数等参数观察能带变化为光子晶体光纤、滤波器等器件的仿真分析提供参考。1. 二维光子晶体能带结构从平面波展开法到可复现的色散曲线做光子晶体仿真的人绕不开一个基本问题给定介电常数分布光在里面的色散关系到底长什么样。二维光子晶体能带结构就是这个问题的答案而平面波展开法PWMPlane Wave Method是求解它最经典、最容易自己动手实现的路径。它不需要商业软件授权不需要复杂网格划分核心就是把 Maxwell 方程在倒格矢空间展开成矩阵本征值问题。我第一次用 PWM 算二维正方晶格光子晶体能带时最大的感受是理论推导看起来吓人但代码量其实不到 200 行。真正花时间的是收敛性验证和参数调试。这篇内容面向想自己写代码算光子晶体能带的工程师和研究生从 Maxwell 方程出发一步步推到可运行的 Python 实现再把踩过的坑摊开讲清楚。2. 平面波展开法的数学框架从 Maxwell 方程到本征值问题2.1 为什么选 PWM 而不是 FDTD 或 FEM光子晶体能带计算有三条主流路线平面波展开法PWM、时域有限差分FDTD、有限元法FEM。FDTD 适合算透射谱和场分布但算能带需要大量时间步进再做傅里叶变换精度和效率都不理想。FEM 精度高但需要网格生成对光子晶体这种周期性结构反而增加了不必要的复杂度。PWM 的优势在于它天然利用周期性把问题转化成一个矩阵本征值问题一次对角化就能得到所有能带。PWM 的适用边界也很明确它假设介电常数是频率无关的实数所以不能直接处理色散材料比如金属在光频段的负介电常数需要额外处理。另外当介电常数对比度极高时PWM 收敛会变慢需要更多平面波。我一般会先确认材料参数是否在 PWM 的舒适区内再决定是否换方法。2.2 从 Maxwell 方程推导到矩阵形式对于二维光子晶体假设结构在 z 方向均匀光沿 xy 平面传播。TE 模电场沿 z 方向和 TM 模磁场沿 z 方向可以分开处理。以 TM 模为例磁场只有 z 分量满足$$\nabla \cdot \left( \frac{1}{\varepsilon(\mathbf{r})} \nabla H_z \right) \frac{\omega^2}{c^2} H_z 0$$由于介电常数 $\varepsilon(\mathbf{r})$ 是周期性的可以展开成傅里叶级数。磁场也按 Bloch 定理展开成平面波叠加。代入后得到$$\sum_{\mathbf{G}} \varepsilon^{-1}(\mathbf{G} - \mathbf{G}) |\mathbf{k} \mathbf{G}|^2 H_{\mathbf{G}} \frac{\omega^2}{c^2} H_{\mathbf{G}}$$这就是一个标准的本征值问题$A \mathbf{h} (\omega/c)^2 \mathbf{h}$其中矩阵 $A$ 的元素由介电常数的傅里叶系数和倒格矢决定。对每个 k 点对角化这个矩阵就能得到本征频率。TE 模的推导类似但介电常数和磁场的角色互换矩阵形式略有不同。实际写代码时我一般把 TE 和 TM 写成两个独立函数避免混淆。2.3 倒格矢截断与收敛性判断理论上平面波数量要取无穷多个实际计算必须截断。常见做法是取一个倒格矢截断半径 $G_{max}$所有满足 $|\mathbf{G}| \leq G_{max}$ 的倒格矢都纳入计算。对于二维正方晶格倒格矢是 $\mathbf{G} (2\pi/a)(m, n)$$m, n$ 为整数。截断半径的选择直接决定计算精度和耗时。我一般从 $G_{max} 4 \times (2\pi/a)$ 开始试逐步增大直到最低几条能带的频率变化小于 1%。对于介电常数对比度在 10:1 以内的结构通常 100~200 个平面波就能收敛。对比度到 100:1 以上时可能需要 500 个以上。注意收敛性验证不能只看最低带高能带收敛更慢。如果你关心第 5 条以上的能带截断半径要相应加大。3. 用 Python 实现二维正方晶格光子晶体能带计算3.1 构建介电常数的傅里叶系数介电常数分布是 PWM 的输入。对于最常见的圆柱孔结构空气孔在介质中介电常数可以写成$$\varepsilon(\mathbf{r}) \varepsilon_b (\varepsilon_a - \varepsilon_b) S(\mathbf{r})$$其中 $S(\mathbf{r})$ 是孔内的指示函数。它的傅里叶系数有解析表达式。对于半径为 $R$ 的圆柱填充比为 $f \pi R^2 / a^2$$$\varepsilon(\mathbf{G}) \begin{cases} \varepsilon_b (\varepsilon_a - \varepsilon_b) f \mathbf{G} 0 \ (\varepsilon_a - \varepsilon_b) f \frac{2 J_1(|\mathbf{G}|R)}{|\mathbf{G}|R} \mathbf{G} \neq 0 \end{cases}$$$J_1$ 是一阶贝塞尔函数。这个解析式避免了数值积分精度和速度都更好。import numpy as np from scipy.special import j1 def eps_fourier(G_vec, R, eps_a, eps_b, a): 计算圆柱孔结构的介电常数傅里叶系数 G_vec: 倒格矢数组, shape (N, 2), 单位 2*pi/a R: 孔半径 eps_a: 孔内介电常数 eps_b: 背景介电常数 a: 晶格常数 f np.pi * R**2 / a**2 # 填充比 G_mag np.sqrt(G_vec[:, 0]**2 G_vec[:, 1]**2) * 2 * np.pi / a eps_G np.zeros(len(G_vec), dtypecomplex) for i, Gm in enumerate(G_mag): if Gm 1e-12: # G 0 eps_G[i] eps_b (eps_a - eps_b) * f else: x Gm * R eps_G[i] (eps_a - eps_b) * f * 2 * j1(x) / x return eps_G这段代码的关键参数是R孔半径和eps_a、eps_b两种材料的介电常数。填充比f由R和a自动决定。注意G_vec的单位是 $2\pi/a$所以计算G_mag时要乘上2*np.pi/a转成物理波矢。3.2 组装本征值矩阵并求解有了介电常数的傅里叶系数就可以组装矩阵了。对于 TM 模矩阵元素为$$A_{\mathbf{G}, \mathbf{G}} |\mathbf{k} \mathbf{G}|^2 \cdot \varepsilon^{-1}(\mathbf{G} - \mathbf{G})$$这里需要介电常数的逆的傅里叶系数。直接对 $\varepsilon(\mathbf{G})$ 求逆再傅里叶变换是不对的正确做法是构建 Toeplitz 矩阵然后数值求逆。我一般用scipy.linalg.toeplitz构建 $\varepsilon$ 矩阵再求逆得到 $\varepsilon^{-1}$ 矩阵。from scipy.linalg import toeplitz, eigvals def compute_bands_tm(k_points, G_vec, eps_G, num_bands8): 计算 TM 模能带 k_points: k 点数组, shape (Nk, 2), 单位 2*pi/a G_vec: 倒格矢数组, shape (N, 2) eps_G: 介电常数傅里叶系数, shape (N,) num_bands: 计算的能带数 N len(G_vec) # 构建介电常数的 Toeplitz 矩阵 # 索引差对应 G - G eps_matrix np.zeros((N, N), dtypecomplex) for i in range(N): for j in range(N): dG G_vec[i] - G_vec[j] # 找到 dG 对应的 eps_G 索引 idx np.where(np.all(np.abs(G_vec - dG) 1e-10, axis1))[0] if len(idx) 0: eps_matrix[i, j] eps_G[idx[0]] eps_inv np.linalg.inv(eps_matrix) bands np.zeros((len(k_points), num_bands)) for ik, k in enumerate(k_points): # 构建 A 矩阵 A np.zeros((N, N), dtypecomplex) for i in range(N): kG k G_vec[i] kG_mag2 np.sum(kG**2) * (2 * np.pi)**2 # 转成物理单位 for j in range(N): A[i, j] kG_mag2 * eps_inv[i, j] # 求解本征值 eigenvalues eigvals(A) # 取实部并排序归一化频率 omega*a/(2*pi*c) freq np.sqrt(np.real(eigenvalues)) / (2 * np.pi) freq np.sort(freq) bands[ik, :] freq[:num_bands] return bands这段代码的逻辑是对每个 k 点构建矩阵 $A$然后求本征值。本征值的平方根就是归一化频率 $\omega a / (2\pi c)$。注意kG_mag2的计算中乘了 $(2\pi)^2$因为k和G_vec的单位都是 $2\pi/a$转成物理波矢要乘 $2\pi/a$平方后就是 $(2\pi/a)^2$。这里为了简洁把 $a$ 归一化为 1。3.3 沿高对称路径扫描 k 点二维正方晶格的第一布里渊区高对称路径是 $\Gamma \to X \to M \to \Gamma$。在倒格矢空间中$\Gamma$ 点$(0, 0)$$X$ 点$(0.5, 0)$$M$ 点$(0.5, 0.5)$沿这条路径取 50~100 个 k 点就能画出完整的能带图。def generate_k_path(n_per_segment50): 生成 Gamma-X-M-Gamma 路径上的 k 点 Gamma np.array([0.0, 0.0]) X np.array([0.5, 0.0]) M np.array([0.5, 0.5]) path [] for start, end in [(Gamma, X), (X, M), (M, Gamma)]: for t in np.linspace(0, 1, n_per_segment, endpointFalse): path.append(start t * (end - start)) path.append(Gamma) # 闭合 return np.array(path)取 150 个 k 点每个 k 点对角化一个 200×200 的矩阵在普通笔记本上大约 10~20 秒完成。如果平面波数量增加到 500时间会到几分钟。我一般先用少量平面波快速验证代码正确性再增加平面波数量做正式计算。4. 参数设置与收敛性验证哪些参数真正影响结果4.1 平面波数量与截断半径的选取平面波数量是 PWM 最重要的参数。它由截断半径 $G_{max}$ 决定。对于二维正方晶格倒格矢是 $(2\pi/a)(m, n)$$m, n$ 为整数。取 $G_{max} N_{max} \cdot (2\pi/a)$则所有满足 $m^2 n^2 \leq N_{max}^2$ 的整数对都纳入计算。$N_{max}$平面波数量近似适用场景329快速验证代码581低对比度、低能带8197常规计算12441高对比度、高能带16781收敛性验证我一般从 $N_{max} 8$ 开始然后增加到 12 和 16比较最低几条能带的频率变化。如果从 12 到 16 的变化小于 0.5%就认为收敛了。4.2 介电常数对比度对收敛速度的影响介电常数对比度越高介电常数傅里叶系数衰减越慢需要的平面波越多。对于 $\varepsilon 12$ 和 $\varepsilon 1$ 的对比常见于硅和空气$N_{max} 8$ 通常够用。但如果用 $\varepsilon 100$ 以上的材料可能需要 $N_{max} 16$ 甚至更高。一个实用的判断方法是计算介电常数傅里叶系数随 $|\mathbf{G}|$ 的衰减曲线。如果 $|\varepsilon(\mathbf{G})|$ 在 $G_{max}$ 处还没有降到峰值的 1% 以下说明截断半径不够。4.3 用带隙宽度做收敛性判据最直接的收敛性判据是看带隙宽度。光子晶体的带隙是能带计算的核心输出。如果增加平面波数量后带隙宽度的变化小于 1%就可以认为计算收敛了。def check_convergence(R, eps_a, eps_b, a, Nmax_list): 检查不同截断半径下的带隙宽度 results [] for Nmax in Nmax_list: # 生成倒格矢 G_vec [] for m in range(-Nmax, Nmax1): for n in range(-Nmax, Nmax1): if m**2 n**2 Nmax**2: G_vec.append([m, n]) G_vec np.array(G_vec, dtypefloat) eps_G eps_fourier(G_vec, R, eps_a, eps_b, a) k_path generate_k_path(30) bands compute_bands_tm(k_path, G_vec, eps_G, num_bands6) # 找第一条和第二条带之间的带隙 gap np.min(bands[:, 1]) - np.max(bands[:, 0]) results.append((Nmax, len(G_vec), gap)) return results这段代码会输出不同 $N_{max}$ 下的带隙宽度。如果 $N_{max}$ 从 8 增加到 12 时带隙变化很小说明 8 已经够用。我一般会跑三组数据做对比而不是只跑一组就下结论。5. 避坑与排查PWM 能带计算中最容易翻车的五个地方5.1 介电常数傅里叶系数在 G0 处算错现象能带整体偏移低频部分明显不对但形状看起来还像那么回事。原因$\mathbf{G} 0$ 处的傅里叶系数是介电常数的平均值不是背景介电常数。对于圆柱孔结构$\varepsilon(\mathbf{G}0) \varepsilon_b (\varepsilon_a - \varepsilon_b) f$其中 $f$ 是填充比。很多人直接写 $\varepsilon_b$导致平均介电常数错误。解决在代码中显式判断 $|\mathbf{G}| 10^{-12}$ 的情况用填充比公式计算。这个错误很隐蔽因为能带形状不会完全崩掉只是数值偏移。5.2 倒格矢索引与矩阵元素对应错误现象能带出现不合理的简并或交叉或者某些 k 点的频率明显异常。原因矩阵元素 $A_{\mathbf{G}, \mathbf{G}}$ 依赖 $\varepsilon^{-1}(\mathbf{G} - \mathbf{G})$。如果倒格矢数组的索引和矩阵行列的对应关系搞错或者 $\mathbf{G} - \mathbf{G}$ 的查找逻辑有 bug矩阵就会错位。解决写一个测试用例用均匀介质$\varepsilon_a \varepsilon_b$验证。均匀介质的能带是解析的$\omega c|\mathbf{k} \mathbf{G}| / \sqrt{\varepsilon}$。如果代码算出的结果和解析解一致说明矩阵组装没问题。5.3 本征值排序后能带对应关系混乱现象能带图出现交叉但物理上不应该交叉。原因eigvals返回的本征值顺序是任意的直接sort后不同 k 点的同一索引可能对应不同的物理能带。在能带交叉点附近排序会导致能带“交换身份”。解决对于能带图通常不需要追踪能带身份直接画所有本征值即可。但如果要计算带隙需要确保比较的是正确的能带。我一般用np.sort后直接画图带隙通过相邻能带的 min/max 判断不依赖能带索引的连续性。5.4 单位不统一导致频率归一化错误现象能带频率数值看起来合理但和文献对比差一个常数因子。原因PWM 中常用的归一化频率是 $\omega a / (2\pi c)$但有些人用 $\omega a / c$或者忘记把 $2\pi/a$ 乘进去。单位混乱是血泪经验里最常见的问题。解决在代码开头明确注释单位约定。我一般固定用 $a1$$c1$$k$ 和 $G$ 的单位是 $2\pi/a$最终频率是 $\omega a / (2\pi c)$。这样和大多数文献的横轴一致。5.5 高对称路径取点太少导致带隙误判现象带隙宽度算出来偏大或偏小和文献对不上。原因带隙的极值可能出现在高对称路径的中间位置而不是高对称点上。如果 k 点取太少可能错过极值点。解决每个路径段至少取 50 个点。如果带隙很窄小于 1%建议取 100 个点以上。我一般先用 50 个点快速看趋势确认有带隙后再用 200 个点精确定位。6. 进阶技巧用对称性加速计算与验证结果6.1 利用点群对称性减少计算量二维正方晶格的点群是 $C_{4v}$包含 8 个对称操作。利用对称性可以把第一布里渊区的不可约区域缩小到原来的 1/8。具体做法是只取不可约区域内的 k 点然后通过对称操作生成完整能带。对于 $N_{max} 12$441 个平面波的情况计算量可以减少到原来的 1/8速度提升非常明显。实现上我一般先生成不可约区域的 k 点计算完后用对称操作镜像到整个布里渊区。需要注意的是对称操作作用在 k 点上时倒格矢也要相应变换确保矩阵组装正确。6.2 用均匀介质解析解做代码验证在正式计算光子晶体之前先用均匀介质验证代码。均匀介质的能带是解析的$$\omega \frac{c|\mathbf{k} \mathbf{G}|}{\sqrt{\varepsilon}}$$如果代码对均匀介质算出的结果和解析解一致误差在机器精度内说明矩阵组装、本征值求解、频率归一化都没问题。这个验证步骤只需要几分钟但能省下大量调试时间。def validate_uniform(eps, Nmax5): 用均匀介质验证代码正确性 G_vec [] for m in range(-Nmax, Nmax1): for n in range(-Nmax, Nmax1): if m**2 n**2 Nmax**2: G_vec.append([m, n]) G_vec np.array(G_vec, dtypefloat) # 均匀介质eps_G[0] eps, 其余为 0 eps_G np.zeros(len(G_vec), dtypecomplex) eps_G[0] eps k_path generate_k_path(20) bands compute_bands_tm(k_path, G_vec, eps_G, num_bands4) # 解析解 for ik, k in enumerate(k_path[:5]): kG k G_vec kG_mag np.sqrt(np.sum(kG**2, axis1)) * 2 * np.pi analytic np.sort(kG_mag / np.sqrt(eps) / (2 * np.pi)) print(fk{k}, 数值{bands[ik,:4]}, 解析{analytic[:4]})运行这段代码如果数值解和解析解在 1e-10 以内一致就可以放心做正式计算了。6.3 带隙宽度对填充比的扫描光子晶体带隙宽度随填充比变化存在一个最优填充比。我一般会扫描填充比从 0.1 到 0.5找带隙最大的点。这个扫描用 $N_{max} 8$ 快速做找到最优区域后再用 $N_{max} 12$ 精确计算。填充比带隙宽度归一化频率计算时间$N_{max}8$0.150.03215 秒0.200.04815 秒0.250.05515 秒0.300.05115 秒0.350.04215 秒从表中可以看出填充比在 0.25 附近带隙最大。这个结论和文献一致说明代码可靠。6.4 我自己的使用习惯我现在做二维光子晶体能带计算固定流程是先用均匀介质验证代码再用 $N_{max} 8$ 快速扫描参数空间找到感兴趣的区域后用 $N_{max} 12$ 或 16 做精确计算最后用对称性加速。整个过程从零开始写代码到出结果大约半天时间。PWM 的代码一旦写对后续就是调参数和验证不需要反复改核心逻辑。最大的教训是不要跳过均匀介质验证。我曾经因为一个索引 bug算了三天才发现矩阵组装错了所有结果都要重来。如果一开始花十分钟做均匀介质验证这个问题当场就能发现。希望帮到你。本文还有配套的精品资源点击获取