平面波展开法MATLAB实现：二维光子晶体完全禁带判定与调参避坑指南-育师

简介：这份资源面向光子晶体光学特性研究与 MATLAB 数值计算的学习者，核心是二维平面波扩展法（PWEM）分析光子晶体完全禁带的程序集合。包内共 27 个文件，全部为 .m 脚本，压缩包约 78KB，涵盖方形、三角、蜂窝、椭圆、矩形、柱状等多类晶格结构的禁带计算入口，以及若干 test 测试脚本，便于逐模块运行与对照验证。程序围绕 Bloch 定理构建平面波基函数、组装哈密顿量矩阵并求解特征值，从而得到频谱分布与禁带宽度、中心频率等参数，可用于光子晶体光纤、微腔与激光器等器件的设计分析。已有 260 人学习，适合熟悉 MATLAB、希望快速上手 PWEM 并进一步探索结构参数敏感性的研究人员参考。

1. 平面波展开法到底在算什么：从 PBG 到完全禁带的物理直觉

做光子晶体研究的人，绕不开一个核心问题：给定介电常数分布，这个结构到底有没有完全禁带（PBG）？平面波展开法（PWEM）就是回答这个问题最经典、最直接的手段。它的思路并不复杂——把电磁场在倒格矢空间展开成一系列平面波的叠加，把麦克斯韦方程组转化成一个矩阵本征值问题，解出来的本征频率随波矢 k 的变化关系就是能带结构。如果某个频率区间内，对所有 k 方向都不存在传播模式，那就是完全禁带。

这个标题对应的场景很具体：用 MATLAB 实现二维平面波法，分析光子晶体的能带结构，判断是否存在完全禁带。适合谁？做光子晶体、光子集成、微波光子学方向的研究生和工程师，尤其是需要自己写代码而不是依赖商业软件的人。因为只有自己实现一遍，你才能真正理解收敛性、截断阶数、极化模式这些参数到底怎么影响结果。很多论文里给的带隙宽度对不上，十有八九是截断阶数不够或者介电常数傅里叶展开出了问题。

2. 平面波展开法的数学框架与 MATLAB 实现路径

2.1 从麦克斯韦方程组到本征值问题

二维光子晶体的平面波展开法，核心推导起点是频域麦克斯韦方程组。假设结构在 z 方向均匀、电磁波在 xy 平面传播，对于 TE 极化（E 沿 z 方向），电场满足：

∇ × (1/ε(r) ∇ × E) = (ω/c)² E

把 1/ε(r) 和 E(r) 都按倒格矢 G 展开成傅里叶级数，代入后整理，得到标准本征值方程：

Σ_G' |k+G| |k+G'| ε⁻¹(G-G') E(G') = (ω/c)² E(G)

