简介这套MATLAB代码包实现了严格耦合波分析RCWA与平面波展开法两种主流数值算法面向光学工程、电磁场仿真领域的研究人员和学生用于模拟光波在周期性介质结构中的传播、衍射以及光子晶体带隙和电磁场分布的计算。压缩包共24个文件仅17KB其中13个.m脚本构成完整实现涵盖主程序、网格设置、卷积矩阵、分层计算、波矢计算等核心模块另有8个.dat参数文件、1个.txt说明文件、1个README.md使用说明和1个.asv自动保存文件便于读者快速运行和二次修改。已有85人学习下载。代码可帮助研究者快速验证理论模型掌握傅里叶级数展开、耦合波方程组构建、PML吸收边界设置等关键实现细节并通过MATLAB绘图功能直观展示衍射效率和场分布。适用于光栅衍射效率分析、光学滤波器与分束器设计、光子晶体研究等场景是一套兼顾教学与科研的实用工具。1. 严格耦合波分析和平面波展开法的MATLAB实现从光栅衍射到光子晶体能带严格耦合波分析和平面波展开法RCWA、PWE是处理周期结构与光相互作用时绕不开的两套数值方法。RCWA把结构沿纵向分层、横向做傅里叶展开通过特征值问题算出光栅的反射、透射和各级衍射效率PWE则把光子晶体的电磁场展开到倒格子空间直接求解色散关系得到能带和带隙。很多初学者一上来就照着论文搭矩阵结果要么衍射效率不守恒要么能带里冒出一条理直气壮的假带隙最后只能归咎于“数值玄学”。这篇笔记按照我实际调试一轮的经验从数学背景、MATLAB代码到调参避坑把这两条路完整走一遍。适合正在做光栅、超表面、光子晶体能带计算想用MATLAB快速跑通首个版本的科研和工程人员。2. 开跑前必读傅里叶展开、特征值问题和两种方法的选型逻辑2.1 严格耦合波分析到底在解什么方程RCWA的思路是把一个连续介质结构“切”成若干水平层每一层内的介电常数只随横向坐标变化在纵向保持不变。只要横向介电常数满足周期性就能展成傅里叶级数光栅区内的电磁场也写成一组空间谐波。把这些展开式代进麦克斯韦方程组偏微分方程就退化成一个大代数方程组每层内部横向场分布对应一个特征值问题解出来的特征值和特征向量分别对应本征模沿纵向的传播常数和横向场形状。折射率突变表现在傅里叶系数的高阶级次不为零会让不同谐波之间发生耦合这也是“耦合波”三个字的来源。这个过程中的核心参数有两个傅里叶展开的谐波级次 N以及结构纵向分层的层数。N 直接决定矩阵尺寸是2N1×2N1N 越大越能刻画折射率突变但也会带来高阶级次的收敛问题。层数则决定了结构能否被近似成阶梯形比如一个表面浮雕光栅切 20 层和切 200 层结果可以差很多。RCWA 之所以是“严格”的是因为它保留全矢量电磁场关系不做标量近似代价是切层、建矩阵、做边界匹配每一步都要处理数值稳定性。2.2 平面波展开法到底在解什么方程PWE 的思路和 RCWA 正好反过来它不关心外部入射波只关心无限周期介质里能存在哪些本征模式。介电常数是一个周期函数先把它的倒数展开成傅里叶级数然后把电磁场展开成一系列平面波的叠加代进波动方程后同样得到一个代数本征值问题。以二维光子晶体 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}} ]其中 (\kappa) 是介电常数倒数的傅里叶系数(\mathbf{G}) 是倒格子矢量(\mathbf{k}) 是布里渊区内的波矢。对于每一个固定的 (\mathbf{k})解一次本征值问题得到的本征值就是对应模式的频率平方把 (\mathbf{k}) 沿高对称点路径扫过去就得到能带。PWE 的实现门槛比 RCWA 略低因为不需要分层也不需要边界匹配难点全在傅里叶系数和倒格矢构建上。尤其是空气孔和背景介质交替的结构圆孔的傅里叶系数有解析表达式可以直接用贝塞尔函数算不必做数值积分。PWE 的代价是假设结构无限周期、材料无损耗因此它天然不适合少数几层膜或有限尺寸结构。2.3 什么时候用RCWA、什么时候用PWE选型看你要回答的问题。RCWA 擅长回答“一束光打在这个有限厚度周期结构上反射和透射各是多少”适合光栅、超表面、抗反射结构这一类有明确入射条件的散射问题。PWE 擅长回答“这种无限周期介质里能传播哪些模式、有没有带隙”适合光子晶体、声子晶体这一类本征模式问题。对比项RCWAPWE典型输出各级衍射效率、反射率、透射率能带结构、带隙位置和宽度结构假设有限厚度、横向周期、纵向分层无限周期、纵向无边界入射条件需要定义平面波入射角、偏振不需要外部入射只扫波矢材料损耗支持复折射率一般先假设无损耗矩阵规模与谐波级次有关与平面波截断数有关我做选择时还有一个习惯如果结构在传播方向上只有一层或两层优先 RCWA因为它直给出效率曲线如果关心的是 “某个波矢下模式存不存在”直接上 PWE。两个方法的理论基础最后都落在傅里叶展开和特征值上所以拿到任何一个的编号问题另一个排查起来也不会太陌生。3. RCWA的MATLAB实现一维光栅衍射效率从建矩阵到出曲线3.1 设置入射条件与离散级次先从一个我能完全掌控的最小模型开始一维二元光栅TE 偏振无吸收介质。这种结构介电常数剖面是周期方波傅里叶系数有闭合表达式适合验证代码逻辑。参数全部用纳米为单位入射角用弧度。% 严格耦合波分析RCWA最小示例一维二元光栅TE偏振 lambda 633; % 入射波长单位 nm period 600; % 光栅周期 duty 0.5; % 占空比脊宽度 / 周期 depth 400; % 光栅深度 theta 0; % 入射角0 表示正入射 n_i 1.0; % 入射区折射率空气 n_g 1.5; % 光栅脊折射率 n_s 1.5; % 基底折射率 N 15; % 空间谐波截断级次 m (-N:N).; % 级次索引总数 M 2*N1 M length(m);这里 N 的初始值先不要设太大15 级对应 31 个谐波矩阵是 31×31一次特征值求解几乎是瞬时的。先跑通整个流程确认能量守恒没问题再把 N 往上加。级次索引 m 是傅里叶展开的阶数它在后面既决定横向波矢匹配也决定介电常数傅里叶系数。3.2 组装配对矩阵与特征值求解RCWA 最核心的一步是把介电常数函数变成 Toeplitz 矩阵。Toeplitz 矩阵的第 i 行第 j 列只取决于级次差刚好对应傅里叶级数卷积的性质。二元光栅的傅里叶系数可以用解析公式不用做数值积分。k0 2*pi/lambda; % 横向波矢匹配条件入射波矢与各衍射级次 kx_m n_i*k0*sin(theta) - m*2*pi/period; Kx diag(kx_m/k0); % 介电常数的傅里叶系数 eps_avg duty*n_g^2 (1-duty)*n_i^2; delta_eps n_g^2 - n_i^2; eps_coeff zeros(M,1); for idx 1:M mm m(idx); if mm 0 eps_coeff(idx) eps_avg; else eps_coeff(idx) delta_eps * sin(pi*duty*mm)/(pi*mm); end end E toeplitz(eps_coeff); % 耦合矩阵 E(m,n) eps_{m-n} % TE 偏振特征值问题 A Kx^2 - E; [W, D] eig(A); q2 diag(D); q sqrt(q2); Q diag(q); V W*Q;toeplitz函数接受一个列向量作为第一列和第一行的关联序列正好把eps_coeff展开成完整的差卷积矩阵。A Kx^2 - E是 TE 偏振下的标准特征矩阵特征向量 W 的每一列代表一个本征模在各级次上的投影Q 的对角元素是本征模的纵向传播常数。注意这里的横向波矢匹配条件已经包含入射角 θ如果周期和波长接近kx_m/k0的绝对值可能大于 1对应倏逝波特征值出现复数不是程序错误而是物理现象。3.3 组装边界条件并算衍射效率有了光栅层内的本征模还需要把入射区、光栅层、基底区三个区域的场在前后两个界面处拼接。标准做法是把反射幅度、光栅层前向波、背向波、透射幅度都作为未知量列一个块状线性方程组一次性解出来。% 入射区和基底区的纵向波矢归一化 kz_i sqrt(n_i^2 - (kx_m/k0).^2); kz_s sqrt(n_s^2 - (kx_m/k0).^2); % 让辐射波满足边界条件反射波向 -z透射波向 z kz_i real(kz_i) 1i*abs(imag(kz_i)); kz_s real(kz_s) 1i*abs(imag(kz_s)); YI diag(kz_i); YS diag(kz_s); % 纵向传播相位 Pp exp( 1i*k0*Q*depth ); Pm exp( -1i*k0*Q*depth ); % 未知量排列R(反射)F(前向波)B(背向波)T(透射) Z zeros(M); I eye(M); e0 zeros(M,1); e0(N1) 1; % 入射波只存在于 0 级 Amat [ I, -W, -W, Z; YI, -V, V, Z; Z, W*Pp, W*Pm, -I; Z, V*Pp, -V*Pm, -YS ]; b [ -e0; -YI*e0; zeros(M,1); zeros(M,1) ]; x Amat \ b; Rf x(1:M); Tf x(3*M1:4*M); % 衍射效率透射沿 z反射沿 -z取绝对值保证能量为正 DE_r -abs(Rf).^2 .* real(kz_i) ./ real(kz_i(1)); DE_t abs(Tf).^2 .* real(kz_s) ./ real(kz_i(1)); fprintf(反射效率总和: %.6f\n, sum(DE_r)); fprintf(透射效率总和: %.6f\n, sum(DE_t));这段代码是典型的矩阵组装法适合深度不超过几个波长的介质光栅。边界条件的物理意义是界面处电场切向分量连续、磁场切向分量连续入射场的 0 级幅度预先置为 1其余未知量由方程组决定。上面的DE_r带了负号因为反射波的能量通量沿 -z 方向而参考基准是入射波沿 z 的通量如果符号约定和 V 矩阵不一致效率可能出现负值这是常见现象调整A或V的符号即可。更稳妥的做法是在无耗结构里检查sum(DE_r)sum(DE_t)是否等于 1这是每一版代码都要跑的硬指标。4. PWE的MATLAB实现二维光子晶体的能带结构与带隙扫描4.1 倒格子基矢与平面波截断PWE 的第一步是生成倒格子矢量集合这一步决定了后续矩阵的维度和频率精确度。以正方晶格为例实空间周期常数设为 a倒格子基矢就是 (2\pi/a) 的正交向量。平面波截断级数 N 表示每个方向上取 -N 到 N 的倒格矢平面波总数是 ((2N1)^2)。% 平面波展开法PWE最小示例正方晶格圆柱空气孔TM 模式 a 1.0; % 晶格常数归一化 r 0.3*a; % 空气孔半径 f pi*r^2/a^2; % 空气孔填充率 eps_b 12.0; % 背景介质介电常数 eps_a 1.0; % 空气孔介电常数 N 4; % 平面波截断级数 n -N:N; [nx, ny] meshgrid(n, n); Gx nx*2*pi/a; Gy ny*2*pi/a; G [Gx(:), Gy(:)]; NG length(G); % 平面波总数 (2*N1)^2N 4 时 NG 81矩阵尺寸是 81×81特征值求解很快。N 增大到 8 时 NG 289内存和耗时都会明显上升。对于大多数能带结构验证N 4 到 6 足够看禁带趋势N 8 以上适合确认收敛后的精确带边位置。这里我把所有长度都用 a 的倍数归一化后面输出的频率自动是无量纲的 (a/\lambda)。4.2 组装本征矩阵并求解色散关系PWE 的关键是介电常数倒数傅里叶系数 (\kappa)。对圆孔结构可以预先计算所有倒格矢差对应的结构因子 S再组合成 (\kappa) 矩阵。随后沿高对称点路径 (\Gamma-X-M-\Gamma) 逐点扫描 k组装 TM 模式的本征矩阵并求特征值。% 预计算任意两个倒格矢之差 Gdiff reshape(G, [NG,1,2]) - reshape(G, [1,NG,2]); Gd_norm sqrt(sum(Gdiff.^2, 3)); % 圆孔结构因子 S(G) S zeros(NG, NG); mask Gd_norm 1e-12; S(~mask) f; S(mask) 2*f * besselj(1, Gd_norm(mask)*r) ./ (Gd_norm(mask)*r); % 介电常数倒数的傅里叶系数 eta 1/eps_b (1/eps_a - 1/eps_b) * S; % 布里渊区高对称点 Gamma [0,0]*pi/a; X [1,0]*pi/a; M [1,1]*pi/a; npoints 40; seg1 [linspace(Gamma(1),X(1),npoints); linspace(Gamma(2),X(2),npoints)].; seg2 [linspace(X(1),M(1),npoints); linspace(X(2),M(2),npoints)].; seg3 [linspace(M(1),Gamma(1),npoints); linspace(M(2),Gamma(2),npoints)].; kpath [seg1; seg2; seg3]; % 逐 k 点求解特征值 nband 8; % 需要的最低频带数 eigs_all zeros(size(kpath,1), nband); for ik 1:size(kpath,1) k kpath(ik,:); Mmat zeros(NG, NG); for iG1 1:NG kG1 k G(iG1,:); for iG2 1:NG kG2 k G(iG2,:); Mmat(iG1,iG2) eta(iG1,iG2) * dot(kG1, kG2); end end Mmat 0.5*(MmatMmat.); % 强制对称避免数值噪声影响 [~, D] eig(Mmat); ev real(diag(D)); ev(ev 0) 0; % 极小负特征值按零处理 ev_sorted sort(ev); eigs_all(ik,:) ev_sorted(1:nband); end freq sqrt(eigs_all) / (2*pi) * a;这里的双重循环在 NG 较小时完全够用。N 4 时每个 k 点要循环 81×81 次路径上一共 120 个 k 点几秒钟就能跑完。当 N 提高到 6 或 8双重循环会明显变慢我一般用向量的差预先算好 eta 矩阵再用稀疏矩阵或分块计算优化而不是直接套双重循环否则内存会先撑不住。注意ev(ev 0) 0这只是清理数值上可能出现的微小负特征值不是物理修正如果出现大量负特征值说明介电常数系数或倒格矢构造有误。4.3 能带绘制与带隙确认能带图是最直观的验证结果。横轴用线性索引而不是真实 k 距离因为不同段的长度不一样但画图时想要每个高对称点等距分布这是常见做法。xk 1:size(kpath,1); tick_pos [1, npoints, 2*npoints, 3*npoints]; tick_label {\Gamma, X, M, \Gamma}; plot(xk, freq, .-, MarkerSize, 3); grid on; xlim([1, size(kpath,1)]); set(gca, XTick, tick_pos, XTickLabel, tick_label); ylabel(归一化频率 a/\lambda);画完后先看最低频率带是否从 Gamma 点附近的零开始。TM 模式的最低带在长波极限下应该趋近于零这是能带理论的基本预期。如果最低点明显偏高多半是结构因子 S 或介电常数系数写错。接着看相邻带之间是否出现完全分离的空白区间那就是带隙。带隙宽度要到高频带和低频带都收敛后才能下结论否则很容易把数值误差看成带隙。5. RCWA与PWE的调参避坑与常见问题排查5.1 RCWA衍射效率不收敛级次截断的玄学现象N 从 5 加到 310 级衍射效率还在 0.7 到 0.9 之间来回漂看不出稳定趋势。原因二元光栅的折射率突变在傅里叶空间里衰减很慢高阶级次贡献不可忽略当光栅深度和波长可比拟时截断误差会直接反映在效率上。解决不要凭经验定 N先画一条收敛曲线。把 N 从 2 扫到 20看目标效率变化什么时候小于 1e-3。这里的判断标准和你要找共振还是算背景效率有关算超表面共振往往需要 N 更大。5.2 PWE平面波截断不够高频带变成锯齿现象能带低频区很光滑高频区出现折返、锯齿甚至出现上下带交叉。原因平面波数量不足时高频模式对应的倒格矢信息被截断矩阵只能在有限基里近似描述真实模式。解决把 N 从 4 提到 6 或 8比较最低几条带的变化量。如果低频带位置只移动不到 1%就可以认为低频结论可靠。高频带本来就不应该直接用于带隙判断带隙上下沿通常涉及较低频带所以优先关心带边收敛。5.3 假带隙结构因子或k路径错位现象算出来的带隙和文献或解析预期明显对不上比如理论上闭合的带隙在这里却张开了。原因一个常见问题是空气孔半径和填充率换算错误例如把半径直接当直径用另一个问题是 k 路径坐标没有用弧度单位导致倒格矢缩放全部偏掉。解决先用一个已知例子校准代码。我习惯把 r 0.3a 的正方晶格空气孔 TM 能带跑一遍和标准数值结果核对带边位置。再有就是打印 S(0) 是不是等于填充率 f圆孔在 G0 处的结构因子一定有这个性质如果这里错了后面全部白算。5.4 TM偏振的RCWA矩阵奇异介电常数矩阵取逆的经典坑现象TE 偏振能跑通切到 TM 偏振后矩阵出现 NaN或者特征值在折射率对比度增大时突然爆炸。原因TM 偏振的本构关系涉及电位移矢量的分量数值格式里需要用介电常数倒数的 Toeplitz 矩阵而不是先把 E 矩阵求逆再用。很多人理解成“取逆”结果得到完全不同的矩阵结构收敛性非常差。解决构建介电常数倒数 (1/\varepsilon(x)) 的傅里叶系数用toeplitz组装成 (E_{\mathrm{inv}}) 矩阵再参与 TM 的特征矩阵运算。如果整个代码里只出现了一个inv(E)基本可以肯定是踩了这个坑。TM 的收敛通常比 TE 慢要多留几个级次的余量。5.5 能量守恒检查反射加透射不等于1现象无耗介质光栅算出 sum(DE_r) sum(DE_t) 1.15 或 0.82怎么调参数都救不回来。原因特征值平方根的符号选取不统一尤其是倏逝波的虚部符号错误会让边界条件方程组把能量算错方向也可能是光栅深度太大时Pa 和 Pm 指数项相差悬殊线性求解器出现数值病态。解决先检查传播级次的数量确认效率求和时只计入real(kz_i)0和real(kz_s)0的通道。再检查 q 的符号统一保证虚部为正或实部为正让前向波和背向波的相位方向一致。深度超过几个波长时不要用直接矩阵求解改用 S 矩阵递推这样能避免指数爆炸。每一步改动后立刻跑能量守恒通常一两轮就能把符号和方向的问题揪出来。6. 让计算更稳更快的进阶做法参数扫描、并行化与输出标准化6.1 对衍射级次做自适应截断每次手动调 N 很浪费。我一般把 RCWA 包成一个函数然后在外面做自适应循环从 N 2 开始算每次把 N 翻倍或加 5看目标效率变化是否小于设定阈值满足就停止。这种方法比直接设一个很大的 N 省时间也能避免不同参数下精度不一致的问题。N 2; DE_old 0; while N 30 [~, DE_t] rcwa1d(lambda, period, duty, depth, N, theta); if abs(DE_t(1) - DE_old) 1e-4 break; end DE_old DE_t(1); N N 2; end收敛阈值不能设得太夸张设备加工误差带来的效率变化一般也在百分之一以上追到 1e-4 已经足够判断趋势。6.2 能带扫描的高效循环与并行化PWE 的 k 点之间可以并行因为每个 k 点的特征值求解互相独立。先确保矩阵构建是纯函数再把外层循环换成parfor注意每个线程只写eigs_all(ik,:)这一行。如果本机没有并行工具箱可以把 k 路径拆成三段时间分别求解再拼起来内存压力反而更小。parfor ik 1:size(kpath,1) eigs_all(ik,:) pwe2d_tm_one_k(kpath(ik,:), G, eta, NG, nband); endpwe2d_tm_one_k只负责某个 k 点的矩阵组装和特征值求解保证循环体内部没有共享变量。参数扫描线下发任务时优先把 k 点和结构参数放在最外层避免频繁重建倒格矢。6.3 从MATLAB脚本变成可复用函数库脚本跑通后第一件事是把参数入口收敛成函数命名要能看出物理量。function [DE_r, DE_t, angles] rcwa1d(varargin) % varargin: lambda, period, duty, depth, N, theta, n_i, n_g, n_s end function [freq, kdist] pwe2d_tm(a, r, eps_b, eps_a, N, kpath) % 输出频带与归一化 k 坐标 end我后期还会强制返回一个结果结构体包含所有输入参数、收敛级次、能量守恒偏差这样每次仿真记录都能回溯。自己有过一次教训同一套光栅参数换了入射角几天后想复现却想不起当时用的级次只能重新扫一遍浪费不少时间。从那以后我每次都把收敛历史一起存下来。希望帮到你。本文还有配套的精品资源点击获取
简介这套MATLAB代码包实现了严格耦合波分析RCWA与平面波展开法两种主流数值算法面向光学工程、电磁场仿真领域的研究人员和学生用于模拟光波在周期性介质结构中的传播、衍射以及光子晶体带隙和电磁场分布的计算。压缩包共24个文件仅17KB其中13个.m脚本构成完整实现涵盖主程序、网格设置、卷积矩阵、分层计算、波矢计算等核心模块另有8个.dat参数文件、1个.txt说明文件、1个README.md使用说明和1个.asv自动保存文件便于读者快速运行和二次修改。已有85人学习下载。代码可帮助研究者快速验证理论模型掌握傅里叶级数展开、耦合波方程组构建、PML吸收边界设置等关键实现细节并通过MATLAB绘图功能直观展示衍射效率和场分布。适用于光栅衍射效率分析、光学滤波器与分束器设计、光子晶体研究等场景是一套兼顾教学与科研的实用工具。1. 严格耦合波分析和平面波展开法的MATLAB实现从光栅衍射到光子晶体能带严格耦合波分析和平面波展开法RCWA、PWE是处理周期结构与光相互作用时绕不开的两套数值方法。RCWA把结构沿纵向分层、横向做傅里叶展开通过特征值问题算出光栅的反射、透射和各级衍射效率PWE则把光子晶体的电磁场展开到倒格子空间直接求解色散关系得到能带和带隙。很多初学者一上来就照着论文搭矩阵结果要么衍射效率不守恒要么能带里冒出一条理直气壮的假带隙最后只能归咎于“数值玄学”。这篇笔记按照我实际调试一轮的经验从数学背景、MATLAB代码到调参避坑把这两条路完整走一遍。适合正在做光栅、超表面、光子晶体能带计算想用MATLAB快速跑通首个版本的科研和工程人员。2. 开跑前必读傅里叶展开、特征值问题和两种方法的选型逻辑2.1 严格耦合波分析到底在解什么方程RCWA的思路是把一个连续介质结构“切”成若干水平层每一层内的介电常数只随横向坐标变化在纵向保持不变。只要横向介电常数满足周期性就能展成傅里叶级数光栅区内的电磁场也写成一组空间谐波。把这些展开式代进麦克斯韦方程组偏微分方程就退化成一个大代数方程组每层内部横向场分布对应一个特征值问题解出来的特征值和特征向量分别对应本征模沿纵向的传播常数和横向场形状。折射率突变表现在傅里叶系数的高阶级次不为零会让不同谐波之间发生耦合这也是“耦合波”三个字的来源。这个过程中的核心参数有两个傅里叶展开的谐波级次 N以及结构纵向分层的层数。N 直接决定矩阵尺寸是2N1×2N1N 越大越能刻画折射率突变但也会带来高阶级次的收敛问题。层数则决定了结构能否被近似成阶梯形比如一个表面浮雕光栅切 20 层和切 200 层结果可以差很多。RCWA 之所以是“严格”的是因为它保留全矢量电磁场关系不做标量近似代价是切层、建矩阵、做边界匹配每一步都要处理数值稳定性。2.2 平面波展开法到底在解什么方程PWE 的思路和 RCWA 正好反过来它不关心外部入射波只关心无限周期介质里能存在哪些本征模式。介电常数是一个周期函数先把它的倒数展开成傅里叶级数然后把电磁场展开成一系列平面波的叠加代进波动方程后同样得到一个代数本征值问题。以二维光子晶体 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}} ]其中 (\kappa) 是介电常数倒数的傅里叶系数(\mathbf{G}) 是倒格子矢量(\mathbf{k}) 是布里渊区内的波矢。对于每一个固定的 (\mathbf{k})解一次本征值问题得到的本征值就是对应模式的频率平方把 (\mathbf{k}) 沿高对称点路径扫过去就得到能带。PWE 的实现门槛比 RCWA 略低因为不需要分层也不需要边界匹配难点全在傅里叶系数和倒格矢构建上。尤其是空气孔和背景介质交替的结构圆孔的傅里叶系数有解析表达式可以直接用贝塞尔函数算不必做数值积分。PWE 的代价是假设结构无限周期、材料无损耗因此它天然不适合少数几层膜或有限尺寸结构。2.3 什么时候用RCWA、什么时候用PWE选型看你要回答的问题。RCWA 擅长回答“一束光打在这个有限厚度周期结构上反射和透射各是多少”适合光栅、超表面、抗反射结构这一类有明确入射条件的散射问题。PWE 擅长回答“这种无限周期介质里能传播哪些模式、有没有带隙”适合光子晶体、声子晶体这一类本征模式问题。对比项RCWAPWE典型输出各级衍射效率、反射率、透射率能带结构、带隙位置和宽度结构假设有限厚度、横向周期、纵向分层无限周期、纵向无边界入射条件需要定义平面波入射角、偏振不需要外部入射只扫波矢材料损耗支持复折射率一般先假设无损耗矩阵规模与谐波级次有关与平面波截断数有关我做选择时还有一个习惯如果结构在传播方向上只有一层或两层优先 RCWA因为它直给出效率曲线如果关心的是 “某个波矢下模式存不存在”直接上 PWE。两个方法的理论基础最后都落在傅里叶展开和特征值上所以拿到任何一个的编号问题另一个排查起来也不会太陌生。3. RCWA的MATLAB实现一维光栅衍射效率从建矩阵到出曲线3.1 设置入射条件与离散级次先从一个我能完全掌控的最小模型开始一维二元光栅TE 偏振无吸收介质。这种结构介电常数剖面是周期方波傅里叶系数有闭合表达式适合验证代码逻辑。参数全部用纳米为单位入射角用弧度。% 严格耦合波分析RCWA最小示例一维二元光栅TE偏振 lambda 633; % 入射波长单位 nm period 600; % 光栅周期 duty 0.5; % 占空比脊宽度 / 周期 depth 400; % 光栅深度 theta 0; % 入射角0 表示正入射 n_i 1.0; % 入射区折射率空气 n_g 1.5; % 光栅脊折射率 n_s 1.5; % 基底折射率 N 15; % 空间谐波截断级次 m (-N:N).; % 级次索引总数 M 2*N1 M length(m);这里 N 的初始值先不要设太大15 级对应 31 个谐波矩阵是 31×31一次特征值求解几乎是瞬时的。先跑通整个流程确认能量守恒没问题再把 N 往上加。级次索引 m 是傅里叶展开的阶数它在后面既决定横向波矢匹配也决定介电常数傅里叶系数。3.2 组装配对矩阵与特征值求解RCWA 最核心的一步是把介电常数函数变成 Toeplitz 矩阵。Toeplitz 矩阵的第 i 行第 j 列只取决于级次差刚好对应傅里叶级数卷积的性质。二元光栅的傅里叶系数可以用解析公式不用做数值积分。k0 2*pi/lambda; % 横向波矢匹配条件入射波矢与各衍射级次 kx_m n_i*k0*sin(theta) - m*2*pi/period; Kx diag(kx_m/k0); % 介电常数的傅里叶系数 eps_avg duty*n_g^2 (1-duty)*n_i^2; delta_eps n_g^2 - n_i^2; eps_coeff zeros(M,1); for idx 1:M mm m(idx); if mm 0 eps_coeff(idx) eps_avg; else eps_coeff(idx) delta_eps * sin(pi*duty*mm)/(pi*mm); end end E toeplitz(eps_coeff); % 耦合矩阵 E(m,n) eps_{m-n} % TE 偏振特征值问题 A Kx^2 - E; [W, D] eig(A); q2 diag(D); q sqrt(q2); Q diag(q); V W*Q;toeplitz函数接受一个列向量作为第一列和第一行的关联序列正好把eps_coeff展开成完整的差卷积矩阵。A Kx^2 - E是 TE 偏振下的标准特征矩阵特征向量 W 的每一列代表一个本征模在各级次上的投影Q 的对角元素是本征模的纵向传播常数。注意这里的横向波矢匹配条件已经包含入射角 θ如果周期和波长接近kx_m/k0的绝对值可能大于 1对应倏逝波特征值出现复数不是程序错误而是物理现象。3.3 组装边界条件并算衍射效率有了光栅层内的本征模还需要把入射区、光栅层、基底区三个区域的场在前后两个界面处拼接。标准做法是把反射幅度、光栅层前向波、背向波、透射幅度都作为未知量列一个块状线性方程组一次性解出来。% 入射区和基底区的纵向波矢归一化 kz_i sqrt(n_i^2 - (kx_m/k0).^2); kz_s sqrt(n_s^2 - (kx_m/k0).^2); % 让辐射波满足边界条件反射波向 -z透射波向 z kz_i real(kz_i) 1i*abs(imag(kz_i)); kz_s real(kz_s) 1i*abs(imag(kz_s)); YI diag(kz_i); YS diag(kz_s); % 纵向传播相位 Pp exp( 1i*k0*Q*depth ); Pm exp( -1i*k0*Q*depth ); % 未知量排列R(反射)F(前向波)B(背向波)T(透射) Z zeros(M); I eye(M); e0 zeros(M,1); e0(N1) 1; % 入射波只存在于 0 级 Amat [ I, -W, -W, Z; YI, -V, V, Z; Z, W*Pp, W*Pm, -I; Z, V*Pp, -V*Pm, -YS ]; b [ -e0; -YI*e0; zeros(M,1); zeros(M,1) ]; x Amat \ b; Rf x(1:M); Tf x(3*M1:4*M); % 衍射效率透射沿 z反射沿 -z取绝对值保证能量为正 DE_r -abs(Rf).^2 .* real(kz_i) ./ real(kz_i(1)); DE_t abs(Tf).^2 .* real(kz_s) ./ real(kz_i(1)); fprintf(反射效率总和: %.6f\n, sum(DE_r)); fprintf(透射效率总和: %.6f\n, sum(DE_t));这段代码是典型的矩阵组装法适合深度不超过几个波长的介质光栅。边界条件的物理意义是界面处电场切向分量连续、磁场切向分量连续入射场的 0 级幅度预先置为 1其余未知量由方程组决定。上面的DE_r带了负号因为反射波的能量通量沿 -z 方向而参考基准是入射波沿 z 的通量如果符号约定和 V 矩阵不一致效率可能出现负值这是常见现象调整A或V的符号即可。更稳妥的做法是在无耗结构里检查sum(DE_r)sum(DE_t)是否等于 1这是每一版代码都要跑的硬指标。4. PWE的MATLAB实现二维光子晶体的能带结构与带隙扫描4.1 倒格子基矢与平面波截断PWE 的第一步是生成倒格子矢量集合这一步决定了后续矩阵的维度和频率精确度。以正方晶格为例实空间周期常数设为 a倒格子基矢就是 (2\pi/a) 的正交向量。平面波截断级数 N 表示每个方向上取 -N 到 N 的倒格矢平面波总数是 ((2N1)^2)。% 平面波展开法PWE最小示例正方晶格圆柱空气孔TM 模式 a 1.0; % 晶格常数归一化 r 0.3*a; % 空气孔半径 f pi*r^2/a^2; % 空气孔填充率 eps_b 12.0; % 背景介质介电常数 eps_a 1.0; % 空气孔介电常数 N 4; % 平面波截断级数 n -N:N; [nx, ny] meshgrid(n, n); Gx nx*2*pi/a; Gy ny*2*pi/a; G [Gx(:), Gy(:)]; NG length(G); % 平面波总数 (2*N1)^2N 4 时 NG 81矩阵尺寸是 81×81特征值求解很快。N 增大到 8 时 NG 289内存和耗时都会明显上升。对于大多数能带结构验证N 4 到 6 足够看禁带趋势N 8 以上适合确认收敛后的精确带边位置。这里我把所有长度都用 a 的倍数归一化后面输出的频率自动是无量纲的 (a/\lambda)。4.2 组装本征矩阵并求解色散关系PWE 的关键是介电常数倒数傅里叶系数 (\kappa)。对圆孔结构可以预先计算所有倒格矢差对应的结构因子 S再组合成 (\kappa) 矩阵。随后沿高对称点路径 (\Gamma-X-M-\Gamma) 逐点扫描 k组装 TM 模式的本征矩阵并求特征值。% 预计算任意两个倒格矢之差 Gdiff reshape(G, [NG,1,2]) - reshape(G, [1,NG,2]); Gd_norm sqrt(sum(Gdiff.^2, 3)); % 圆孔结构因子 S(G) S zeros(NG, NG); mask Gd_norm 1e-12; S(~mask) f; S(mask) 2*f * besselj(1, Gd_norm(mask)*r) ./ (Gd_norm(mask)*r); % 介电常数倒数的傅里叶系数 eta 1/eps_b (1/eps_a - 1/eps_b) * S; % 布里渊区高对称点 Gamma [0,0]*pi/a; X [1,0]*pi/a; M [1,1]*pi/a; npoints 40; seg1 [linspace(Gamma(1),X(1),npoints); linspace(Gamma(2),X(2),npoints)].; seg2 [linspace(X(1),M(1),npoints); linspace(X(2),M(2),npoints)].; seg3 [linspace(M(1),Gamma(1),npoints); linspace(M(2),Gamma(2),npoints)].; kpath [seg1; seg2; seg3]; % 逐 k 点求解特征值 nband 8; % 需要的最低频带数 eigs_all zeros(size(kpath,1), nband); for ik 1:size(kpath,1) k kpath(ik,:); Mmat zeros(NG, NG); for iG1 1:NG kG1 k G(iG1,:); for iG2 1:NG kG2 k G(iG2,:); Mmat(iG1,iG2) eta(iG1,iG2) * dot(kG1, kG2); end end Mmat 0.5*(MmatMmat.); % 强制对称避免数值噪声影响 [~, D] eig(Mmat); ev real(diag(D)); ev(ev 0) 0; % 极小负特征值按零处理 ev_sorted sort(ev); eigs_all(ik,:) ev_sorted(1:nband); end freq sqrt(eigs_all) / (2*pi) * a;这里的双重循环在 NG 较小时完全够用。N 4 时每个 k 点要循环 81×81 次路径上一共 120 个 k 点几秒钟就能跑完。当 N 提高到 6 或 8双重循环会明显变慢我一般用向量的差预先算好 eta 矩阵再用稀疏矩阵或分块计算优化而不是直接套双重循环否则内存会先撑不住。注意ev(ev 0) 0这只是清理数值上可能出现的微小负特征值不是物理修正如果出现大量负特征值说明介电常数系数或倒格矢构造有误。4.3 能带绘制与带隙确认能带图是最直观的验证结果。横轴用线性索引而不是真实 k 距离因为不同段的长度不一样但画图时想要每个高对称点等距分布这是常见做法。xk 1:size(kpath,1); tick_pos [1, npoints, 2*npoints, 3*npoints]; tick_label {\Gamma, X, M, \Gamma}; plot(xk, freq, .-, MarkerSize, 3); grid on; xlim([1, size(kpath,1)]); set(gca, XTick, tick_pos, XTickLabel, tick_label); ylabel(归一化频率 a/\lambda);画完后先看最低频率带是否从 Gamma 点附近的零开始。TM 模式的最低带在长波极限下应该趋近于零这是能带理论的基本预期。如果最低点明显偏高多半是结构因子 S 或介电常数系数写错。接着看相邻带之间是否出现完全分离的空白区间那就是带隙。带隙宽度要到高频带和低频带都收敛后才能下结论否则很容易把数值误差看成带隙。5. RCWA与PWE的调参避坑与常见问题排查5.1 RCWA衍射效率不收敛级次截断的玄学现象N 从 5 加到 310 级衍射效率还在 0.7 到 0.9 之间来回漂看不出稳定趋势。原因二元光栅的折射率突变在傅里叶空间里衰减很慢高阶级次贡献不可忽略当光栅深度和波长可比拟时截断误差会直接反映在效率上。解决不要凭经验定 N先画一条收敛曲线。把 N 从 2 扫到 20看目标效率变化什么时候小于 1e-3。这里的判断标准和你要找共振还是算背景效率有关算超表面共振往往需要 N 更大。5.2 PWE平面波截断不够高频带变成锯齿现象能带低频区很光滑高频区出现折返、锯齿甚至出现上下带交叉。原因平面波数量不足时高频模式对应的倒格矢信息被截断矩阵只能在有限基里近似描述真实模式。解决把 N 从 4 提到 6 或 8比较最低几条带的变化量。如果低频带位置只移动不到 1%就可以认为低频结论可靠。高频带本来就不应该直接用于带隙判断带隙上下沿通常涉及较低频带所以优先关心带边收敛。5.3 假带隙结构因子或k路径错位现象算出来的带隙和文献或解析预期明显对不上比如理论上闭合的带隙在这里却张开了。原因一个常见问题是空气孔半径和填充率换算错误例如把半径直接当直径用另一个问题是 k 路径坐标没有用弧度单位导致倒格矢缩放全部偏掉。解决先用一个已知例子校准代码。我习惯把 r 0.3a 的正方晶格空气孔 TM 能带跑一遍和标准数值结果核对带边位置。再有就是打印 S(0) 是不是等于填充率 f圆孔在 G0 处的结构因子一定有这个性质如果这里错了后面全部白算。5.4 TM偏振的RCWA矩阵奇异介电常数矩阵取逆的经典坑现象TE 偏振能跑通切到 TM 偏振后矩阵出现 NaN或者特征值在折射率对比度增大时突然爆炸。原因TM 偏振的本构关系涉及电位移矢量的分量数值格式里需要用介电常数倒数的 Toeplitz 矩阵而不是先把 E 矩阵求逆再用。很多人理解成“取逆”结果得到完全不同的矩阵结构收敛性非常差。解决构建介电常数倒数 (1/\varepsilon(x)) 的傅里叶系数用toeplitz组装成 (E_{\mathrm{inv}}) 矩阵再参与 TM 的特征矩阵运算。如果整个代码里只出现了一个inv(E)基本可以肯定是踩了这个坑。TM 的收敛通常比 TE 慢要多留几个级次的余量。5.5 能量守恒检查反射加透射不等于1现象无耗介质光栅算出 sum(DE_r) sum(DE_t) 1.15 或 0.82怎么调参数都救不回来。原因特征值平方根的符号选取不统一尤其是倏逝波的虚部符号错误会让边界条件方程组把能量算错方向也可能是光栅深度太大时Pa 和 Pm 指数项相差悬殊线性求解器出现数值病态。解决先检查传播级次的数量确认效率求和时只计入real(kz_i)0和real(kz_s)0的通道。再检查 q 的符号统一保证虚部为正或实部为正让前向波和背向波的相位方向一致。深度超过几个波长时不要用直接矩阵求解改用 S 矩阵递推这样能避免指数爆炸。每一步改动后立刻跑能量守恒通常一两轮就能把符号和方向的问题揪出来。6. 让计算更稳更快的进阶做法参数扫描、并行化与输出标准化6.1 对衍射级次做自适应截断每次手动调 N 很浪费。我一般把 RCWA 包成一个函数然后在外面做自适应循环从 N 2 开始算每次把 N 翻倍或加 5看目标效率变化是否小于设定阈值满足就停止。这种方法比直接设一个很大的 N 省时间也能避免不同参数下精度不一致的问题。N 2; DE_old 0; while N 30 [~, DE_t] rcwa1d(lambda, period, duty, depth, N, theta); if abs(DE_t(1) - DE_old) 1e-4 break; end DE_old DE_t(1); N N 2; end收敛阈值不能设得太夸张设备加工误差带来的效率变化一般也在百分之一以上追到 1e-4 已经足够判断趋势。6.2 能带扫描的高效循环与并行化PWE 的 k 点之间可以并行因为每个 k 点的特征值求解互相独立。先确保矩阵构建是纯函数再把外层循环换成parfor注意每个线程只写eigs_all(ik,:)这一行。如果本机没有并行工具箱可以把 k 路径拆成三段时间分别求解再拼起来内存压力反而更小。parfor ik 1:size(kpath,1) eigs_all(ik,:) pwe2d_tm_one_k(kpath(ik,:), G, eta, NG, nband); endpwe2d_tm_one_k只负责某个 k 点的矩阵组装和特征值求解保证循环体内部没有共享变量。参数扫描线下发任务时优先把 k 点和结构参数放在最外层避免频繁重建倒格矢。6.3 从MATLAB脚本变成可复用函数库脚本跑通后第一件事是把参数入口收敛成函数命名要能看出物理量。function [DE_r, DE_t, angles] rcwa1d(varargin) % varargin: lambda, period, duty, depth, N, theta, n_i, n_g, n_s end function [freq, kdist] pwe2d_tm(a, r, eps_b, eps_a, N, kpath) % 输出频带与归一化 k 坐标 end我后期还会强制返回一个结果结构体包含所有输入参数、收敛级次、能量守恒偏差这样每次仿真记录都能回溯。自己有过一次教训同一套光栅参数换了入射角几天后想复现却想不起当时用的级次只能重新扫一遍浪费不少时间。从那以后我每次都把收敛历史一起存下来。希望帮到你。本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
