简介RCWA严格耦合波分析是光栅衍射特性仿真中的核心数值方法这份资源以MATLAB实现了一维光栅TE模式下的RCWA计算程序适合从事微纳光学、衍射光学元件设计的研究人员或相关专业学生用于快速求解光栅的衍射效率并验证理论模型。压缩包共10个文件其中6个.m脚本可实现主程序RCWA.m、特征值计算eig_M.m、波矢求解kvect.m/kvect1.m与介电常数转换eps_conv.m等核心功能另有4个.asv为MATLAB自动备份版本便于追溯修改过程整包仅5KB体量轻巧、便于直接读取和二次开发。目前已有1151人学习/下载。通过该程序可计算不同入射角与波长下各衍射级的效率并能结合MATLAB的矩阵运算与绘图功能分析TE偏振光在一维光栅中的传播行为适合作为光栅严格耦合波分析的入门模板或教学示例也可在此基础上扩展至二维光栅及TM模式研究。1. 光栅衍射效率计算为什么绕不开严格耦合波分析做衍射光栅的工程师手里大概率囤过不止一个 RCWA 压缩包。类似“RCWA程序.rar”这样的文件很多是从老旧的下载站转存来的打开之后是命名混乱的 m 文件和几张 demo 图。RCWA 的全称是 Rigorous Coupled-Wave Analysis中文语境里通常叫严格耦合波分析它不是一个软件而是一类算法给定周期结构与入射平面波直接求出 0 级、±1 级、±2 级等各衍射级次的效率与相位。FDTD 那种在网格上逐步推进的方式当然也能算但亚波长光栅按纳米尺度加密网格后内存和时间开销会迅速失控RCWA 利用周期性做谱展开把麦克斯韦方程压缩成矩阵特征值问题计算速度往往快一到两个数量级。这篇文章就围绕标题里的四个关键词——MATLAB、光栅、严格耦合、衍射效率把 RCWA 的原理、参数和 MATLAB 落地路径写清楚。2. 严格耦合波分析建模从周期结构到特征值问题2.1 周期性光栅与介电常数谐波展开光栅本质上是一个空间周期性变化的介电常数结构。设光栅周期为 Λ沿 x 方向排列z 方向是光的传播方向。矩形光栅在一个周期内的介电常数可以写成epsilon(x) n_air^2 (n_grating^2 - n_air^2) * rect(x / (duty * Λ))其中 duty 是占空比rect 表示矩形窗口函数。这个分段常数函数可以直接展开成 Fourier 级数epsilon(x) Σ eps_m * exp(j * 2π * m * x / Λ)求和范围取 m -N 到 N。对矩形光栅eps_m 有解析形式m0 时 eps_0 duty * n_grating^2 (1-duty) * n_air^2m≠0 时 eps_m (n_grating^2 - n_air^2) * sin(π * m * duty) / (π * m)。注意这里的“光栅常数”指的就是 Λ单位通常是纳米但在 MATLAB 里所有长度建议统一用米否则后面构造波矢时特征值数量级会很乱。对梯形光栅、正弦光栅或任意浮雕轮廓解析公式就不好使了。工程上更通用的做法是对一个周期内的折射率分布直接采样用离散 Fourier 变换得到 eps_m这与用解析公式构建的矩阵在 RCWA 计算流程里没有区别。很多现成 RCWA 程序只支持矩形光栅是因为解析公式写起来最快并不代表算法本身受限。2.2 空间谐波与 Floquet 波矢入射平面波的波矢为 k0 2π/λ横向分量为 k0 * n0 * sinθ。进入周期结构后第 m 级空间谐波的横向波矢由 Floquet 定理给出k_xm k0 * (n0 * sinθ - m * λ / Λ)这个负号对应的是衍射方程 mλ Λ * (n0*sinθ - sinθ_m) 的另一种写法。在 MATLAB 里构造横向波矢向量的写法是k0 2*pi / lambda; m -N:N; kx k0 * (n0 * sind(theta) - m * lambda / period);这里sind直接用角度制避免每次手动转弧度。m * lambda / period这一项体现了周期和波长之比对衍射级次分布的决定作用周期越小相邻级次的横向波矢间隔越大能传播的衍射级越少。这也是为什么亚波长光栅通常只有 0 级和 ±1 级能出现而粗光栅会看到一堆高级次。2.3 特征值方程与“严格”二字的来源把光栅层内的电磁场按空间谐波展开后代入 Maxwell 方程组TE 偏振下会得到一个矩阵形式的 Helmholtz 方程d^2 U / dz^2 (Kx^2 - E) * U其中 Kx 是 k_xm / k0 构成的对角阵E 是介电常数谐波构成的 Toeplitz 矩阵。这个方程的解写成 U(z) V * exp(q * z)代入后变成标准特征值问题eig(Kx^2 - E) q_i^2RCWA 的“严格”就在这里它没有对折射率差、刻蚀深度或占空比做任何小量近似误差来源只有 Fourier 级数的截断阶数 N。只要 N 取得足够大结果就收敛到 Maxwell 方程组的精确解。这一点与标量衍射理论有本质区别标量方法假设光栅周期远大于波长在亚波长尺度下会明显失真。2.4 多层结构、S 矩阵与 RCWA 的适用边界实际光栅刻蚀深度较大时需要把光栅在 z 方向切片每一层分别构建特征值问题再在层界面处匹配电磁场切向分量。早期程序多用传递矩阵法逐层传播但深光栅中 exp(q * depth) 会迅速溢出数值上非常不稳定。现代 RCWA 实现普遍改成 S 矩阵把每层看成二端口网络逐层散射叠加这是让 RCWA 在深宽比 10:1 以上的结构里仍然可用的关键。对比项RCWAFDTD周期结构适配天然满足需要额外设置周期边界内存瓶颈M^2 量级M2N1网格数 × 时间步数斜入射处理修改 k_xm 即可需要处理斜入射周期边界非周期结构不适用加 PML 即可典型单点耗时秒级到分钟级分钟级到小时级这张表不是要说明 RCWA 全面优于 FDTD而是明确它的定位周期结构、平面波入射、关注衍射效率分布这三条同时满足时RCWA 是第一选择。3. 手写 RCWA 的 MATLAB 骨架从参数输入到效率输出3.1 先给一组能跑的物理参数以一个典型的矩形透射光栅为例周期 500 nm入射波长 633 nm占空比 0.5刻蚀深度 250 nm光栅材料折射率 1.5衬底同为 1.5入射角 0 度TE 偏振。把这组参数写进 MATLAB 的参数字典里% rcwa_demo_params.m lambda 633e-9; % 入射波长, m period 500e-9; % 光栅周期, 也就是光栅常数 duty 0.5; % 占空比, 光栅脊宽/周期 depth 250e-9; % 刻蚀深度, m n0 1.0; % 入射介质折射率, 空气 n1 1.5; % 光栅材料折射率 n2 1.5; % 衬底折射率, 这里假设与光栅同材料 theta 0; % 入射角, 单位度 N 15; % 谐波截断阶数, 总级次为 2*N1 TE true; % trueTE, falseTM参数说明N 是 RCWA 里最需要关注的量后面单独讲怎么选。duty对衍射效率分布影响极大0.5 时偶次谐波系数为零这是一个很方便的调试锚点。所有长度统一用米避免在特征值计算中出现量纲混乱。3.2 构建介电常数的 Toeplitz 卷积矩阵RCWA 的核心矩阵 E 是一个 Toeplitz 矩阵每一行由 eps_m 平移得到。MATLAB 里可以直接用toeplitz一行构建但要注意向量的排列顺序% rcwa_build_E.m m -N:N; eps_g n1^2; eps_b n0^2; % 计算矩形光栅的傅里叶系数 eps_m zeros(1, 2*N1); for idx 1:length(m) if m(idx) 0 eps_m(idx) duty*eps_g (1-duty)*eps_b; else eps_m(idx) (eps_g-eps_b) * sin(pi*m(idx)*duty) / (pi*m(idx)); end end % Toeplitz 矩阵化 E toeplitz(eps_m(N1:end), eps_m(N1:-1:1)); if ~TE % TM 偏振下使用 Lalanne 提出的倒数规则 E inv(toeplitz(1./eps_m(N1:end), 1./eps_m(N1:-1:1))); end逻辑说明toeplitz的第一个参数是第一列第二个参数是第一行。eps_m(N1:end)对应从 0 阶开始到最高阶eps_m(N1:-1:1)是从 0 阶到最低阶这样拼出来的矩阵才满足卷积结构而不是转置关系。TM 偏振下不能直接用 eps_m 的 Toeplitz 矩阵而要把 1/eps_m 展开成 Toeplitz 后取逆这是 RCWA 实现中最容易踩的坑直接决定 TM 效率曲线是否收敛。3.3 特征值求解与纵向传播常数有了 E 矩阵后特征值问题只剩下几行代码% rcwa_eig.m k0 2*pi / lambda; theta_rad theta * pi / 180; kx k0 * (n0*sin(theta_rad) - m*lambda/period); Kx diag(kx / k0); A Kx^2 - E; [V, Q2] eig(A); q sqrt(diag(Q2)) * k0; % 符号修正实部取正, 虚部取负, 保证场往正z方向衰减 for idx 1:length(q) if real(q(idx)) 0 q(idx) -q(idx); end if imag(q(idx)) 0 q(idx) -q(idx); end end参数说明Kx^2 - E是 TE 模的方程TM 模只需要换掉 E 的构建方式特征值求解代码完全复用。q的符号修正是关键细节如果省掉边界匹配矩阵会得到指数增长的虚假解效率可能出现负值。MATLAB 的eig输出的特征值顺序不固定但特征向量矩阵 V 与特征值是一一对应的后续矩阵装配必须使用同样的列顺序。3.4 边界匹配和衍射效率输出特征解本身不代表效率效率来自边界匹配。以单层 TE 光栅为例在 z0 和 zdepth 界面处匹配电场和磁场的切向分量得到线性方程组% rcwa_match_bc.m 示意, 完整实现建议使用 S 矩阵 Zg diag(q / k0); M [V V; V*Zg -V*Zg]; % 光栅层内的基函数及其导数 b [eye(2*N1); diag(1i * (kx - k0*n0*cos(theta_rad)) / k0)]; coef M \ b; Rt coef(1:2*N1, 1:2*N1); % 反射振幅 Tt coef(2*N1:end, 1:2*N1); % 透射振幅边界匹配是 RCWA 里最容易出错的部分不同教材对基函数排列顺序的定义不同正负号一错效率就会乱。我一般建议新手不要从零写这一步直接找一份成熟的 S 矩阵实现作为参考把自己构建的 E 和 q 塞进去即可。效率计算到最后一步还要按 Poynting 矢量归一化DE_m |Tt(m)|^2 * real(q_sub / (k0 * n0 * cos(theta_rad)))这里的 q_sub 是衬底里的纵向波矢不能直接用光栅层的 q否则透射效率之和不会守恒。3.5 输出自检表跑通第一版后建议把以下三项检查放在脚本末尾检查项判断方法能量守恒sum(反射效率)sum(透射效率) 应等于 1 ± 1e-3对称性占空比 0.5 垂直入射时m 与 -m 级效率应相等零级趋势占空比扫描时 0 级效率变化应平滑无突变尖峰如果能量不守恒先怀疑截断阶数不够再检查边界匹配矩阵符号最后检查 Toeplitz 的列顺序。这个排查顺序能解决九成以上的 RCWA 数值问题。4. 光栅衍射效率的参数整定截断阶数、占空比与偏振4.1 衍射级截断阶数 N 怎么取N 选小了高级次衍射级贡献被截断效率偏低N 选大了矩阵规模按 (2N1)^2 增长单点计算从毫秒级涨到秒级。经验上可以按周期与波长之比来定初始值周期/波长 Λ/λ推荐 N理由 0.810~15亚波长区高阶谐波衰减快0.8 ~ 220~30有多个传播衍射级参与能量分配 240 以上高级次衍射级仍携带可观能量这个表不是硬性规定RCWA 只有一条判断标准效率随 N 增大不再显著变化。实际操作中我会从 N15 起步跑一次后把 N 翻倍再跑两次结果的同一级次效率差小于 0.001 就继续用较小的 N。注意衍射级数多不代表 N 必须大真正决定 N 的是折射率对比度和刻蚀深度高折射率差结构需要的 N 明显更多。4.2 占空比和光栅周期对效率分布的影响占空比改变的是 eps_m 的包络函数。占空比 0.5 时sin(π * m * 0.5) 在偶数 m 处取零因此所有偶次衍射级消失这是验证程序是否正确的一个很省力的特征如果垂直入射、占空比 0.5 时算出了 2 级衍射效率说明傅里叶系数有问题。周期对效率的影响更直接。固定波长时周期越大传播衍射级越多能量分散到更多级次中单级效率通常下降。做优化时我一般把周期和占空比放在同一个二维扫描里先粗扫定位峰值区域再在峰值附近加密扫描而不是单独调一个变量。RCWA 单点计算成本低这种参数扫描在 MATLAB 里跑几百个点也就几十秒。4.3 深光栅和高折射率差的数值稳定处理深度增加时exp(q * depth) 中的大实数会指数增长直接构造传递矩阵很快就溢出。常见做法是每层都用 S 矩阵级联把指数因子限制在 [0,1] 范围内。MATLAB 里判断数值是否损坏的快捷方法是看透射效率是否出现负值负效率几乎都是特征值符号选错或 S 矩阵拼接顺序错误。另外一个容易忽略的问题是层数划分。对矩形光栅单层模型在解析上是严格的但如果占空比沿深度变化比如梯形轮廓至少要切 10 层以上。切层不是越多越好RCWA 每层都要做一次特征值分解50 层的计算量已经比较可观工程上通常先做收敛性测试确定最小层数。4.4 斜入射与偏振切换时的参数整定斜入射时入射角 θ 进入 k_xm 公式整个衍射级次的横向波矢不再对称1 级和 -1 级的效率开始分化。这个现象在光栅光谱仪设计中是核心指标RCWA 里只需要改theta一个参数不需要动其他结构。注意斜入射时同一个 N 下需要的谐波截断可能更多因为衍射级次数量随入射角变化建议在斜入射计算里直接用比垂直入射高 50% 的 N。偏振切换时TE 和 TM 的 E 矩阵构建方式不同这是 RCWA 程序里最容易写错的地方。TE 用 eps_m 的 Toeplitz 直接参与特征值方程TM 必须用 1/eps_m 展开再取逆。如果只追求快速定性结果可以写两个分支if TE E toeplitz(eps_m(N1:end), eps_m(N1:-1:1)); else E inv(toeplitz(1./eps_m(N1:end), 1./eps_m(N1:-1:1))); end对比 TE 和 TM 效率曲线是判断程序是否正确的另一个手段在深亚波长光栅里两者差异往往超过 20 个百分点如果算出来几乎一样程序大概率有 bug。4.5 与实验值、FDTD 结果对标时先检查什么RCWA 仿真和实验值对不上时先怀疑的不是算法而是模型。实际光栅的侧壁角不是 90 度底部有圆角刻蚀深度有偏差这些几何误差对效率的影响经常比理论精度大一个量级。用分光计或 SEM 标定的光栅周期和实际浮雕轮廓远比把 RCWA 的 N 从 30 调到 50 更能缩小偏差。和 FDTD 对标的场景则相反两者用同一组矩形参数时差异应该极小如果差异超过 2%优先检查 RCWA 的截断阶数和 FDTD 的网格精度而不是模型偏差。5. 收敛性检查的技巧能量守恒偏差与截断扫描RCWA 结果可不可信不看程序跑没跑通而看收敛性有没有被验证。最常用的硬指标有两个能量守恒偏差和效率随 N 的收敛曲线。能量守恒是所有无损耗光栅必须满足的条件MATLAB 里用一行脚本就能算energy_balance sum(DE_reflect) sum(DE_transmit) - 1;DE_reflect和DE_transmit是所有传播衍射级的效率向量。能量偏差在 1e-3 量级可以接受超过 1e-2 说明 N 不足或是边界匹配有误。注意如果光栅材料有吸收能量守恒公式要写成 sum(DE_reflect) sum(DE_transmit) absorption 1千万不要直接拿透反射效率加和去判断。第二个技巧是固定物理参数只扫描 NN_list 5:5:40; DE0 zeros(size(N_list)); DE1 zeros(size(N_list)); for k 1:length(N_list) DE rcwa_1d(633e-9, 500e-9, 0.5, 250e-9, ... 1.5, 1.5, 0, TE, N_list(k)); DE0(k) DE(0); DE1(k) DE(1); % 这里按代码内部约定取 1 级 end plot(N_list, DE0, -o, N_list, DE1, -s); xlabel(截断阶数 N); ylabel(衍射效率); legend(0级, 1级);把 0 级和 1 级效率画在同一张图上最直观的判断标准是曲线的尾部变成水平线或者在小范围内振荡且振幅小于 0.001。如果 N 从 5 加到 40 效率还在单调漂移说明问题不止是截断要回头检查介电常数矩阵构建逻辑。N 的扫描间隔用奇数偶数无所谓但建议从 N5 这种很小的值开始能看到明显的欠收敛状态反而有助于理解截断误差的特征。最后一个对日常调试很实用的技巧把 rcwa_1d 封装成纯函数后在参数扫描循环里对每个参数点都自动做一次能量守恒校验任何异常点立刻打警告。做参数优化时异常值往往不是优化算法的问题而是 RCWA 本身在某个占空比或深度下收敛变慢这时单独对该点加密 N 即可不用改全局设置。把上述收敛核验命令写进一个convergence_test.m下次跑参数扫描之前先执行一遍比直接信任默认 N 值省下来的调试时间往往不止一小时。本文还有配套的精品资源点击获取 MATLAB实现光栅衍射效率的严格耦合波分析
简介RCWA严格耦合波分析是光栅衍射特性仿真中的核心数值方法这份资源以MATLAB实现了一维光栅TE模式下的RCWA计算程序适合从事微纳光学、衍射光学元件设计的研究人员或相关专业学生用于快速求解光栅的衍射效率并验证理论模型。压缩包共10个文件其中6个.m脚本可实现主程序RCWA.m、特征值计算eig_M.m、波矢求解kvect.m/kvect1.m与介电常数转换eps_conv.m等核心功能另有4个.asv为MATLAB自动备份版本便于追溯修改过程整包仅5KB体量轻巧、便于直接读取和二次开发。目前已有1151人学习/下载。通过该程序可计算不同入射角与波长下各衍射级的效率并能结合MATLAB的矩阵运算与绘图功能分析TE偏振光在一维光栅中的传播行为适合作为光栅严格耦合波分析的入门模板或教学示例也可在此基础上扩展至二维光栅及TM模式研究。1. 光栅衍射效率计算为什么绕不开严格耦合波分析做衍射光栅的工程师手里大概率囤过不止一个 RCWA 压缩包。类似“RCWA程序.rar”这样的文件很多是从老旧的下载站转存来的打开之后是命名混乱的 m 文件和几张 demo 图。RCWA 的全称是 Rigorous Coupled-Wave Analysis中文语境里通常叫严格耦合波分析它不是一个软件而是一类算法给定周期结构与入射平面波直接求出 0 级、±1 级、±2 级等各衍射级次的效率与相位。FDTD 那种在网格上逐步推进的方式当然也能算但亚波长光栅按纳米尺度加密网格后内存和时间开销会迅速失控RCWA 利用周期性做谱展开把麦克斯韦方程压缩成矩阵特征值问题计算速度往往快一到两个数量级。这篇文章就围绕标题里的四个关键词——MATLAB、光栅、严格耦合、衍射效率把 RCWA 的原理、参数和 MATLAB 落地路径写清楚。2. 严格耦合波分析建模从周期结构到特征值问题2.1 周期性光栅与介电常数谐波展开光栅本质上是一个空间周期性变化的介电常数结构。设光栅周期为 Λ沿 x 方向排列z 方向是光的传播方向。矩形光栅在一个周期内的介电常数可以写成epsilon(x) n_air^2 (n_grating^2 - n_air^2) * rect(x / (duty * Λ))其中 duty 是占空比rect 表示矩形窗口函数。这个分段常数函数可以直接展开成 Fourier 级数epsilon(x) Σ eps_m * exp(j * 2π * m * x / Λ)求和范围取 m -N 到 N。对矩形光栅eps_m 有解析形式m0 时 eps_0 duty * n_grating^2 (1-duty) * n_air^2m≠0 时 eps_m (n_grating^2 - n_air^2) * sin(π * m * duty) / (π * m)。注意这里的“光栅常数”指的就是 Λ单位通常是纳米但在 MATLAB 里所有长度建议统一用米否则后面构造波矢时特征值数量级会很乱。对梯形光栅、正弦光栅或任意浮雕轮廓解析公式就不好使了。工程上更通用的做法是对一个周期内的折射率分布直接采样用离散 Fourier 变换得到 eps_m这与用解析公式构建的矩阵在 RCWA 计算流程里没有区别。很多现成 RCWA 程序只支持矩形光栅是因为解析公式写起来最快并不代表算法本身受限。2.2 空间谐波与 Floquet 波矢入射平面波的波矢为 k0 2π/λ横向分量为 k0 * n0 * sinθ。进入周期结构后第 m 级空间谐波的横向波矢由 Floquet 定理给出k_xm k0 * (n0 * sinθ - m * λ / Λ)这个负号对应的是衍射方程 mλ Λ * (n0*sinθ - sinθ_m) 的另一种写法。在 MATLAB 里构造横向波矢向量的写法是k0 2*pi / lambda; m -N:N; kx k0 * (n0 * sind(theta) - m * lambda / period);这里sind直接用角度制避免每次手动转弧度。m * lambda / period这一项体现了周期和波长之比对衍射级次分布的决定作用周期越小相邻级次的横向波矢间隔越大能传播的衍射级越少。这也是为什么亚波长光栅通常只有 0 级和 ±1 级能出现而粗光栅会看到一堆高级次。2.3 特征值方程与“严格”二字的来源把光栅层内的电磁场按空间谐波展开后代入 Maxwell 方程组TE 偏振下会得到一个矩阵形式的 Helmholtz 方程d^2 U / dz^2 (Kx^2 - E) * U其中 Kx 是 k_xm / k0 构成的对角阵E 是介电常数谐波构成的 Toeplitz 矩阵。这个方程的解写成 U(z) V * exp(q * z)代入后变成标准特征值问题eig(Kx^2 - E) q_i^2RCWA 的“严格”就在这里它没有对折射率差、刻蚀深度或占空比做任何小量近似误差来源只有 Fourier 级数的截断阶数 N。只要 N 取得足够大结果就收敛到 Maxwell 方程组的精确解。这一点与标量衍射理论有本质区别标量方法假设光栅周期远大于波长在亚波长尺度下会明显失真。2.4 多层结构、S 矩阵与 RCWA 的适用边界实际光栅刻蚀深度较大时需要把光栅在 z 方向切片每一层分别构建特征值问题再在层界面处匹配电磁场切向分量。早期程序多用传递矩阵法逐层传播但深光栅中 exp(q * depth) 会迅速溢出数值上非常不稳定。现代 RCWA 实现普遍改成 S 矩阵把每层看成二端口网络逐层散射叠加这是让 RCWA 在深宽比 10:1 以上的结构里仍然可用的关键。对比项RCWAFDTD周期结构适配天然满足需要额外设置周期边界内存瓶颈M^2 量级M2N1网格数 × 时间步数斜入射处理修改 k_xm 即可需要处理斜入射周期边界非周期结构不适用加 PML 即可典型单点耗时秒级到分钟级分钟级到小时级这张表不是要说明 RCWA 全面优于 FDTD而是明确它的定位周期结构、平面波入射、关注衍射效率分布这三条同时满足时RCWA 是第一选择。3. 手写 RCWA 的 MATLAB 骨架从参数输入到效率输出3.1 先给一组能跑的物理参数以一个典型的矩形透射光栅为例周期 500 nm入射波长 633 nm占空比 0.5刻蚀深度 250 nm光栅材料折射率 1.5衬底同为 1.5入射角 0 度TE 偏振。把这组参数写进 MATLAB 的参数字典里% rcwa_demo_params.m lambda 633e-9; % 入射波长, m period 500e-9; % 光栅周期, 也就是光栅常数 duty 0.5; % 占空比, 光栅脊宽/周期 depth 250e-9; % 刻蚀深度, m n0 1.0; % 入射介质折射率, 空气 n1 1.5; % 光栅材料折射率 n2 1.5; % 衬底折射率, 这里假设与光栅同材料 theta 0; % 入射角, 单位度 N 15; % 谐波截断阶数, 总级次为 2*N1 TE true; % trueTE, falseTM参数说明N 是 RCWA 里最需要关注的量后面单独讲怎么选。duty对衍射效率分布影响极大0.5 时偶次谐波系数为零这是一个很方便的调试锚点。所有长度统一用米避免在特征值计算中出现量纲混乱。3.2 构建介电常数的 Toeplitz 卷积矩阵RCWA 的核心矩阵 E 是一个 Toeplitz 矩阵每一行由 eps_m 平移得到。MATLAB 里可以直接用toeplitz一行构建但要注意向量的排列顺序% rcwa_build_E.m m -N:N; eps_g n1^2; eps_b n0^2; % 计算矩形光栅的傅里叶系数 eps_m zeros(1, 2*N1); for idx 1:length(m) if m(idx) 0 eps_m(idx) duty*eps_g (1-duty)*eps_b; else eps_m(idx) (eps_g-eps_b) * sin(pi*m(idx)*duty) / (pi*m(idx)); end end % Toeplitz 矩阵化 E toeplitz(eps_m(N1:end), eps_m(N1:-1:1)); if ~TE % TM 偏振下使用 Lalanne 提出的倒数规则 E inv(toeplitz(1./eps_m(N1:end), 1./eps_m(N1:-1:1))); end逻辑说明toeplitz的第一个参数是第一列第二个参数是第一行。eps_m(N1:end)对应从 0 阶开始到最高阶eps_m(N1:-1:1)是从 0 阶到最低阶这样拼出来的矩阵才满足卷积结构而不是转置关系。TM 偏振下不能直接用 eps_m 的 Toeplitz 矩阵而要把 1/eps_m 展开成 Toeplitz 后取逆这是 RCWA 实现中最容易踩的坑直接决定 TM 效率曲线是否收敛。3.3 特征值求解与纵向传播常数有了 E 矩阵后特征值问题只剩下几行代码% rcwa_eig.m k0 2*pi / lambda; theta_rad theta * pi / 180; kx k0 * (n0*sin(theta_rad) - m*lambda/period); Kx diag(kx / k0); A Kx^2 - E; [V, Q2] eig(A); q sqrt(diag(Q2)) * k0; % 符号修正实部取正, 虚部取负, 保证场往正z方向衰减 for idx 1:length(q) if real(q(idx)) 0 q(idx) -q(idx); end if imag(q(idx)) 0 q(idx) -q(idx); end end参数说明Kx^2 - E是 TE 模的方程TM 模只需要换掉 E 的构建方式特征值求解代码完全复用。q的符号修正是关键细节如果省掉边界匹配矩阵会得到指数增长的虚假解效率可能出现负值。MATLAB 的eig输出的特征值顺序不固定但特征向量矩阵 V 与特征值是一一对应的后续矩阵装配必须使用同样的列顺序。3.4 边界匹配和衍射效率输出特征解本身不代表效率效率来自边界匹配。以单层 TE 光栅为例在 z0 和 zdepth 界面处匹配电场和磁场的切向分量得到线性方程组% rcwa_match_bc.m 示意, 完整实现建议使用 S 矩阵 Zg diag(q / k0); M [V V; V*Zg -V*Zg]; % 光栅层内的基函数及其导数 b [eye(2*N1); diag(1i * (kx - k0*n0*cos(theta_rad)) / k0)]; coef M \ b; Rt coef(1:2*N1, 1:2*N1); % 反射振幅 Tt coef(2*N1:end, 1:2*N1); % 透射振幅边界匹配是 RCWA 里最容易出错的部分不同教材对基函数排列顺序的定义不同正负号一错效率就会乱。我一般建议新手不要从零写这一步直接找一份成熟的 S 矩阵实现作为参考把自己构建的 E 和 q 塞进去即可。效率计算到最后一步还要按 Poynting 矢量归一化DE_m |Tt(m)|^2 * real(q_sub / (k0 * n0 * cos(theta_rad)))这里的 q_sub 是衬底里的纵向波矢不能直接用光栅层的 q否则透射效率之和不会守恒。3.5 输出自检表跑通第一版后建议把以下三项检查放在脚本末尾检查项判断方法能量守恒sum(反射效率)sum(透射效率) 应等于 1 ± 1e-3对称性占空比 0.5 垂直入射时m 与 -m 级效率应相等零级趋势占空比扫描时 0 级效率变化应平滑无突变尖峰如果能量不守恒先怀疑截断阶数不够再检查边界匹配矩阵符号最后检查 Toeplitz 的列顺序。这个排查顺序能解决九成以上的 RCWA 数值问题。4. 光栅衍射效率的参数整定截断阶数、占空比与偏振4.1 衍射级截断阶数 N 怎么取N 选小了高级次衍射级贡献被截断效率偏低N 选大了矩阵规模按 (2N1)^2 增长单点计算从毫秒级涨到秒级。经验上可以按周期与波长之比来定初始值周期/波长 Λ/λ推荐 N理由 0.810~15亚波长区高阶谐波衰减快0.8 ~ 220~30有多个传播衍射级参与能量分配 240 以上高级次衍射级仍携带可观能量这个表不是硬性规定RCWA 只有一条判断标准效率随 N 增大不再显著变化。实际操作中我会从 N15 起步跑一次后把 N 翻倍再跑两次结果的同一级次效率差小于 0.001 就继续用较小的 N。注意衍射级数多不代表 N 必须大真正决定 N 的是折射率对比度和刻蚀深度高折射率差结构需要的 N 明显更多。4.2 占空比和光栅周期对效率分布的影响占空比改变的是 eps_m 的包络函数。占空比 0.5 时sin(π * m * 0.5) 在偶数 m 处取零因此所有偶次衍射级消失这是验证程序是否正确的一个很省力的特征如果垂直入射、占空比 0.5 时算出了 2 级衍射效率说明傅里叶系数有问题。周期对效率的影响更直接。固定波长时周期越大传播衍射级越多能量分散到更多级次中单级效率通常下降。做优化时我一般把周期和占空比放在同一个二维扫描里先粗扫定位峰值区域再在峰值附近加密扫描而不是单独调一个变量。RCWA 单点计算成本低这种参数扫描在 MATLAB 里跑几百个点也就几十秒。4.3 深光栅和高折射率差的数值稳定处理深度增加时exp(q * depth) 中的大实数会指数增长直接构造传递矩阵很快就溢出。常见做法是每层都用 S 矩阵级联把指数因子限制在 [0,1] 范围内。MATLAB 里判断数值是否损坏的快捷方法是看透射效率是否出现负值负效率几乎都是特征值符号选错或 S 矩阵拼接顺序错误。另外一个容易忽略的问题是层数划分。对矩形光栅单层模型在解析上是严格的但如果占空比沿深度变化比如梯形轮廓至少要切 10 层以上。切层不是越多越好RCWA 每层都要做一次特征值分解50 层的计算量已经比较可观工程上通常先做收敛性测试确定最小层数。4.4 斜入射与偏振切换时的参数整定斜入射时入射角 θ 进入 k_xm 公式整个衍射级次的横向波矢不再对称1 级和 -1 级的效率开始分化。这个现象在光栅光谱仪设计中是核心指标RCWA 里只需要改theta一个参数不需要动其他结构。注意斜入射时同一个 N 下需要的谐波截断可能更多因为衍射级次数量随入射角变化建议在斜入射计算里直接用比垂直入射高 50% 的 N。偏振切换时TE 和 TM 的 E 矩阵构建方式不同这是 RCWA 程序里最容易写错的地方。TE 用 eps_m 的 Toeplitz 直接参与特征值方程TM 必须用 1/eps_m 展开再取逆。如果只追求快速定性结果可以写两个分支if TE E toeplitz(eps_m(N1:end), eps_m(N1:-1:1)); else E inv(toeplitz(1./eps_m(N1:end), 1./eps_m(N1:-1:1))); end对比 TE 和 TM 效率曲线是判断程序是否正确的另一个手段在深亚波长光栅里两者差异往往超过 20 个百分点如果算出来几乎一样程序大概率有 bug。4.5 与实验值、FDTD 结果对标时先检查什么RCWA 仿真和实验值对不上时先怀疑的不是算法而是模型。实际光栅的侧壁角不是 90 度底部有圆角刻蚀深度有偏差这些几何误差对效率的影响经常比理论精度大一个量级。用分光计或 SEM 标定的光栅周期和实际浮雕轮廓远比把 RCWA 的 N 从 30 调到 50 更能缩小偏差。和 FDTD 对标的场景则相反两者用同一组矩形参数时差异应该极小如果差异超过 2%优先检查 RCWA 的截断阶数和 FDTD 的网格精度而不是模型偏差。5. 收敛性检查的技巧能量守恒偏差与截断扫描RCWA 结果可不可信不看程序跑没跑通而看收敛性有没有被验证。最常用的硬指标有两个能量守恒偏差和效率随 N 的收敛曲线。能量守恒是所有无损耗光栅必须满足的条件MATLAB 里用一行脚本就能算energy_balance sum(DE_reflect) sum(DE_transmit) - 1;DE_reflect和DE_transmit是所有传播衍射级的效率向量。能量偏差在 1e-3 量级可以接受超过 1e-2 说明 N 不足或是边界匹配有误。注意如果光栅材料有吸收能量守恒公式要写成 sum(DE_reflect) sum(DE_transmit) absorption 1千万不要直接拿透反射效率加和去判断。第二个技巧是固定物理参数只扫描 NN_list 5:5:40; DE0 zeros(size(N_list)); DE1 zeros(size(N_list)); for k 1:length(N_list) DE rcwa_1d(633e-9, 500e-9, 0.5, 250e-9, ... 1.5, 1.5, 0, TE, N_list(k)); DE0(k) DE(0); DE1(k) DE(1); % 这里按代码内部约定取 1 级 end plot(N_list, DE0, -o, N_list, DE1, -s); xlabel(截断阶数 N); ylabel(衍射效率); legend(0级, 1级);把 0 级和 1 级效率画在同一张图上最直观的判断标准是曲线的尾部变成水平线或者在小范围内振荡且振幅小于 0.001。如果 N 从 5 加到 40 效率还在单调漂移说明问题不止是截断要回头检查介电常数矩阵构建逻辑。N 的扫描间隔用奇数偶数无所谓但建议从 N5 这种很小的值开始能看到明显的欠收敛状态反而有助于理解截断误差的特征。最后一个对日常调试很实用的技巧把 rcwa_1d 封装成纯函数后在参数扫描循环里对每个参数点都自动做一次能量守恒校验任何异常点立刻打警告。做参数优化时异常值往往不是优化算法的问题而是 RCWA 本身在某个占空比或深度下收敛变慢这时单独对该点加密 N 即可不用改全局设置。把上述收敛核验命令写进一个convergence_test.m下次跑参数扫描之前先执行一遍比直接信任默认 N 值省下来的调试时间往往不止一小时。本文还有配套的精品资源点击获取