这里 ε⁻¹(G-G') 是介电常数倒数的傅里叶系数，k 是第一布里渊区内的波矢，G 是倒格矢。对于 TM 极化（H 沿 z 方向），方程形式略有不同，但思路一致。

关键点在于：ε⁻¹(G) 的傅里叶系数不能直接用 ε(G) 的倒数。对于圆形介质柱光子晶体，ε⁻¹(G) 有解析表达式，涉及贝塞尔函数和填充比。这是很多人第一次写代码时翻车的地方——直接对 ε(r) 做 FFT 然后取倒数，结果带隙位置偏了一大截。

2.2 用 MATLAB 构建倒格矢与介电常数矩阵

下面是我常用的倒格矢生成和介电常数傅里叶系数计算的核心代码。以正方晶格、圆形介质柱为例：

% 参数设置 a = 1; % 晶格常数，归一化 r = 0.2*a; % 介质柱半径 eps_bg = 1; % 背景介电常数（空气） eps_rod = 12; % 介质柱介电常数 N = 15; % 截断阶数，即每个方向取 -N 到 N % 生成倒格矢 % 正方晶格倒格矢：G = m*b1 + n*b2 % b1 = (2*pi/a, 0), b2 = (0, 2*pi/a) m = -N:N; n = -N:N; [M, Nn] = meshgrid(m, n); Gx = (2*pi/a) * M(:); Gy = (2*pi/a) * Nn(:); G = [Gx, Gy]; numG = size(G, 1); % 总平面波数 % 计算介电常数倒数的傅里叶系数 % 对于圆形柱，有解析公式 f = pi*r^2/a^2; % 填充比 eps_inv = zeros(numG, numG); for i = 1:numG for j = 1:numG dG = G(i,:) - G(j,:); dG_mag = norm(dG); if dG_mag < 1e-10 % G=0 分量 eps_inv(i,j) = f/eps_rod + (1-f)/eps_bg; else % G≠0 分量，圆形柱解析式 eps_inv(i,j) = (1/eps_rod - 1/eps_bg) * ... f * 2 * besselj(1, dG_mag*r) / (dG_mag*r); end end end

这段代码的逻辑说明：倒格矢的生成决定了平面波基组的完备性，N 越大基组越完备，但矩阵维度是 (2N+1)²，N=15 时矩阵已经 961×961，再大计算量急剧上升。介电常数傅里叶系数的计算是精度关键，圆形柱的解析式来自贝塞尔函数的积分性质，比数值 FFT 稳定得多。参数 r/a 即填充比直接决定带隙宽度，一般 r/a 在 0.15 到 0.3 之间容易出现完全禁带。

2.3 扫描布里渊区得到能带结构

有了介电常数矩阵，接下来对第一布里渊区高对称点路径上的每个 k 点构造本征值问题并求解：

% 高对称点路径：Γ -> X -> M -> Γ（正方晶格） Gamma = [0, 0]; X = [pi/a, 0]; M = [pi/a, pi/a]; num_k = 40; % 每段取40个点 % 生成 k 点路径 k_path = []; k_path = [k_path; linspace(Gamma(1), X(1), num_k)', linspace(Gamma(2), X(2), num_k)']; k_path = [k_path; linspace(X(1), M(1), num_k)', linspace(X(2), M(2), num_k)']; k_path = [k_path; linspace(M(1), Gamma(1), num_k)', linspace(M(2), Gamma(2), num_k)']; num_k_total = size(k_path, 1); freqs = zeros(num_k_total, numG); for ik = 1:num_k_total k = k_path(ik, :); % 构造本征值矩阵 A = zeros(numG, numG); for i = 1:numG for j = 1:numG kG_i = k + G(i,:); kG_j = k + G(j,:); A(i,j) = norm(kG_i) * norm(kG_j) * eps_inv(i,j); end end % 求解本征值 eigenvalues = eig(A); freqs(ik, :) = sqrt(abs(eigenvalues)) * a / (2*pi); end % 绘图 figure; plot(1:num_k_total, freqs, 'b.', 'MarkerSize', 2); xlabel('波矢 k 路径'); ylabel('归一化频率 (\omega a / 2\pi c)'); title('二维光子晶体能带结构（TE 极化）'); % 标注高对称点 xticks([1, num_k, 2*num_k, 3*num_k]); xticklabels({'\Gamma', 'X', 'M', '\Gamma'}); ylim([0, 1.2]); grid on;

逻辑说明：对每个 k 点，矩阵 A 的元素由 |k+G| 和介电常数傅里叶系数共同决定。求解 eig 得到本征值后开方得到频率。绘图时把所有 k 点的本征频率叠加，就能看到能带的分布。如果某段频率范围内没有任何本征值曲线穿过，那就是禁带。参数 num_k 决定路径采样密度，40 一般够用，太少会漏掉带隙的精细结构。

3. 完全禁带判定的三个关键参数与调参策略

3.1 截断阶数 N 怎么选才不翻车

截断阶数 N 是平面波展开法里最玄学的参数。N 太小，基组不完备，算出来的带隙位置和宽度都不可信；N 太大，矩阵维度爆炸，计算时间不可接受。我的血泪经验是：先做收敛性测试。固定结构参数，从 N=5 开始逐步增加，观察带隙宽度和中心频率的变化。当 N 增加到某个值后，结果变化小于 1%，就认为收敛了。

对于介电常数对比度在 10 左右的光子晶体，N=12 到 15 通常够用。但如果介质柱和背景的介电常数比超过 20，收敛会慢很多，可能需要 N=20 以上。注意，N 增大时矩阵维度是 (2N+1)²，N=15 是 961 维，N=20 是 1681 维，eig 求解时间差了好几倍。

提示：每次改结构参数后都要重新做收敛性测试，不能拿之前的最优 N 直接套用。

3.2 极化模式的选择与带隙差异

TE 极化和 TM 极化的能带结构可以完全不同。对于圆形介质柱在空气中的正方晶格，TE 极化（电场沿柱轴方向）通常更容易出现完全禁带，而 TM 极化往往没有完全禁带或者带隙很窄。这不是代码写错了，是物理本质。

如果你只算了一种极化就下结论说“没有完全禁带”，那是不严谨的。正确做法是两种极化都算，分别判断。有些结构在 TE 下有大带隙，在 TM 下完全没有，这种结构对偏振敏感，实际应用中要特别注意。

代码上只需要改本征值矩阵的构造方式：TE 极化用 ε⁻¹ 的傅里叶系数，TM 极化用 ε 的傅里叶系数直接构造。别搞混了，否则结果完全对不上。

3.3 填充比 r/a 对带隙宽度的影响规律

填充比是设计光子晶体时最直接的可调参数。对于正方晶格圆形柱结构，TE 极化下，完全禁带宽度随 r/a 的变化呈先增后减的趋势。r/a 太小，介电常数调制弱，带隙打不开；r/a 太大，介质柱之间耦合过强，能带折叠效应导致带隙消失。

我一般会先扫一遍 r/a 从 0.1 到 0.4，步长 0.02，快速找到带隙最大的填充比区间，再在这个区间内精细扫描。下面是一个扫参的代码框架：

r_over_a = 0.1:0.02:0.4; gap_width = zeros(size(r_over_a)); gap_center = zeros(size(r_over_a)); for ir = 1:length(r_over_a) r = r_over_a(ir) * a; % 重新计算 eps_inv 和能带 % ...（省略重复代码，同2.2和2.3节） % 提取带隙 % 假设 freqs 已经算好 % 找禁带：对所有 k 点，统计频率间隙 % 简化处理：取所有本征频率排序后找最大间隙 all_freqs = sort(freqs(:)); diffs = diff(all_freqs); [max_gap, idx] = max(diffs); gap_width(ir) = max_gap; gap_center(ir) = (all_freqs(idx) + all_freqs(idx+1)) / 2; end figure; yyaxis left; plot(r_over_a, gap_width, 'b-o'); ylabel('带隙宽度 (\omega a / 2\pi c)'); yyaxis right; plot(r_over_a, gap_center, 'r-s'); ylabel('带隙中心频率 (\omega a / 2\pi c)'); xlabel('填充比 r/a'); title('带隙宽度与中心频率随填充比的变化'); grid on;

这段代码的逻辑是：对每个填充比，重新计算能带，然后从所有本征频率中找最大间隙作为带隙宽度的估计。注意这种简化方法只适用于全局带隙的粗略估计，严格判断完全禁带需要对每个 k 点分别检查频率间隙是否在所有 k 方向上都存在。参数步长 0.02 是经验值，太粗会漏掉极值点，太细计算量翻倍。

4. 避坑与排查：平面波展开法最常见的五个翻车现场

4.1 带隙位置整体偏移：介电常数傅里叶系数算错了

现象：算出来的带隙中心频率和文献对不上，偏移量超过 10%。

原因：最常见的是直接对 ε(r) 做数值 FFT 然后取倒数得到 ε⁻¹(G)，这是错误的。正确的做法是对 1/ε(r) 做傅里叶展开，或者用解析公式。对于圆形柱，解析公式涉及贝塞尔函数；对于方形柱或其他形状，需要用数值积分仔细处理。

解决：检查你的 ε⁻¹(G) 计算方式。如果是圆形柱，用 2.2 节的解析式；如果是任意形状，用足够密的网格做数值积分，确保 G=0 分量和 G≠0 分量都正确。

4.2 能带图出现虚假简并：截断阶数不足

现象：能带图上出现大量平直的或者不光滑的曲线，看起来像简并但实际上是数值噪声。

原因：截断阶数 N 太小，基组不完备，高频率部分的平面波没有被充分展开。频率越高，需要的平面波越多。

解决：增大 N 做收敛性测试。如果增大 N 后曲线变光滑且带隙稳定，说明之前确实不够。一般低频部分收敛快，高频部分收敛慢，如果只关心第一带隙，N=12 可能够用，但如果要看第二带隙，N 要更大。

4.3 完全禁带消失：k 点路径采样太稀

现象：明明文献说这个结构有完全禁带，你的代码算出来没有。

原因：k 点路径采样太稀，漏掉了带隙边缘的极值点。完全禁带要求对所有 k 方向都成立，如果采样点不够密，可能刚好跳过了带隙最窄的方向。

解决：增加每段的 k 点数量，从 40 增加到 80 甚至 100。另外，除了高对称点路径，最好在布里渊区边界附近加密采样，因为带隙极值通常出现在高对称点或边界上。

4.4 计算时间爆炸：矩阵维度失控

现象：N=15 时程序跑了几十分钟还没出结果。

原因：MATLAB 的 eig 对 961×961 的矩阵求解本身就需要时间，如果嵌套循环里还有大量重复计算，时间会成倍增加。另外，如果对每个 k 点都重新构造完整的矩阵并求所有本征值，计算量是 num_k × numG³ 量级。

解决：利用矩阵的稀疏性或者对称性加速。对于平面波展开法，矩阵 A 是稠密的，但可以用 eigs 求部分本征值而不是全部。如果只关心低频部分，用 eigs(A, 20, 'sm') 求最小的 20 个本征值，速度比 eig 快很多。另外，把 eps_inv 的计算提到 k 循环外面，避免重复计算。

4.5 极化模式搞混：TE 和 TM 的结果张冠李戴

现象：算出来的 TE 能带和文献的 TM 能带形状一样。

原因：TE 和 TM 的本征值矩阵构造方式不同，TE 用 ε⁻¹ 的傅里叶系数，TM 用 ε 的傅里叶系数。如果代码里写反了，结果就完全对不上。

解决：仔细检查矩阵构造部分。TE 极化（E 沿 z 方向）的方程是 Σ |k+G||k+G'| ε⁻¹(G-G') E(G') = (ω/c)² E(G)；TM 极化（H 沿 z 方向）的方程是 Σ (k+G)·(k+G') ε⁻¹(G-G') H(G') = (ω/c)² H(G)。注意点积和模长的区别。

