简介这是一份面向凝聚态物理、电子结构与材料计算初学者的平面波展开法PWE能带计算源码用于分析周期结构如AB复合结构中电子的能量-动量关系并绘制能带曲线。包内共2个文件均为MATLAB脚本.m分别对应周期势场定义与参数设置、平面波基底构建与哈密顿量矩阵求解等核心步骤压缩包仅1KB体量极小适合直接阅读代码、按需修改参数并重新运行。已有574人学习浏览适合正在学习能带理论、需要快速上手周期结构数值计算的学生与科研人员。通过这套代码可掌握从倒格矢建立、平面波截断数选择到特征值求解与能带绘制的完整数值流程也能进一步理解PWE方法在材料导电性、磁性和光学性质预测中的实际用途。1. PWE能带计算为什么先用平面波展开而不是直接开COMSOL拿到一个周期结构——光子晶体、声子晶体、或任意介电常数周期性排布的单元——第一件事永远是问它有没有带隙带隙开在哪几个频率这个问题的标准答案不是立刻建模仿真而是先把平面波展开PWE算一遍。PWE直接在倒空间把电磁场写成布洛赫平面波叠加把偏微分方程变成标准矩阵特征值问题几分钟内就能扫完整条能带曲线。对于二维正方晶格这类规则结构几十行Python就够对于参数扫描和趋势判断它的速度和稳定性远超通用有限元。这套方法适合谁两类人一类是做超材料、波导、滤波结构设计需要在数百种几何参数里挑带隙最大的候选另一类是已经有COMSOL等工具、但想找个独立手段交叉验证色散曲线的人。PWE不是万能的它假设结构无限周期、材料无损耗但这恰好让它成为“先算趋势、再精算细节”的第一站。2. PWE的数学骨架从Bloch波到标准特征值问题2.1 倒格子与平面波展开为什么G向量决定了一切周期结构最核心的性质是介电常数满足ε(r R) ε(r)R 是任意格矢。在这个背景下电磁场模式满足Bloch定理场可以写成调幅平面波 e^{i(kG)·r} 的叠加k 是限制在第一布里渊区内的波矢G 是倒格矢。所谓平面波展开就是把真实的场分解成一套完备的平面波基组然后用这套基组去离散麦克斯韦方程组。倒格矢 G 的集合决定了计算精度。二维正方晶格晶格常数 a正格子基矢 a1 a(1,0)a2 a(0,1)倒格子基矢 b1 (2π/a)(1,0)b2 (2π/a)(0,1)任意 G n1·b1 n2·b2。截断数 N 表示 n1、n2 从 -N 取到 N总平面波数就是 (2N1)²。N5 时 121 个平面波N7 时 225 个。这个数字直接决定矩阵维度也决定了对介电常数空间调制的分辨能力。介电常数本身也要展开成傅里叶级数。关键在于方程里出现的是 ε 的逆的傅里叶系数 ε⁻¹(G)而不是对 ε(G) 取倒数。两者的傅里叶系数完全不同混用是新手最容易犯的错。对于一个圆孔/圆柱结构ε⁻¹(G) 有解析表达式这极大简化了实现。2.2 TM/TE两种偏振的矩阵形式别把两个方程搞混二维周期结构的光学模式按偏振分成两组方程长得像但矩阵结构完全不同。TM 偏振电场沿 z 方向磁场在面内满足标量亥姆霍兹方程 ∇²E_z (ω/c)² ε(r) E_z 0TE 偏振磁场沿 z 方向满足 ∇·[ε⁻¹(r) ∇H_z] (ω/c)² H_z 0。注意 TE 方程里的 ε⁻¹ 是在算子内部这导致两个方程的离散形式不一样。把 E_z(r) Σ_G A_G e^{i(kG)·r} 代入 TM 方程利用卷积定理后得到标准矩阵特征值问题Σ_{G} ε⁻¹(G−G) |kG|² A_{G} (ω/c)² A_G矩阵元是 M_TM[i,j] ε⁻¹(Gᵢ−Gⱼ)·|kGⱼ|²。TE 偏振的推导类似但向量恒等式会额外产生一个 (kGᵢ)·(kGⱼ) 的内积因子M_TE[i,j] ε⁻¹(Gᵢ−Gⱼ)·(kGᵢ)·(kGⱼ)两个方程的共同点是左边矩阵由 ε⁻¹ 的傅里叶系数控制右边特征值就是归一化频率的平方。由于 ε⁻¹(G) ε⁻¹(−G)TE 矩阵恰好是对称的可以直接用scipy.linalg.eigh加速TM 矩阵则不对称需要用通用特征值求解器。提示TM 方程里 |kG|² 只作用在列下标上矩阵不对称是数学上的自然结果不是代码写错了。2.3 为什么从倒空间切入比实空间更快PWE 的吸引力在“一次组装、全带扫描”。实空间有限元方法如 COMSOL 弱形式需要先划分网格、再逐点扫描 k 向量并求解每个 k 点都是一次完整的有限元装配而 PWE 里 k 只出现在对角线项 |kG|² 和内积 (kGᵢ)·(kGⱼ) 中扫一个 k 点只需重新计算一个一维数组并更新矩阵然后调用一次特征值求解器。二维结构一个 k 点的矩阵装配加求解大约几十毫秒扫 60 个 k 点也就几秒。代价是截断效应。N 越大越能分辨介电常数边缘的突变但矩阵维度随 N 平方增长二维结构 N10 时矩阵已是 441×441三维结构则是 (2N1)³ 维度涨得极快。所以 PWE 快速筛选选二维、低频段足够三维或高精度计算再转到有限元。3. 用Python把PWE跑起来正方晶格光子晶体能带的最小实现3.1 搭骨架晶格、倒格矢与截断G向量下面这段代码先构建倒格矢和截断 G 向量集合是整个 PWE 的地基。所有后续矩阵装配都是在这个 G 列表上做循环。import numpy as np from scipy.special import jv from scipy.linalg import eigh def build_G_vectors(N, a): 生成二维正方晶格截断后的倒格矢列表。 N: 截断数G n1*b1 n2*b2n1,n2 取 -N..N a: 晶格常数 返回形状为 (num_G, 2) 的数组 b1 2 * np.pi / a * np.array([1.0, 0.0]) b2 2 * np.pi / a * np.array([0.0, 1.0]) G_list [] for n1 in range(-N, N 1): for n2 in range(-N, N 1): G_list.append(n1 * b1 n2 * b2) return np.array(G_list) # 例子晶格常数 1截断数 5共 121 个 G 向量 a 1.0 N 5 G build_G_vectors(N, a) print(G 向量数量:, len(G))这段代码没有难点但n1 * b1 n2 * b2的逐项相加是后续所有差值计算的基础。注意 G 列表的顺序会直接影响矩阵填充的顺序但不会影响特征值结果——求解器对基函数的排列顺序不敏感。真正影响结果的是 N 的取值后面第 4 章会专门讲怎么取。3.2 圆孔模型的解析傅里叶系数为什么不用FFT也能算对于均匀背景介质中挖圆柱孔的结构ε⁻¹ 的傅里叶系数有闭式解用贝塞尔函数 J₁ 表示。这是 PWE 里最优雅的一步——不需要任何实空间网格直接精确给出任意 G 向量差的系数。def eps_inv_fourier(G, eps_b, eps_h, r, f, a): 圆孔/圆柱结构的 1/eps 傅里叶系数解析表达式。 eps_b: 背景介质相对介电常数 eps_h: 孔内介质相对介电常数 r: 圆孔半径 f: 填充率 pi * r^2 / a^2 out np.zeros(len(G), dtypecomplex) for i, g in enumerate(G): gn np.linalg.norm(g) if gn 1e-12: # G0 项体积加权平均 out[i] f / eps_h (1 - f) / eps_b else: # 圆孔形状因子2 * J1(x) / x shape 2.0 * jv(1, gn * r) / (gn * r) out[i] (1.0 / eps_h - 1.0 / eps_b) * f * shape return out def build_eps_inv_matrix(G, eps_b, eps_h, r, a): 构造 eps_inv(G_i - G_j) 的完整矩阵形状 (num_G, num_G) diff G[:, None, :] - G[None, :, :] # (num_G, num_G, 2) diff_flat diff.reshape(-1, 2) coeff_flat eps_inv_fourier(diff_flat, eps_b, eps_h, r, np.pi * r**2 / a**2, a) return coeff_flat.reshape(len(G), len(G))参数说明jv(1, x)是第一类一阶贝塞尔函数shape 2*J1(x)/x是二维圆盘傅里叶变换的标准形式在 x→0 时趋向 1因此 G0 单独处理。build_eps_inv_matrix用广播把两两差一次性算出来代码短但内存随 G 数量平方增长——N7 时 225×225 没有问题N15 时 961×961 也才几 MB压力不大。这个解析公式的物理意义很清晰(1/eps_h - 1/eps_b)是介电常数反差反差越大非零 G 分量的系数越大带隙越容易打开。f是填充率控制孔的大小。这两个参数就是后面调带隙的核心旋钮。3.3 装矩阵求特征值一个函数扫出一条带有了 ε⁻¹ 矩阵剩下就是把 TM/TE 矩阵按公式装配并求解。下面的函数接受一个 k 点返回该 k 点的归一化频率数组。def solve_kpoint(k, G, eps_inv_mat, polarizationTM): 给定波矢 k求解该点的本征频率。 k: 二维波矢单位与 G 一致使用 2pi/a 的自然单位 polarization: TM 或 TE kG k G # (num_G, 2) kG2 np.sum(kG**2, axis1) # |kG|^2 if polarization TM: # M[i,j] eps_inv(Gi-Gj) * |kGj|^2 M eps_inv_mat * kG2[None, :] else: # TE: M[i,j] eps_inv(Gi-Gj) * (kGi)·(kGj) inner kG kG.T # (num_G, num_G) M eps_inv_mat * inner # 求解特征值取实部去掉负值数值噪声 eigvals np.linalg.eigvals(M).real eigvals eigvals[eigvals 1e-10] eigvals.sort() # 归一化频率omega * a / (2*pi*c) sqrt(lambda) * a / (2*pi) freqs np.sqrt(eigvals) * a / (2 * np.pi) return freqsTM 和 TE 的装配差异只在一行TM 用kG2[None, :]广播到每一列TE 用内积矩阵kG kG.T。np.linalg.eigvals不返回特征向量速度快适合只画能带图如果要后续做带跟踪或波函数分析就得改用np.linalg.eig。注意这里特征值排序是按频率大小排的。相邻 k 点之间同一物理带可能因简并或排序翻转而“跳带”后面避坑章节专门处理。3.4 扫完整条能带从 Γ 到 X 再到 M最后把 k 点路径列出来逐点求解并绘图。正方晶格第一布里渊区的高对称点是 Γ(0,0)、X(π/a,0)、M(π/a,π/a)路径通常取 Γ→X→M→Γ。import matplotlib.pyplot as plt def band_path_points(N_k30): 生成 Gamma-X-M-Gamma 路径上的 k 点序列 Gamma np.array([0.0, 0.0]) X np.array([np.pi/a, 0.0]) M np.array([np.pi/a, np.pi/a]) seg1 [Gamma (X - Gamma) * t / N_k for t in range(N_k)] seg2 [X (M - X) * t / N_k for t in range(N_k)] seg3 [M (Gamma - M) * t / N_k for t in range(N_k)] return np.array(seg1 seg2 seg3) k_path band_path_points(30) num_bands 8 # 只画前 8 条带 bands [] for k in k_path: freqs solve_kpoint(k, G, build_eps_inv_matrix(G, 12.0, 1.0, 0.3, a), TM) # 取前 num_bands 条 bands.append(freqs[:num_bands]) bands np.array(bands) # 横轴映射用累计线段长度做刻度 def k_distance(kpath): dist [0.0] for i in range(1, len(kpath)): dist.append(dist[-1] np.linalg.norm(kpath[i] - kpath[i-1])) return np.array(dist) x_axis k_distance(k_path) for band_idx in range(num_bands): plt.plot(x_axis, bands[:, band_idx], b-, lw1) plt.xlabel(k path) plt.ylabel(frequency (a/λ)) plt.xlim(0, x_axis[-1]) plt.show()参数说明eps_b12.0是硅的典型介电常数eps_h1.0是空气孔r0.3a是常见空气孔半径。填充率f π×0.3² ≈ 0.283。num_bands8只关心低频段的前几条带高频段平面波截断误差大画多了反而误导。横轴用累积距离而不是 k 点序号这样三条边的长度比例与真实布里渊区对应能带宽度可读。这段代码跑完你会看到在第 45 条带附近出现带隙——这就是 PWE 的核心输出。整个流程从构建 G 到出图不过几十行也印证了我前面说的“秒级出全带图”。4. PWE参数怎么设截断数、填充率与COMSOL弱形式对结果的影响4.1 三个必调参数N、填充率、折射率对比度PWE 的结果质量由三个参数主导平面波截断数 N、结构填充率 f、介电常数对比度 Δε。N 决定求解精度f 和 Δε 决定物理上有没有带隙。先说 N——它是最容易出问题的数值参数。N3 时矩阵 49×49能带大概形状能看但带隙边界可以偏 5% 以上N7 时 225 个平面波低频带隙误差收敛到 1% 以内再往上收益递减N10 对二维结构是性价比拐点。判断收敛的标准做法是固定其他参数把 N 翻一倍看带隙上下沿漂移多少漂移小于 1% 才认为截断收敛。填充率 f 是物理参数直接决定带隙的宽度和位置。圆孔空气孔在硅背景里f 从 0.2 到 0.45 之间通常存在一个最佳值使带隙最大——太小介电常数调制弱太大孔快连成片结构退化成孤立的介质柱阵列带隙反而关闭。扫描 f 时每步重新算一次 ε⁻¹ 矩阵PWE 的优势在这里体现几十个 f 值每个扫一遍能带总耗时不过几分钟。折射率对比度 Δε ε_b/ε_h 决定带隙能否打开。对比度接近 1 时无论怎么调 f 都没有带隙因为这是微扰极限布拉格散射太弱。一般经验是 TM 偏振需要对比度大于 5 才有可能在正方晶格中出现完整带隙TE 偏振要求更高。所以做设计时先检查材料选型——如果两种介质折射率差不够后面全是白算。4.2 用COMSOL弱形式复核PWE结果PWE 的“快”是一把双刃剑无限周期假设意味着它看不到缺陷态、看不到有限尺寸效应。所以我的工作流总是 PWE 先扫、COMSOL 后验。在 COMSOL 里用弱形式 PDE 模块写色散光子晶体的特征值问题本质是把 TM 方程变成弱形式积分∫∇v·∇E_z dΩ − (ω/c)² ∫v·ε(r)E_z dΩ 0再加上布洛赫周期边界条件 k 作为参数扫描。这和 PWE 是同一个物理问题的两套离散方式正好互相印证。对比项PWECOMSOL 弱形式离散空间倒空间G 截断实空间有限元网格精度控制N510最大网格尺寸 λ/10λ/20边界条件自动满足周期需要显式加 Bloch 周期条件单个 k 点耗时毫秒级秒到分钟级缺陷/超胞不支持天然支持典型用途参数扫描、趋势筛选精算验证、缺陷态设计交叉验证时有个细节COMSOL 的特征频率解出来是实际频率 f_COMSOLPWE 解出来是归一化频率 a/λ。换算公式是 f (a/λ)·c/a把 PWE 结果的纵轴乘上 c/a 再和 COMSOL 对比单位完全一致。我第一次对比时忘了乘 c/a两条曲线差了三倍白白折腾半天。COMSOL 弱形式最大的坑在网格色散曲线对网格密度不敏感但场分布对网格很敏感带隙边界附近尤其如此。我一般会在 PWE 预测的带隙边缘频率附近做一次网格细化确认带隙上下沿不随网格变化。如果 PWE 和 COMSOL 在某个 k 点偏差超过 2%优先怀疑 PWE 截断不够而不是 COMSOL 网格不够——PWE 对大对比度收敛慢是通病。4.3 参数扫描自动化把 PWE 变成设计工具PWE 的批量能力让它天然适合做“反设计”的前端。比如要找使带隙最大的空气孔半径 r我一般写一个两层循环外层遍历 r ∈ [0.2a, 0.45a]步长 0.01a内层对每个 r 扫一遍能带记录第 4 带上沿带隙下界和第 5 带下沿带隙上界之间的宽度。整个过程用上面solve_kpoint函数直接套200 次参数扫描大概十几分钟输出一张“带隙宽度 vs r”曲线设计目标一目了然。这个脚本需要小改一处solve_kpoint每次调用都要重建 ε⁻¹ 矩阵而它只依赖结构参数r、f、ε不依赖 k。把build_eps_inv_matrix提到循环外每个 r 只建一次矩阵然后几十个 k 点共享能再快 510 倍。这也是 PWE 适合参数扫描的根本原因——矩阵一次装配全程复用。5. PWE实战避坑5个让能带曲线翻车的细节5.1 特征值排序跳变导致能带图断带现象能带曲线在某个 k 点附近突然出现一条“断带”或两条带交叉后上下位置互换曲线看着像锯齿。原因每个 k 点独立调用特征值求解器特征值按数值大小排序。在简并点附近两本征值几乎相等求解器返回的顺序在相邻 k 点间可能随机翻转。这不是物理现象是数值排序的假象。解决如果只是画能带图判断带隙接受局部交叉即可带隙上下沿不受影响。如果要做带跟踪比如提取某条带的有效质量就必须保留特征向量用相邻 k 点的波函数内积判定归属。具体做法是np.linalg.eig返回特征向量矩阵计算当前 k 点某带的波函数与上一 k 点所有带的投影取投影最大的作为同一带。这个逻辑在简并点仍然不完美需要辅以对称性分类。5.2 TE/TM方程搞混导致带结构完全错误现象算出的能带和文献对照第一带斜率不对带隙位置也完全不同。原因TE/TM 的矩阵元就差一个因子——TM 是 |kG|²TE 是 (kGᵢ)·(kGⱼ)。两种偏振在低频段的色散关系差异很大用错方程特征值会差到 30% 以上。我见过有人把 TM 方程里的 |kG|² 写成了 |kGᵢ|²矩阵转置后特征值密集区域完全变形。解决写代码前先在纸上把方程推导一遍确认求和下标。一个简单的脚本级自检令 ε(r) 为常数 ε₀两种偏振都应该退化为 ω c|k|/√ε₀ 的直线。在代码里设eps_b eps_h如果输出的前几条带是直线且斜率 1/√ε₀方程就没错。这个测试应该写进单元测试每次改代码重跑。5.3 实空间FFT生成傅里叶系数时网格太粗现象改用 FFT 方法计算 ε⁻¹(G) 之后带隙宽度随 N 震荡怎么都不收敛。原因实空间离散网格把圆孔边界变成阶梯状。网格越粗阶梯误差越大高频傅里叶系数偏离真实值。解析圆孔公式没有这个问题但一旦换成任意形状结构如十字孔、C 形单元只能走 FFT于是阶梯误差进来了。解决FFT 网格数至少取 (2N1) 的 10 倍以上。N7 时网格至少 150×150常见做法是 256×256。另一个更稳的方案是超采样对每个实空间网格单元判断圆孔交叠面积用面积加权平均介电常数而不是用中心点判定。这一步能让收敛速度快一倍。5.4 大折射率对比度下的Gibbs振荡现象介电常数对比度大于 10比如硅/空气 12:1时增加 N 后带隙边界不但不收敛还会上下跳动。原因阶跃介电函数的傅里叶系数按 1/|G| 衰减截断后在阶跃边缘产生 Gibbs 振荡。这个振荡在实空间表现为介电常数的“振铃”直接影响高频带的准确性并通过矩阵耦合污染低频带。解决工程上两个选择。第一对傅里叶系数乘一个平滑因子 σ_N sinc(|G|π/N)等效于实空间做低通滤波能显著抑制振铃但会轻微展宽带隙只适合趋势扫描。第二也是我更推荐的避免在 PWE 里处理超高对比度转到实空间 FEM。PWE 的优势在中低对比度的快速筛选对比度大于 15 时它的收敛性已经不如网格自适应细化后的 COMSOL。5.5 单位漏除2π导致能带曲线整体偏移现象算出的归一化频率比文献值大 6.28 倍或者横轴刻度怎么都对不上。原因特征值 λ (ω/c)²ω √λ归一化频率 a/λ 需要 ωa/(2πc) √λ·a/(2π)。有人开方后直接画 √λ忘了除 2π。更隐蔽的版本是晶格常数 a 用了 1但倒格矢算成了 1/a 而不是 2π/a导致整个能带横坐标缩放错误。解决固定一个已知结论做验证——均匀介质 ε12 时第一带应该穿过 Γ 点频率为 0在 X 点频率为 a/λ √(12)·1/(2π·) 让我算一下ω c|k|/√ε_b归一化 a/λ ωa/(2πc) (a|k|)/(2π√ε_b)。X 点 |k|π/a所以 a/λ 1/(2√12) ≈ 0.144。如果算出来是这个数单位就对了。我把这条写进了自动化测试脚本每次改代码跑一遍被这个坑绊过的次数降到零。6. 能带曲线算完还能做什么群速度提取与对称性加速验证能带曲线上手以后第一件值得做的事是从数据里提取群速度 v_g dω/dk。这条信息对波导设计和慢光器件特别关键。实现方式是对能带数据做数值差分对每个 k 点取前后两个 k 点的频率差除以波矢差。需要注意差分要在带跟踪之后做否则排序跳变会把 v_g 算成负值或巨值。我一般用三阶中心差分公式提高精度代码就三行v_g[k] (f[k1] - f[k-1]) / (k_dist[k1] - k_dist[k-1])边界点用一阶向前/向后差分。算出群速度后能带平的地方就是慢光区域直线斜的地方是宽带导通区这个信息比单纯看带隙位置更能指导设计。第二个实用技巧是利用结构对称性缩减 k 点扫描范围。正方晶格的完整布里渊区是正方形但不可约布里渊区只有它的八分之一——Γ→X→M 的三角形区域。所有物理上不同的模式都包含在这个区域内其他区域通过镜像对称操作可以复制。扫 Γ→X→M 的三角形路径已经覆盖了全部信息这也是我第 3 章代码里把路径设为 Γ→X→M→Γ 的原因后半个三角形其实是前半个的镜像能带必须完全重合正好充当自检信号。如果你发现 Γ→X→M 段和 M→Γ 段不对称说明算法或 k 点路径有错。第三个进阶用法是对称性标记。对于旋转对称的圆孔结构TM 能带在 Γ 点有双重简并TE 能带在 M 点可能出现四重简并。这些简并点对带隙的开关有决定性影响——如果带隙边界恰好在高对称点带隙宽度对结构参数特别敏感微调 r 就可能关闭带隙。做参数扫描时我总是额外关注高对称点附近的带隙边界而不是只看带隙区间的中间值。我做能带计算这两年养成的习惯是任何新结构的第一次计算都先用 PWE 把 Γ→X→M 全路径扫一遍然后把带隙上下沿的 k 点位置记下来再用 COMSOL 只精算那两个 k 点附近——这样既保证了效率又防止 PWE 的无限周期假设掩盖了有限结构的边界效应。希望这个工作流对你也有用。本文还有配套的精品资源点击获取
PWE平面波展开法：用Python计算光子晶体能带全指南
简介这是一份面向凝聚态物理、电子结构与材料计算初学者的平面波展开法PWE能带计算源码用于分析周期结构如AB复合结构中电子的能量-动量关系并绘制能带曲线。包内共2个文件均为MATLAB脚本.m分别对应周期势场定义与参数设置、平面波基底构建与哈密顿量矩阵求解等核心步骤压缩包仅1KB体量极小适合直接阅读代码、按需修改参数并重新运行。已有574人学习浏览适合正在学习能带理论、需要快速上手周期结构数值计算的学生与科研人员。通过这套代码可掌握从倒格矢建立、平面波截断数选择到特征值求解与能带绘制的完整数值流程也能进一步理解PWE方法在材料导电性、磁性和光学性质预测中的实际用途。1. PWE能带计算为什么先用平面波展开而不是直接开COMSOL拿到一个周期结构——光子晶体、声子晶体、或任意介电常数周期性排布的单元——第一件事永远是问它有没有带隙带隙开在哪几个频率这个问题的标准答案不是立刻建模仿真而是先把平面波展开PWE算一遍。PWE直接在倒空间把电磁场写成布洛赫平面波叠加把偏微分方程变成标准矩阵特征值问题几分钟内就能扫完整条能带曲线。对于二维正方晶格这类规则结构几十行Python就够对于参数扫描和趋势判断它的速度和稳定性远超通用有限元。这套方法适合谁两类人一类是做超材料、波导、滤波结构设计需要在数百种几何参数里挑带隙最大的候选另一类是已经有COMSOL等工具、但想找个独立手段交叉验证色散曲线的人。PWE不是万能的它假设结构无限周期、材料无损耗但这恰好让它成为“先算趋势、再精算细节”的第一站。2. PWE的数学骨架从Bloch波到标准特征值问题2.1 倒格子与平面波展开为什么G向量决定了一切周期结构最核心的性质是介电常数满足ε(r R) ε(r)R 是任意格矢。在这个背景下电磁场模式满足Bloch定理场可以写成调幅平面波 e^{i(kG)·r} 的叠加k 是限制在第一布里渊区内的波矢G 是倒格矢。所谓平面波展开就是把真实的场分解成一套完备的平面波基组然后用这套基组去离散麦克斯韦方程组。倒格矢 G 的集合决定了计算精度。二维正方晶格晶格常数 a正格子基矢 a1 a(1,0)a2 a(0,1)倒格子基矢 b1 (2π/a)(1,0)b2 (2π/a)(0,1)任意 G n1·b1 n2·b2。截断数 N 表示 n1、n2 从 -N 取到 N总平面波数就是 (2N1)²。N5 时 121 个平面波N7 时 225 个。这个数字直接决定矩阵维度也决定了对介电常数空间调制的分辨能力。介电常数本身也要展开成傅里叶级数。关键在于方程里出现的是 ε 的逆的傅里叶系数 ε⁻¹(G)而不是对 ε(G) 取倒数。两者的傅里叶系数完全不同混用是新手最容易犯的错。对于一个圆孔/圆柱结构ε⁻¹(G) 有解析表达式这极大简化了实现。2.2 TM/TE两种偏振的矩阵形式别把两个方程搞混二维周期结构的光学模式按偏振分成两组方程长得像但矩阵结构完全不同。TM 偏振电场沿 z 方向磁场在面内满足标量亥姆霍兹方程 ∇²E_z (ω/c)² ε(r) E_z 0TE 偏振磁场沿 z 方向满足 ∇·[ε⁻¹(r) ∇H_z] (ω/c)² H_z 0。注意 TE 方程里的 ε⁻¹ 是在算子内部这导致两个方程的离散形式不一样。把 E_z(r) Σ_G A_G e^{i(kG)·r} 代入 TM 方程利用卷积定理后得到标准矩阵特征值问题Σ_{G} ε⁻¹(G−G) |kG|² A_{G} (ω/c)² A_G矩阵元是 M_TM[i,j] ε⁻¹(Gᵢ−Gⱼ)·|kGⱼ|²。TE 偏振的推导类似但向量恒等式会额外产生一个 (kGᵢ)·(kGⱼ) 的内积因子M_TE[i,j] ε⁻¹(Gᵢ−Gⱼ)·(kGᵢ)·(kGⱼ)两个方程的共同点是左边矩阵由 ε⁻¹ 的傅里叶系数控制右边特征值就是归一化频率的平方。由于 ε⁻¹(G) ε⁻¹(−G)TE 矩阵恰好是对称的可以直接用scipy.linalg.eigh加速TM 矩阵则不对称需要用通用特征值求解器。提示TM 方程里 |kG|² 只作用在列下标上矩阵不对称是数学上的自然结果不是代码写错了。2.3 为什么从倒空间切入比实空间更快PWE 的吸引力在“一次组装、全带扫描”。实空间有限元方法如 COMSOL 弱形式需要先划分网格、再逐点扫描 k 向量并求解每个 k 点都是一次完整的有限元装配而 PWE 里 k 只出现在对角线项 |kG|² 和内积 (kGᵢ)·(kGⱼ) 中扫一个 k 点只需重新计算一个一维数组并更新矩阵然后调用一次特征值求解器。二维结构一个 k 点的矩阵装配加求解大约几十毫秒扫 60 个 k 点也就几秒。代价是截断效应。N 越大越能分辨介电常数边缘的突变但矩阵维度随 N 平方增长二维结构 N10 时矩阵已是 441×441三维结构则是 (2N1)³ 维度涨得极快。所以 PWE 快速筛选选二维、低频段足够三维或高精度计算再转到有限元。3. 用Python把PWE跑起来正方晶格光子晶体能带的最小实现3.1 搭骨架晶格、倒格矢与截断G向量下面这段代码先构建倒格矢和截断 G 向量集合是整个 PWE 的地基。所有后续矩阵装配都是在这个 G 列表上做循环。import numpy as np from scipy.special import jv from scipy.linalg import eigh def build_G_vectors(N, a): 生成二维正方晶格截断后的倒格矢列表。 N: 截断数G n1*b1 n2*b2n1,n2 取 -N..N a: 晶格常数 返回形状为 (num_G, 2) 的数组 b1 2 * np.pi / a * np.array([1.0, 0.0]) b2 2 * np.pi / a * np.array([0.0, 1.0]) G_list [] for n1 in range(-N, N 1): for n2 in range(-N, N 1): G_list.append(n1 * b1 n2 * b2) return np.array(G_list) # 例子晶格常数 1截断数 5共 121 个 G 向量 a 1.0 N 5 G build_G_vectors(N, a) print(G 向量数量:, len(G))这段代码没有难点但n1 * b1 n2 * b2的逐项相加是后续所有差值计算的基础。注意 G 列表的顺序会直接影响矩阵填充的顺序但不会影响特征值结果——求解器对基函数的排列顺序不敏感。真正影响结果的是 N 的取值后面第 4 章会专门讲怎么取。3.2 圆孔模型的解析傅里叶系数为什么不用FFT也能算对于均匀背景介质中挖圆柱孔的结构ε⁻¹ 的傅里叶系数有闭式解用贝塞尔函数 J₁ 表示。这是 PWE 里最优雅的一步——不需要任何实空间网格直接精确给出任意 G 向量差的系数。def eps_inv_fourier(G, eps_b, eps_h, r, f, a): 圆孔/圆柱结构的 1/eps 傅里叶系数解析表达式。 eps_b: 背景介质相对介电常数 eps_h: 孔内介质相对介电常数 r: 圆孔半径 f: 填充率 pi * r^2 / a^2 out np.zeros(len(G), dtypecomplex) for i, g in enumerate(G): gn np.linalg.norm(g) if gn 1e-12: # G0 项体积加权平均 out[i] f / eps_h (1 - f) / eps_b else: # 圆孔形状因子2 * J1(x) / x shape 2.0 * jv(1, gn * r) / (gn * r) out[i] (1.0 / eps_h - 1.0 / eps_b) * f * shape return out def build_eps_inv_matrix(G, eps_b, eps_h, r, a): 构造 eps_inv(G_i - G_j) 的完整矩阵形状 (num_G, num_G) diff G[:, None, :] - G[None, :, :] # (num_G, num_G, 2) diff_flat diff.reshape(-1, 2) coeff_flat eps_inv_fourier(diff_flat, eps_b, eps_h, r, np.pi * r**2 / a**2, a) return coeff_flat.reshape(len(G), len(G))参数说明jv(1, x)是第一类一阶贝塞尔函数shape 2*J1(x)/x是二维圆盘傅里叶变换的标准形式在 x→0 时趋向 1因此 G0 单独处理。build_eps_inv_matrix用广播把两两差一次性算出来代码短但内存随 G 数量平方增长——N7 时 225×225 没有问题N15 时 961×961 也才几 MB压力不大。这个解析公式的物理意义很清晰(1/eps_h - 1/eps_b)是介电常数反差反差越大非零 G 分量的系数越大带隙越容易打开。f是填充率控制孔的大小。这两个参数就是后面调带隙的核心旋钮。3.3 装矩阵求特征值一个函数扫出一条带有了 ε⁻¹ 矩阵剩下就是把 TM/TE 矩阵按公式装配并求解。下面的函数接受一个 k 点返回该 k 点的归一化频率数组。def solve_kpoint(k, G, eps_inv_mat, polarizationTM): 给定波矢 k求解该点的本征频率。 k: 二维波矢单位与 G 一致使用 2pi/a 的自然单位 polarization: TM 或 TE kG k G # (num_G, 2) kG2 np.sum(kG**2, axis1) # |kG|^2 if polarization TM: # M[i,j] eps_inv(Gi-Gj) * |kGj|^2 M eps_inv_mat * kG2[None, :] else: # TE: M[i,j] eps_inv(Gi-Gj) * (kGi)·(kGj) inner kG kG.T # (num_G, num_G) M eps_inv_mat * inner # 求解特征值取实部去掉负值数值噪声 eigvals np.linalg.eigvals(M).real eigvals eigvals[eigvals 1e-10] eigvals.sort() # 归一化频率omega * a / (2*pi*c) sqrt(lambda) * a / (2*pi) freqs np.sqrt(eigvals) * a / (2 * np.pi) return freqsTM 和 TE 的装配差异只在一行TM 用kG2[None, :]广播到每一列TE 用内积矩阵kG kG.T。np.linalg.eigvals不返回特征向量速度快适合只画能带图如果要后续做带跟踪或波函数分析就得改用np.linalg.eig。注意这里特征值排序是按频率大小排的。相邻 k 点之间同一物理带可能因简并或排序翻转而“跳带”后面避坑章节专门处理。3.4 扫完整条能带从 Γ 到 X 再到 M最后把 k 点路径列出来逐点求解并绘图。正方晶格第一布里渊区的高对称点是 Γ(0,0)、X(π/a,0)、M(π/a,π/a)路径通常取 Γ→X→M→Γ。import matplotlib.pyplot as plt def band_path_points(N_k30): 生成 Gamma-X-M-Gamma 路径上的 k 点序列 Gamma np.array([0.0, 0.0]) X np.array([np.pi/a, 0.0]) M np.array([np.pi/a, np.pi/a]) seg1 [Gamma (X - Gamma) * t / N_k for t in range(N_k)] seg2 [X (M - X) * t / N_k for t in range(N_k)] seg3 [M (Gamma - M) * t / N_k for t in range(N_k)] return np.array(seg1 seg2 seg3) k_path band_path_points(30) num_bands 8 # 只画前 8 条带 bands [] for k in k_path: freqs solve_kpoint(k, G, build_eps_inv_matrix(G, 12.0, 1.0, 0.3, a), TM) # 取前 num_bands 条 bands.append(freqs[:num_bands]) bands np.array(bands) # 横轴映射用累计线段长度做刻度 def k_distance(kpath): dist [0.0] for i in range(1, len(kpath)): dist.append(dist[-1] np.linalg.norm(kpath[i] - kpath[i-1])) return np.array(dist) x_axis k_distance(k_path) for band_idx in range(num_bands): plt.plot(x_axis, bands[:, band_idx], b-, lw1) plt.xlabel(k path) plt.ylabel(frequency (a/λ)) plt.xlim(0, x_axis[-1]) plt.show()参数说明eps_b12.0是硅的典型介电常数eps_h1.0是空气孔r0.3a是常见空气孔半径。填充率f π×0.3² ≈ 0.283。num_bands8只关心低频段的前几条带高频段平面波截断误差大画多了反而误导。横轴用累积距离而不是 k 点序号这样三条边的长度比例与真实布里渊区对应能带宽度可读。这段代码跑完你会看到在第 45 条带附近出现带隙——这就是 PWE 的核心输出。整个流程从构建 G 到出图不过几十行也印证了我前面说的“秒级出全带图”。4. PWE参数怎么设截断数、填充率与COMSOL弱形式对结果的影响4.1 三个必调参数N、填充率、折射率对比度PWE 的结果质量由三个参数主导平面波截断数 N、结构填充率 f、介电常数对比度 Δε。N 决定求解精度f 和 Δε 决定物理上有没有带隙。先说 N——它是最容易出问题的数值参数。N3 时矩阵 49×49能带大概形状能看但带隙边界可以偏 5% 以上N7 时 225 个平面波低频带隙误差收敛到 1% 以内再往上收益递减N10 对二维结构是性价比拐点。判断收敛的标准做法是固定其他参数把 N 翻一倍看带隙上下沿漂移多少漂移小于 1% 才认为截断收敛。填充率 f 是物理参数直接决定带隙的宽度和位置。圆孔空气孔在硅背景里f 从 0.2 到 0.45 之间通常存在一个最佳值使带隙最大——太小介电常数调制弱太大孔快连成片结构退化成孤立的介质柱阵列带隙反而关闭。扫描 f 时每步重新算一次 ε⁻¹ 矩阵PWE 的优势在这里体现几十个 f 值每个扫一遍能带总耗时不过几分钟。折射率对比度 Δε ε_b/ε_h 决定带隙能否打开。对比度接近 1 时无论怎么调 f 都没有带隙因为这是微扰极限布拉格散射太弱。一般经验是 TM 偏振需要对比度大于 5 才有可能在正方晶格中出现完整带隙TE 偏振要求更高。所以做设计时先检查材料选型——如果两种介质折射率差不够后面全是白算。4.2 用COMSOL弱形式复核PWE结果PWE 的“快”是一把双刃剑无限周期假设意味着它看不到缺陷态、看不到有限尺寸效应。所以我的工作流总是 PWE 先扫、COMSOL 后验。在 COMSOL 里用弱形式 PDE 模块写色散光子晶体的特征值问题本质是把 TM 方程变成弱形式积分∫∇v·∇E_z dΩ − (ω/c)² ∫v·ε(r)E_z dΩ 0再加上布洛赫周期边界条件 k 作为参数扫描。这和 PWE 是同一个物理问题的两套离散方式正好互相印证。对比项PWECOMSOL 弱形式离散空间倒空间G 截断实空间有限元网格精度控制N510最大网格尺寸 λ/10λ/20边界条件自动满足周期需要显式加 Bloch 周期条件单个 k 点耗时毫秒级秒到分钟级缺陷/超胞不支持天然支持典型用途参数扫描、趋势筛选精算验证、缺陷态设计交叉验证时有个细节COMSOL 的特征频率解出来是实际频率 f_COMSOLPWE 解出来是归一化频率 a/λ。换算公式是 f (a/λ)·c/a把 PWE 结果的纵轴乘上 c/a 再和 COMSOL 对比单位完全一致。我第一次对比时忘了乘 c/a两条曲线差了三倍白白折腾半天。COMSOL 弱形式最大的坑在网格色散曲线对网格密度不敏感但场分布对网格很敏感带隙边界附近尤其如此。我一般会在 PWE 预测的带隙边缘频率附近做一次网格细化确认带隙上下沿不随网格变化。如果 PWE 和 COMSOL 在某个 k 点偏差超过 2%优先怀疑 PWE 截断不够而不是 COMSOL 网格不够——PWE 对大对比度收敛慢是通病。4.3 参数扫描自动化把 PWE 变成设计工具PWE 的批量能力让它天然适合做“反设计”的前端。比如要找使带隙最大的空气孔半径 r我一般写一个两层循环外层遍历 r ∈ [0.2a, 0.45a]步长 0.01a内层对每个 r 扫一遍能带记录第 4 带上沿带隙下界和第 5 带下沿带隙上界之间的宽度。整个过程用上面solve_kpoint函数直接套200 次参数扫描大概十几分钟输出一张“带隙宽度 vs r”曲线设计目标一目了然。这个脚本需要小改一处solve_kpoint每次调用都要重建 ε⁻¹ 矩阵而它只依赖结构参数r、f、ε不依赖 k。把build_eps_inv_matrix提到循环外每个 r 只建一次矩阵然后几十个 k 点共享能再快 510 倍。这也是 PWE 适合参数扫描的根本原因——矩阵一次装配全程复用。5. PWE实战避坑5个让能带曲线翻车的细节5.1 特征值排序跳变导致能带图断带现象能带曲线在某个 k 点附近突然出现一条“断带”或两条带交叉后上下位置互换曲线看着像锯齿。原因每个 k 点独立调用特征值求解器特征值按数值大小排序。在简并点附近两本征值几乎相等求解器返回的顺序在相邻 k 点间可能随机翻转。这不是物理现象是数值排序的假象。解决如果只是画能带图判断带隙接受局部交叉即可带隙上下沿不受影响。如果要做带跟踪比如提取某条带的有效质量就必须保留特征向量用相邻 k 点的波函数内积判定归属。具体做法是np.linalg.eig返回特征向量矩阵计算当前 k 点某带的波函数与上一 k 点所有带的投影取投影最大的作为同一带。这个逻辑在简并点仍然不完美需要辅以对称性分类。5.2 TE/TM方程搞混导致带结构完全错误现象算出的能带和文献对照第一带斜率不对带隙位置也完全不同。原因TE/TM 的矩阵元就差一个因子——TM 是 |kG|²TE 是 (kGᵢ)·(kGⱼ)。两种偏振在低频段的色散关系差异很大用错方程特征值会差到 30% 以上。我见过有人把 TM 方程里的 |kG|² 写成了 |kGᵢ|²矩阵转置后特征值密集区域完全变形。解决写代码前先在纸上把方程推导一遍确认求和下标。一个简单的脚本级自检令 ε(r) 为常数 ε₀两种偏振都应该退化为 ω c|k|/√ε₀ 的直线。在代码里设eps_b eps_h如果输出的前几条带是直线且斜率 1/√ε₀方程就没错。这个测试应该写进单元测试每次改代码重跑。5.3 实空间FFT生成傅里叶系数时网格太粗现象改用 FFT 方法计算 ε⁻¹(G) 之后带隙宽度随 N 震荡怎么都不收敛。原因实空间离散网格把圆孔边界变成阶梯状。网格越粗阶梯误差越大高频傅里叶系数偏离真实值。解析圆孔公式没有这个问题但一旦换成任意形状结构如十字孔、C 形单元只能走 FFT于是阶梯误差进来了。解决FFT 网格数至少取 (2N1) 的 10 倍以上。N7 时网格至少 150×150常见做法是 256×256。另一个更稳的方案是超采样对每个实空间网格单元判断圆孔交叠面积用面积加权平均介电常数而不是用中心点判定。这一步能让收敛速度快一倍。5.4 大折射率对比度下的Gibbs振荡现象介电常数对比度大于 10比如硅/空气 12:1时增加 N 后带隙边界不但不收敛还会上下跳动。原因阶跃介电函数的傅里叶系数按 1/|G| 衰减截断后在阶跃边缘产生 Gibbs 振荡。这个振荡在实空间表现为介电常数的“振铃”直接影响高频带的准确性并通过矩阵耦合污染低频带。解决工程上两个选择。第一对傅里叶系数乘一个平滑因子 σ_N sinc(|G|π/N)等效于实空间做低通滤波能显著抑制振铃但会轻微展宽带隙只适合趋势扫描。第二也是我更推荐的避免在 PWE 里处理超高对比度转到实空间 FEM。PWE 的优势在中低对比度的快速筛选对比度大于 15 时它的收敛性已经不如网格自适应细化后的 COMSOL。5.5 单位漏除2π导致能带曲线整体偏移现象算出的归一化频率比文献值大 6.28 倍或者横轴刻度怎么都对不上。原因特征值 λ (ω/c)²ω √λ归一化频率 a/λ 需要 ωa/(2πc) √λ·a/(2π)。有人开方后直接画 √λ忘了除 2π。更隐蔽的版本是晶格常数 a 用了 1但倒格矢算成了 1/a 而不是 2π/a导致整个能带横坐标缩放错误。解决固定一个已知结论做验证——均匀介质 ε12 时第一带应该穿过 Γ 点频率为 0在 X 点频率为 a/λ √(12)·1/(2π·) 让我算一下ω c|k|/√ε_b归一化 a/λ ωa/(2πc) (a|k|)/(2π√ε_b)。X 点 |k|π/a所以 a/λ 1/(2√12) ≈ 0.144。如果算出来是这个数单位就对了。我把这条写进了自动化测试脚本每次改代码跑一遍被这个坑绊过的次数降到零。6. 能带曲线算完还能做什么群速度提取与对称性加速验证能带曲线上手以后第一件值得做的事是从数据里提取群速度 v_g dω/dk。这条信息对波导设计和慢光器件特别关键。实现方式是对能带数据做数值差分对每个 k 点取前后两个 k 点的频率差除以波矢差。需要注意差分要在带跟踪之后做否则排序跳变会把 v_g 算成负值或巨值。我一般用三阶中心差分公式提高精度代码就三行v_g[k] (f[k1] - f[k-1]) / (k_dist[k1] - k_dist[k-1])边界点用一阶向前/向后差分。算出群速度后能带平的地方就是慢光区域直线斜的地方是宽带导通区这个信息比单纯看带隙位置更能指导设计。第二个实用技巧是利用结构对称性缩减 k 点扫描范围。正方晶格的完整布里渊区是正方形但不可约布里渊区只有它的八分之一——Γ→X→M 的三角形区域。所有物理上不同的模式都包含在这个区域内其他区域通过镜像对称操作可以复制。扫 Γ→X→M 的三角形路径已经覆盖了全部信息这也是我第 3 章代码里把路径设为 Γ→X→M→Γ 的原因后半个三角形其实是前半个的镜像能带必须完全重合正好充当自检信号。如果你发现 Γ→X→M 段和 M→Γ 段不对称说明算法或 k 点路径有错。第三个进阶用法是对称性标记。对于旋转对称的圆孔结构TM 能带在 Γ 点有双重简并TE 能带在 M 点可能出现四重简并。这些简并点对带隙的开关有决定性影响——如果带隙边界恰好在高对称点带隙宽度对结构参数特别敏感微调 r 就可能关闭带隙。做参数扫描时我总是额外关注高对称点附近的带隙边界而不是只看带隙区间的中间值。我做能带计算这两年养成的习惯是任何新结构的第一次计算都先用 PWE 把 Γ→X→M 全路径扫一遍然后把带隙上下沿的 k 点位置记下来再用 COMSOL 只精算那两个 k 点附近——这样既保证了效率又防止 PWE 的无限周期假设掩盖了有限结构的边界效应。希望这个工作流对你也有用。本文还有配套的精品资源点击获取
