简介平面波展开法PWE是计算周期性结构电子能带的常用数值方法。该资源面向物理、材料及凝聚态领域的研究者与学习者围绕 AB 型周期结构提供了一套基于 MATLAB 的能带曲线计算实现。包内共 2 个 M 文件整体仅 1 KB脚本覆盖周期势场定义、平面波基底构建、哈密顿量矩阵求解及能带绘制等关键步骤适合想快速上手 PWE 算法或进行二次开发的读者。能带形状与间隔直接决定材料的导电、绝缘与光学性质因此理解计算过程对材料设计和半导体研究很有价值。该资源已有 576 人浏览学习。通过研读这两个脚本读者可掌握从参数设定、矩阵组装到特征值提取的完整数值流程并能便捷地迁移至其他简单周期系统的能带计算中。1. 算能带曲线先分清目标PWE 不是仿真是把周期结构直接变成矩阵很多人在周期结构里算能带第一反应是开一个 FDTD 或有限元工具在实空间里把结构画出来、撒网格、扫频率、看透射峰。问题是这类方法算的是“一个有限样品在某种入射条件下的响应”不是“无限周期系统在不同 Bloch 波矢下的本征频率”。而平面波展开法PWEPlane Wave Expansion恰恰反着来它不追踪波怎么传播而是把周期结构里的场拆成一组平面波的叠加把“计算能带曲线”变成一个本征值问题直接在倒空间求解。面对光子晶体、声子晶体、超构材料这类周期结构PWE 通常是最快、最容易自洽验证的路径尤其适合你手头已有材料参数、想知道带隙位置和色散关系的场景。这篇文章从原理讲到 1D/2D 的 Python 实现再给出参数设置和踩坑记录适合要自己写代码、又不想把结构当黑匣子交给仿真软件的从业者。2. 先立住理论Bloch 定理、倒空间网格和介质参数展开2.1 Bloch 定理让“无限结构”坍缩成一个原胞周期结构的核心特点是无论原胞在空间里重复多少次场在平移一个晶格矢量之后只会多出一个相位因子。这就是 Bloch 定理在周期性介质中本征场可以写成周期函数与平面波的乘积波矢k是真正的好量子数。所以“无限大”的结构不需要真去建无限大模型。你只需要在一个原胞上定义材料分布然后把场展开成[ \psi_{\mathbf{k}}(\mathbf{r}) e^{i\mathbf{k}\cdot\mathbf{r}} \sum_{\mathbf{G}} u_{\mathbf{G}} e^{i\mathbf{G}\cdot\mathbf{r}} ]其中G是倒格矢。展开之后原本需要处理边界条件的微分方程变成了一个以G为下标的本征值问题。对每一个k解一次矩阵本征问题就能得到一组频率把k沿着布里渊区的高对称线扫过去把得到的频率点连起来就是能带曲线。这个转化是 PWE 的关键k是你主动扫描的变量G是矩阵的维度来源。周期结构不会出现在网格文件里而是体现在介质参数的傅里叶系数中。2.2 1/ε(r) 的傅里叶系数比 ε(r) 本身更值得认真对待实际写 PWE 时多数光子晶体代码并不直接用介电常数 ε(r) 的展开而是用它的倒数 1/ε(r)。原因出在电磁场方程的形式上。对非磁性、无源介质磁场满足[ \nabla \times \left[ \frac{1}{\varepsilon(\mathbf{r})} \nabla \times \mathbf{H} \right] \left( \frac{\omega}{c} \right)^2 \mathbf{H} ]这个方程左边是厄米的离散成矩阵后本征值一定是实数数值上非常好处理。如果改用电场方程出来的矩阵就不是标准厄米形式本征值容易混入虚部谱也不干净。所以最稳妥的做法是先把原胞内每个位置的 1/ε(r) 写出来再对它做傅里叶展开。对于分层结构、圆杆、方杆这类规则几何傅里叶系数有解析式直接代入即可。G0 的系数就是 1/ε(r) 的体积平均值其他 G 的系数体现了结构对比度对比度越高高频分量的贡献越大。2.3 布里渊区路径与 k 抽样算“曲线”而不是算“点”能带图的横轴不是频率而是布里渊区边界上的k路径。二维正方形晶格通常走 Γ → X → M → R三维或 Γ → X → M → Γ二维一维结构走 Γ → X 就够。这里有一个容易忽略的点能带曲线的“带隙”指的是在所有可能的k下都不存在本征频率的区间。如果你只扫 Γ 点看到的只是 Γ 点的态密度不是带隙。所以在计算前先把晶格类型对应的倒格矢和布里渊区高对称点标出来。正方形晶格的倒空间还是正方形六角晶格的倒空间是六角形路径选错算出来的“带隙”很可能是个视觉错觉。3. 用 Python 手写 1D 布拉格堆栈能带全代码与参数3.1 最小 PWE 实现一维周期分层结构的能带曲线先从最简单的 1D 情形开始一维布拉格堆栈即两种介质交替分层周期为 d。设 A 层厚度 f*d、介电常数 ε_aB 层厚度 (1-f)*d、介电常数 ε_b波沿 x 方向传播。标量波动方程展开成平面波后本征方程写成[ (kG)^2 u_G \left( \frac{\omega}{c} \right)^2 \sum_{G} \varepsilon_{G-G} u_{G} ]也就是一个广义本征值问题左边是对角矩阵右边是介电常数傅里叶系数构成的满矩阵。下面的代码可以直接跑出前几条能带。import numpy as np import scipy.linalg as la def eps_coeff_1d(m, eps_a, eps_b, f): 一维分层结构的介电常数傅里叶系数。 m 是倒格矢序数对应 G 2*pi*m / d。 结构定义0~f*d 是 A 层f*d~d 是 B 层。 m np.asarray(m, dtypefloat) c np.zeros_like(m, dtypecomplex) zero (np.abs(m) 1e-12) c[zero] f * eps_a (1 - f) * eps_b nz ~zero c[nz] (eps_a - eps_b) * np.exp(-1j*np.pi*m[nz]*f) * np.sin(np.pi*m[nz]*f) / (np.pi*m[nz]) return c def band_1d(eps_a4.0, eps_b1.0, f0.5, N30, nk200): # 周期取 d1所以倒格矢单位是 2*pi idx np.arange(-N, N1) G 2.0 * np.pi * idx nG len(idx) # 组装介电常数矩阵 B[G, G] eps_{G-G} B np.zeros((nG, nG), dtypecomplex) for p in range(nG): for q in range(nG): B[p, q] eps_coeff_1d(idx[p] - idx[q], eps_a, eps_b, f) # k 路径一维布里渊区 0 ~ pi/d ks np.linspace(0, np.pi, nk) bands [] for k in ks: A np.diag((k G)**2) # 广义本征问题A u (omega/c)^2 B u w la.eigvals(A, B) w np.real(w[np.isreal(w)]) w w[w 1e-8] bands.append(np.sort(w)[:6]) bands np.array(bands) # 归一化频率omega*a / (2*pi*c) freqs np.sqrt(bands) / (2.0*np.pi) return ks, freqs这段代码的逻辑很直白先建倒空间坐标再组装右端矩阵 B然后对每个 k 求广义本征值取前几条实数本征值作为能带。归一化频率除以 2π目的是让横轴成为常见的 ωa/(2πc) 无量纲量方便和文献对比。参数说明N 是倒格矢截断数矩阵维度是 2N1。N30 时矩阵只有 61×61本征求解几乎瞬时完成。f 是 A 层占空比ε_a 和 ε_b 是两种介质的介电常数。把 ε_a 换成 8.9、f 改成 0.5能带里会明显看到第一带隙变宽。3.2 代码逻辑和参数为什么取前几条本征值而不是全部广义本征值求出来的频率是 ω/c 的平方所以代码里取了正实数后开根号。实际结构中介电常数恒为正B 矩阵正定理论上不会有负本征值。但因为矩阵用了复数傅里叶系数数值上会有微小虚部所以先过滤掉了虚部。取前 6 条能带是习惯做法。光子晶体能带图一般关心低频段的前几条带它们决定了带隙的位置和宽度。如果 N 取得太小高序 G 的平面波没被充分表示前几条带可能还凑合但带边缘频率会有毫秒级误差。后面会专门讲收敛性判断。这里有一个值得注意的细节代码里 B 矩阵用双层循环组装N 不大时无所谓。但如果你把 N 加到 200循环就会拖慢。更高效的做法是用 scipy.linalg.toeplitz 直接构造循环矩阵既快又不容易写错。3.3 画能带图与收敛性检查不要只看曲线形状拿到 freqs 后直接 plot 就能看到经典的布拉格堆栈能带。但真正判断代码写没写对要看两点。第一Γ 点k0的频率不应该是零。一维分层结构的低频段会从某个有限频率开始如果看到 k0 处有零频模说明边界条件或者矩阵组装出了问题。第二把 N 从 10 改成 40观察第一条带隙的下边缘频率变化。变化小于 1% 就认为收敛。很多教程里只给你看最终能带图不提收敛性实际上一旦对比度很高、比如 ε_a13、ε_b1需要 50 个以上倒格矢才能让带边缘稳定。4. 扩展到 2D 周期结构正方形晶格介质柱的本征问题4.1 2D 本征问题的矩阵组装TE 和 TM 用的是两套矩阵二维光子晶体里结构沿 z 方向不变波在 xy 平面传播。这时极化方向是个绕不开的判断点电场沿 z 方向E 极化/类 TM和磁场沿 z 方向H 极化/类 TE对应的矩阵完全不一样。PWE 的标准写法是用 1/ε 的傅里叶系数 κ(G)对本征场做平面波展开后TM 类极化的矩阵是[ \sum_{\mathbf{G}} \kappa(\mathbf{G}-\mathbf{G}), (\mathbf{k}\mathbf{G})\cdot(\mathbf{k}\mathbf{G}), u_{\mathbf{G}} \frac{\omega^2}{c^2} u_{\mathbf{G}} ]TE 类极化则是[ \sum_{\mathbf{G}} \kappa(\mathbf{G}-\mathbf{G}), |\mathbf{k}\mathbf{G}|,|\mathbf{k}\mathbf{G}|, u_{\mathbf{G}} \frac{\omega^2}{c^2} u_{\mathbf{G}} ]两者只有核函数不同前者是倒格矢之间的点积后者是模长的乘积。这个差别直接影响带隙位置很多人在 2D 能带上“翻车”就是因为 TE、TM 混用了同一套矩阵。4.2 正方形晶格介质柱的傅里叶系数与求解代码二维结构的核心输入是原胞内 1/ε 的分布。最常见的是背景介质中放一个圆形介质柱柱半径 R、介电常数 ε_r背景 ε_b。圆盘的傅里叶系数有解析式用贝塞尔函数 J1 表达from scipy.special import j1 def inv_eps_fourier_2d(Gx, Gy, er, eb, R, a1.0): 正方形晶格中圆形介质柱的 1/epsilon 傅里叶系数。 Gx, Gy 是倒格矢的两个分量单位是 2*pi/a。 柱中心在原点背景是 eb柱体是 er。 q np.sqrt(Gx**2 Gy**2) f_fill np.pi * R**2 / (a*a) coeff np.zeros_like(q, dtypecomplex) g0 q 1e-12 # G 0体积平均值 coeff[g0] f_fill / er (1 - f_fill) / eb # G ! 0圆盘形状因子 nz ~g0 coeff[nz] (1/er - 1/eb) * f_fill * 2.0 * j1(q[nz]*R) / (q[nz]*R) return coeff def assemble_pwe_2d(kx, ky, Gx, Gy, er, eb, R, modetm): 组装 2D PWE 矩阵。modetm 用点积核modete 用模长核。 nG len(Gx) H np.zeros((nG, nG), dtypecomplex) for i in range(nG): kGi np.array([kx Gx[i], ky Gy[i]]) for j in range(nG): kGj np.array([kx Gx[j], ky Gy[j]]) kappa inv_eps_fourier_2d(Gx[i]-Gx[j], Gy[i]-Gy[j], er, eb, R) if mode tm: H[i, j] kappa * np.dot(kGi, kGj) else: H[i, j] kappa * np.linalg.norm(kGi) * np.linalg.norm(kGj) return H组装完成后对每个 k 点求 H 的本征值取最小几个开根号就是该点的能带频率。注意 2D 矩阵是复数厄米矩阵要使用scipy.linalg.eigh不要用eig否则虚频分支会出来捣乱。4.3 2D 能带参数与参考结果怎么知道算对了以正方形晶格空气孔或介质柱为例取晶格常数 a1介质柱半径 0.2a柱体 ε_r8.9背景 ε_b1。第一带隙通常出现在 TM 极化频率约在 0.30~0.45 c/a 之间。如果算出来第一带隙落在 0.1 附近多半是矩阵里核函数写错或者 k 路径没有归一化。倒格矢的截断直接影响 2D 结果。正方形晶格中G 的选取一般是 [-N, N] 的二维网格总共 (2N1)² 个平面波。N15 时有 961 个平面波矩阵已经接近千维本征求解还不算太吃力N30 就是 3721 维建议改用稀疏本征求解器只求最低几条带。占空比 f 也是关键参数。介质柱直径从 0.3a 变到 0.5a带隙宽度变化很大。你在调参时每改一次半径都应该重新跑一遍收敛性检查不要沿用上一组半径确定的 N。5. 算能带最容易翻车的五个坑现象、原因与解决5.1 平面波数量不够带边缘频率随 N 漂移现象把 N 从 20 改到 40能带图整体形状没变但带隙边界频率移动了 3% 以上。原因高对比度结构里1/ε 分布存在陡峭突变傅里叶系数衰减慢。平面波截断太狠边缘处的高频成分被砍掉。解决固定一个关心的 k 点比如带隙边界所在的 X 点扫描 N10, 20, 40, 80观察最低非零本征值的变化。相对变化小于 0.5% 再画最终的能带图。这个步骤也是血泪经验能带曲线看起来平滑不代表收敛了。5.2 Γ 点出现零频模式现象能带图里 Γ 点频率为 0而且不止一条。原因本征方程允许无旋无散的自由度存在。在一维标量方程里还不太容易遇到2D 或 3D 的矢量 PWE 里未约束的平面波展开会引入额外零模。解决先确认自己用的是磁场方程而不是电场方程其次在 2D 里按极化分开计算用标量场而不是矢量场展开。如果必须要做 3D 矢量 PWE需要构造横向投影算子把 (kG) 方向的纵向分量剔除。5.3 介电常数傅里叶系数符号或归一化错误现象能带图算出来异常要么频率整体偏小一个量级要么带隙位置完全对不上文献。原因最常见的是圆盘形状因子漏了占空比 f或者 G0 的系数没有加上背景介质项。另一个常见错误是把正格矢和倒格矢的单位混了周期 a 忘了归一化。解决单独打印 G0 的系数手动验证它是否等于 f/ε_r (1-f)/ε_b。再检查 G 很小时的数值应该趋近于体积平均值。这两项对了矩阵组装的基本盘就对了。5.4 TE、TM 极化用错同一套矩阵现象正方形晶格介质柱算出来的 TM 带隙和 TE 带隙一模一样。原因组装矩阵时没有区分点积核和模长核。也有人把模式名称记反了E 极化写了模长核。解决先拿一个已知结构做锚点。比如正方晶格介质柱TM 极化通常在 Γ-X 方向有带隙TE 极化的带隙往往更窄或不存在。用这个特性去验证自己的代码。5.5 高对称点简并没出现现象Γ 点或 M 点出现了本应简并的能带却轻微劈裂。原因一种原因是倒格矢截断使用了不对称的矩形区域比如 nx 取 20、ny 取 15破坏了晶格对称性。另一种原因是 k 点没有落在高对称线上采样点偏移了。解决倒格矢截断必须保持晶格自身的旋转对称性。正方形晶格就取 [-N,N]×[-N,N]不要取矩形k 路径的端点和拐点严格落在 Γ、X、M 的坐标上不要用“差不多”的浮点数。6. 验证和进阶用参考值兜底再把 PWE 用到更复杂的周期结构上能带代码写完第一件事不是画漂亮的彩图而是找一组别人已经发表过的结构做对照。一维布拉格堆栈可以用传输矩阵法交叉验证二维正方晶格介质柱可以参考光子晶体教材里的经典能带图。只对标一两个频率点是不够的要把整个带隙边界、简并点位置都比一遍确认前几条带的相对顺序都一致。PWE 的验证里还有一个实用技巧把占空比推到极限看看能不能回归到已知结果。比如让介质柱半径趋近于 0结构退化成均匀介质能带应该是直线色散让占空比等于 0 或 1也一样。这些极值情况能暴露矩阵组装里的隐藏 bug比拿着正常参数反复调更快。算完能带之后下一步往往不是继续加大 N而是算群速度。群速度可以从能带曲线直接求导得到先对离散的能带做平滑再用中心差分求 dω/dk。也可以直接把本征频率和本征矢一起输出利用本征矢计算场分布做模式分析。如果后面要用 PWE 结果去驱动有限尺寸结构的仿真建议把本征频率、本征矢、对应 k 点三个数据同时存下来做成一个轻量级的 HDF5 文件后续不管是画模式场还是做耦合分析都不用重新跑一遍 PWE。到了这一步如果你发现 N 大到矩阵求解变慢或者想算声子晶体、弹性波能带可以考虑把稠密本征求解换成稀疏求解器比如 scipy.sparse.linalg 里的 shift-invert 模式。PWE 的物理本质不变变的只是求本征值的手段。我自己写能带代码这几年最大的体会是平面波展开法看起来很“数学”但真正消耗时间的从来不是推导公式而是矩阵组装里那些坐标归一化、截断对称性、极化选错的问题。这些问题只要在代码里留下清晰的检查和锚点对比后面换任何周期结构都不会再翻车。希望帮到你。本文还有配套的精品资源点击获取