5. 进阶技巧：用对称性加速与带隙自动判定

5.1 利用点群对称性缩减计算量

正方晶格的光子晶体具有 C4v 对称性，倒格矢空间中的平面波基组可以按对称性分类，把矩阵分块对角化。这样矩阵维度可以缩减到原来的几分之一，计算速度大幅提升。具体做法是构造对称化基组，但实现起来比较复杂，适合对性能有极致要求的情况。

更实用的加速方法是：只计算不可约布里渊区内的 k 点。正方晶格的不可约布里渊区是 Γ-X-M-Γ 围成的三角形区域，只需要在这个区域内采样，然后通过对称性展开到整个布里渊区。这样 k 点数量可以减少到原来的 1/8 左右。

5.2 自动判定完全禁带的代码实现

手动看能带图判断禁带容易出错，尤其是带隙比较窄的时候。下面是一个自动判定的代码片段：

% 自动判定完全禁带 % freqs: num_k_total x numG 的本征频率矩阵 % 思路：对所有 k 点，找共同的频率间隙 num_bands = 10; % 只关心前10条能带 freqs_trim = freqs(:, 1:num_bands); % 对每个 k 点，相邻能带之间的间隙 % 完全禁带要求：存在一个频率区间，对所有 k 点都落在某两条能带之间 % 简化方法：取所有 k 点的第 n 和第 n+1 条能带 % 检查 max(第n条) < min(第n+1条) 是否成立 gap_found = false; for n = 1:num_bands-1 band_n_max = max(freqs_trim(:, n)); band_n1_min = min(freqs_trim(:, n+1)); if band_n_max < band_n1_min gap_width = band_n1_min - band_n_max; gap_center = (band_n_max + band_n1_min) / 2; fprintf('发现完全禁带：第%d和第%d条能带之间\n', n, n+1); fprintf('带隙宽度：%.4f，中心频率：%.4f\n', gap_width, gap_center); gap_found = true; end end if ~gap_found fprintf('未发现完全禁带\n'); end

逻辑说明：完全禁带的严格定义是存在一个频率区间，使得对所有 k 点，该区间内都没有本征频率。上面的代码通过检查相邻能带的全局最大值和最小值来判断：如果第 n 条能带在所有 k 点的最大频率小于第 n+1 条能带在所有 k 点的最小频率，那么这两条能带之间就是完全禁带。参数 num_bands 控制检查的能带数量，一般前 10 条够用，因为高阶带隙在实际中很少被利用。