严格耦合波分析（RCWA）如何准确计算一维光栅衍射效率？

简介：RCWA（严格耦合波分析）计算一维光栅衍射效率的资源包，面向光学与光子学领域的研究人员、工程师及高年级学生，适用于光谱仪、激光器、光通信等场景中的光栅设计与性能评估。压缩包共2个文件，包含1个MATLAB脚本和1篇caj格式的学术论文，整体大小仅4.84MB。MATLAB脚本覆盖了设置光栅周期、填充因子、材料折射率、入射波长与极化状态等参数，并输出各衍射级效率的完整求解流程；论文则从理论层面系统阐述RCWA在相位光栅衍射特性研究中的具体应用，包括方程推导和案例分析，帮助读者在动手实践的同时深入理解数值结果背后的物理机理。相比时域有限差分等方法，RCWA在分析和优化一维光栅结构时更具计算效率，尤其适用于高阶衍射模式。目前已有1945人学习下载，对优化光栅参数、提升器件性能具有直接参考价值。

1. 严格耦合波分析：为什么一维光栅的衍射效率不能靠标量近似来算

我之前做一维光栅衍射效率，最开始直接用光栅方程算出衍射角，再用标量衍射积分给各级效率。周期 5 微米、深度 200 纳米的感觉“大光栅”似乎够了，但换成 TM 偏振的金属光栅，平板仿真和实验对不上，效率差了 30% 以上，还出现一种“幽灵峰”，让所有级次能量总和超过 1。后来才意识到：标量近似把每个周期当成孤立散射中心，忽略了光栅内部的多次散射和各级之间的相干。RCWA 用傅里叶级数展开周期性介电常数，把麦克斯韦方程组变成可解的代数问题，给出的“透射系数”才是可以比较的物理效率。这个算法不是黑匣子，而是一个把波动方程沿周期方向降维的框架。

2. 从波动方程到耦合波方程组：光栅区如何被“降维”成一个矩阵

2.1 标量近似的边界：什么条件下计算开始“翻车”

光栅方程给角度，不给定效率。标量衍射积分用相位调制 exp(j φ(x)) 直接算远场，意味着假设光栅深度很浅，光只经过一次，散射体彼此独立。这个假设在光栅深度超过 0.1λ 时就开始不可靠：光在槽内多次反射，传播衍射级与倏逝级之间互相转化，而标量模型里根本没有倏逝级的位置。判断是否需要用严格办法，最简单的经验是看相对折射率差：空气-硅光栅在近红外下折射率差超过 2，即使深度只有 λ/10，TM 偏振下标量近似的误差也是肉眼可见的；换成金属光栅，误差可以到百分之几十。另一个信号是输入不同偏振时结果几乎没有区别——真光栅一定偏振敏感，标量模型却没有区分偏振的能力，这个现象一出，基本可以断定当前计算模型已经在硬撑。

2.2 空间谐波展开：周期性介电常数怎么进入矩阵

RCWA 的核心是区域分解：把空间分成入射均匀区、光栅周期性区和透射均匀区，只在光栅区展开场，在均匀区保留平面波。一维光栅沿 x 方向周期为 Λ，介电常数满足 ε(x+Λ)=ε(x)，于是把 ε(x) 写成傅里叶级数：

ε(x) = Σ_h ε_h exp(j h K x)，其中 K = 2π/Λ

入射区域里，各级平面波的切向波矢满足 Floquet 条件：

k_x,m = k₀ sinθ - m K

这个 m 被称为衍射级，但在光栅内部它同时也是场的一个空间谐波。光栅区不是均匀介质，这些波之间会发生“耦合”：第 m 级方程里会出现所有 h 阶介电常数展开与第 m-h 级场的乘积，这正是“严格耦合波”里“耦合波”三个字的来源。若研究的是一维矩形光栅，介电常数是占空比为 duty 的两值分布，傅里叶系数能用闭式表达：

ε₀ = f·ε_g + (1-f)·ε_a，ε_h = (ε_g - ε_a)·sin(π h f)/(π h)，h≠0

由这些系数构造 Toeplitz 矩阵 [ε]，元素只依赖于行号与列号之差。周期性介电常数进入方程之后，求解对象从偏微分方程变成常微分方程组，维度就是截断后的衍射级数。下面代码展示了这一步对应的具体操作。

# 只演示核心矩阵构建，其余边界匹配见后续章节 import numpy as np order = 15 # 截断阶数，总级数 2*order+1 m = np.arange(-order, order + 1) period = 1.2e-3 # 光栅周期，单位 mm lam = 0.633e-3 # 入射波长，单位 mm theta = 0.0 # 入射角，正入射 duty = 0.5 # 光栅占空比 k0 = 2 * np.pi / lam K = 2 * np.pi / period kx_m = k0 * np.sin(theta) + m * K # 各空间谐波的切向波矢 def eps_hat(h, eps_g, eps_a, f): """矩形光栅介电常数的傅里叶系数""" if h == 0: return f * eps_g + (1 - f) * eps_a return (eps_g - eps_a) * np.sin(np.pi * h * f) / (np.pi * h) eps_toep = np.zeros((len(m), len(m)), dtype=complex) # Toeplitz 矩阵：行对应第 m 级方程，列对应第 h 级分量 for i in range(len(m)): for j in range(len(m)): eps_toep[i, j] = eps_hat(m[i] - m[j], 12.25, 1.0, duty) Kx = np.diag(kx_m / k0) # 归一化切向波矢 A_TE = Kx @ Kx - eps_toep # TE 偏振层内矩阵 val, vec = np.linalg.eig(A_TE) # 特征值即传播常数平方

这段代码里，eps_hat 返回的 h=0 项是光栅区和空气的等效平均介电常数，h≠0 项描述周期调制强度。注意当 duty=0.5 时只有偶次谐波强度随 h 衰减，而奇次项不为零，这正是方波傅里叶展开的典型行为。A_TE 的维度是 31×31，特征分解直接给出一层内可传播模式的分布。实际工程中应预先检查矩阵条件数，若最大最小奇异值相差超过 1e12，说明截断阶数选择或归一化有问题。

2.3 TE 与 TM 两套方程：差别在连续性边界，不在“是否严格”

TE 偏振下只需求 E_y 分量，光栅区的亥姆霍兹方程代入展开后得到向量微分方程：

S''(z) = (K_x² - k₀²[ε]) S(z)

这正是上节中 A_TE 的来源。矩阵化以后，问题退化为求 A 的特征值和特征向量。TM 偏振则不同：E_x 和 H_y 并存，而 E_x 在介电常数突变的界面上必须满足切向连续条件。如果仍直接用 ε 的 Toeplitz 矩阵，切向连续性条件会随着傅里叶截断产生寄生振荡，这就是“TM 收敛慢”的经典原因。标准的工程处理是用 1/ε 的傅里叶系数参与构造，也就是把介电常数矩阵换成 ε 的逆 Toeplitz 形式，再用公式：

P = K_x · [1/ε] · K_x - k₀² I

构造 TM 层的特征系统。教材里常把“TE 快、TM 慢”归因于场分量数量不同，实质是 TM 必须同时满足 E_x 和 H_y 的连续性，而这两个量在界面处对 ε 不连续的敏感程度完全不一样。到这里，光栅区已被转化为若干个常微分方程系统，下一步就看这些模式如何穿过一层一层的光栅。

3. 分层与 S 矩阵：决定结果是否“能跑出稳定数字”的环节

3.1 任意轮廓的阶梯化：分层数目怎么选

一维光栅真正只有矩形槽时才可以用单层模型。只要槽壁是斜面、正弦轮廓或刻有圆角，就必须沿 z 方向把光栅切成若干薄层，每一层里介电常数只随 x 变化、与 z 无关。这一近似把任意轮廓变成一组阶梯，层数越多越接近原轮廓，但矩阵尺寸并不会随层数增大，只是界面匹配次数增加。实际选择有两个约束：一是每个台阶在 z 方向的厚度要小到能让轮廓误差小于 λ/20；二是台阶沿 z 方向的边不能正好落在入射光干涉极值的位置，否则会产生人为共振。对倾斜角比较大的三角光栅，常见做法是先把轮廓按等弧长采样，再在每层内取平均介电常数，而不是简单用每个采样点的瞬时值。

3.2 层内特征模式与场展开

在每一层内，上文得到的矩阵 A 的特征值就是层内波导模式的传播常数平方，特征向量则是该模式在各空间谐波上的振幅权重。TE 偏振下场可以写成：

E_y(x,z) = Σ_i Σ_m W_im [c_i⁺ exp(-k₀ q_i z) + c_i⁻ exp(k₀ q_i z)] exp(-j k_x,m x)

其中 q_i = sqrt(λ_i)，取实部为正的支。q_i² 为负时 q_i 是虚数，对应沿 z 方向振荡的传播模式；q_i² 为正时 q_i 是实数，对应衰减模式。衰减模式看起来不向远处走，但它在层与层的界面上仍然有场分布，并借耦合影响其他模式的振幅，所以矩阵里必须保留，不能当“无效”丢掉。

边界匹配只用两个连续性条件：TE 情形是 E_y 和 ∂E_y/∂z 跨界连续；TM 情形是 H_y 和 E_x 连续。每在一个层界面列方程，就会得到一组线性关系。把所有层叠加上去后，未知数就是各模式振幅。这里有一个常见误解：有人以为分层越多越精确，于是直接把每层厚度设成 1 纳米，结果层数上千后矩阵匹配误差反而累积得比轮廓误差还大。分层前应该先看轮廓的曲率变化，平面区域用厚层、拐角区域加密。

3.3 为什么 T 矩阵会爆掉，而 S 矩阵能活下来

层数一多，最常见的数值事故是结果出现 NaN 或 Infinity。初级版本会把所有界面方程组织成传递矩阵 T，用 T 把最底层的振幅一步步推送到最顶层。问题在于，模式振幅表达式中既有 exp(-k₀ q_i z) 也有 exp(+k₀ q_i z)。对于衰减模式，正向和反向两种取向在深度达到几个波长后，数量级可以相差 e^10 以上。T 矩阵要同时携带这两种量，底层小量会被放大，顶层大量会溢出，最终一个 31×31 的矩阵乘两三层就毁掉双精度浮点。

S 矩阵换了一个组织方式：它把每个界面或每层光栅看作一双端口，输入是上下两端的入射振幅，输出是上下两端的出射振幅。合并两个相邻端口时使用 Redheffer 星积，在表达式中避免显式计算 e^{+k₀ q_i d}，衰减模的贡献始终保持为有限值的反射/透射系数。对同一个深光栅，T 矩阵在深度超过 5λ 时几乎必然崩溃，而 S 矩阵可以安静地推到更大的深度范围。

# 子层散射矩阵的级联：Redheffer 星积示意 def star_product(SL, SR, n): # SL、SR 均为 2n×2n 分块矩阵，n 是衍射级数维度 I = np.eye(n) # 关键一步：用左除代替显式求逆，缓解条件数过大时的精度损失 D = np.linalg.solve(I - SL[1, 1] @ SR[0, 0], SL[1, 1]) S11 = SL[0, 0] + SL[0, 1] @ D @ SR[0, 0] @ SL[1, 0] S12 = SL[0, 1] @ D @ SR[0, 1] S21 = SR[1, 0] @ SL[1, 0] + SR[1, 0] @ SL[1, 1] @ D @ SR[0, 0] @ SL[1, 0] S22 = SR[1, 0] @ SL[1, 1] @ D @ SR[0, 1] + SR[1, 1] return np.block([[S11, S12], [S21, S22]])

这段代码里的 D 是整个级联中最容易数值翻车的部分：如果层内包含强衰减模式，I - SL[1,1]@SR[0,0] 的条件数会非常大。用 solve 而不是 inv，至少能把残差控制在浮点误差范围内。实际库实现里这一步还会做一次 QR 分解来进一步稳定，但核心逻辑就是上面这一点：把指数增长项隐藏在散射系数里，而不是显式传给下一层。

4. 衍射效率：从系数到可测功率，以及最容易被忽视的倏逝级

4.1 用 Poynting 通量归一化效率

RCWA 算完场后，光栅顶部和底部的均匀介质里每个衍射级就是一组平面波。要知道衍射效率，不能直接取振幅模平方，因为各级平面波沿 z 方向携带的功率并不相等。设第 m 级切向波矢为 k_x,m，在折射率 n 的均匀介质中

k_z,m = k₀·sqrt(n² - (k_x,m/k₀)²)

功率归一化因子就是 Re(k_z,m / k_z,0)。于是反射效率为：

η_m^ref = |R_m|² · Re(k_z,m^ref / k_z,0^ref)

透射效率同理，把 R_m 换成 T_m，并取透射区的 k_z 因子。所谓“衍射级在某个角度突然消失”，对应的是 k_z² ≤ 0：这个级成为倏逝级，不再向远场传能量，因此对探测面不可见；但它不会从物理模型里消失。所以 RCWA 的效率计算第一步不是算振幅，而是先筛一遍哪些级是传播级。

4.2 传播级、倏逝级与光栅方程的恢复

判断某个衍射级是否传播的判别条件回到 k_z² > 0，化简后就是常见的：

-n ≤ n_i sinθ - m λ/Λ ≤ n

在一维光栅中，满足这个条件的 m 是有限个；其余 m 全部是倏逝级，它们在光栅边界处有指数衰减的场。常见实现里会把衰减模式取 Re(q)>0，保证衰减方向正确，不会因为选了相反的支导致能量骤然变负。这一点也解释了为什么“光栅方程给的角度数量和实验里看到的亮级数完全一致”：亮级就是特征值落在传播区间那几档。RCWA 并不依靠光栅方程来“猜”级次，而是在解完矩阵后用 k_z 的实际符号来区分。

4.3 收敛性与能量守恒：先让计算“自洽”再谈准确性

RCWA 是把无穷维的场展开截断成 2N+1 个谐波，N 是收敛参数。工程中一般从 N=21 起步，逐步加到 41、61、81，看目标衍射效率的变化；当变化小于 0.1% 时认为该级达标。另外一种更全局的判断是能量守恒：在入射介质和透射介质都不吸收的情况下，所有传播级反射效率之和加透射效率之和应等于 1，误差通常在 1e-6 量级。若算出来总和是 1.03，说明要么截断阶数不足，要么功率归一化因子用错了介质折射率，要么在边界匹配时把倏逝级的场振幅也当成功率流来数了。

TM 偏振下这个收敛测试尤其关键：TE 往往 N=31 就平了，TM 却要到 N=61 甚至更高。若你发现 TM 的结果在 N=81 时还在缓慢漂移，先怀疑 1/ε 的傅里叶因子规则没有做对，而不是继续增加 N。

5. 一维光栅 RCWA 常见问题排查：五个高发翻车现场与解决路径

5.1 TE 与 TM 算出来几乎一样：截断阶数太小

现象：把光源分别设为 TE、TM 偏振，算出的衍射效率曲线上只差几个百分点，甚至重叠。原因：截断阶数 N 太小，比如 N=1，整个模型只有 3 个空间谐波，无法描述光栅槽内的倏逝场分布，特征系统退化到类似等效介质的结果，偏振差异自然消失。解决：把 N 逐步加到 31、51，观察效率是否分裂出偏振差异；如果仍然接近，再检查是不是把目标波长上的高阶传播级漏掉了，比如某些级次的 k_z 刚好落在传播阈值附近。

5.2 TM 偏振收敛很慢：只展开了 ε，没展开 1/ε

现象：TE 在 N=31 稳定，TM 加到 N=99 还在漂移，且效率曲线出现高频抖动。原因：TM 偏振需要同时满足 E_x 和 H_y 的连续性，而直接使用 [ε] 的 Toeplitz 矩阵会在介电常数突变处产生吉布斯振荡，导致切向场连续性条件按傅里叶级数永远无法闭合。解决：改用 1/ε 的傅里叶系数参与构造 TM 层的特征系统，再用类似“傅里叶因子规则”的方案处理非对角项。这个修正之后，TM 的收敛速度通常会接近 TE，抖动也会消失。

5.3 深光栅出现 NaN/Inf：T 矩阵的指数放大

现象：光栅深度超过 2λ 后，某个波长点直接算出 NaN，附近点的效率出现无规律跳变。原因：代码里用的是 T 矩阵逐层传递，衰减模式的 exp(+k₀qᵢz) 在深槽里指数爆炸，双精度浮点根本存不下这个中间量。解决：把层与层之间的匹配全部改成 S 矩阵或散射矩阵形式；如果已经是 S 矩阵，检查 star_product 中的矩阵左除是否用了显式求逆，换成带求解器的左除可缓解部分条件数问题。

5.4 三角轮廓随分层数跳变：阶梯化策略不对

现象：三角形或正弦光栅，分层数从 32 改成 64，效率参考点从 0.41 跳到 0.44，再改成 128 又变回 0.42。原因：分层太粗时，三角顶点附近被人为切成了平顶，等于把原本的连续轮廓变成了一段段波导；每段之间的相位差随层数变化，恰好落在干涉敏感区就会放大为效率跳变。解决：让每层厚度小于 λ/20，同时轮廓斜率大的区域额外加密；不要用等厚度分层，要按弧长或曲率自适应切分。若光栅本身有圆角，应在圆角半径区域至少保证 3 层以上。

5.5 能量守恒误差忽正忽负：介质折射率与归一化因子不匹配

现象：无吸收光栅算完效率，反射加透射之和在 0.98 到 1.03 之间来回摆动，且随波长无规律变化。原因：透射区用基底折射率算 k_z，但功率归一化因子却用了入射区折射率；或者把倏逝级的 |T_m|² 也加进能量和。解决：把反射效率和透射效率分开归一化，反射用入射区的 k_z,0，透射用透射区的 k_z,0；另外，统计能量前先滤掉 Re(k_z) ≤ 0 的级次，只保留真正的传播级。

6. 设计取舍与进阶：复杂度、泛化与不完全靠 RCWA 的场景

RCWA 最大的代价是矩阵特征分解的 O(N³) 复杂度。N=21 时还能在笔记本电脑上轻松跑波长扫描，N=61 时每个波长点要做 123×123 复数矩阵特征分解，几百个波长点扫描下来就非常痛苦。一个工程习惯是：先用小 N 扫一遍粗糙曲线，找到共振峰或效率极值所在的波长区间，再在这些关键点用大 N 精算。如果对精度要求没那么苛刻，甚至可以只在波长的几个极值点做收敛性验证，其余点用 N=21 的结果。

从一维推广到二维光栅时，问题会突然大很多。二维方形柱阵列在 y 方向也要展开傅里叶级数，矩阵维度从 2N+1 变成 (2N_x+1)(2N_y+1)，内存按平方增长，计算量增长得更快。如果入射面不落在主截面，还会出现锥形衍射，波矢不再只在面内守恒，公式里会多出一整套 y 方向的边界条件。这也是为什么很多一维 RCWA 代码能跑得很好，一旦改二维就频繁崩溃——不一定是物理模型错了，而是矩阵规模触发了数值稳定性问题。

RCWA 并不适合所有光栅。当光栅周期大于 10λ 且结构深宽比很大时，传播级数量多得吓人，傅里叶展开逼近真实轮廓变得越来越吃力；此时 FDTD 或边界元方法可能更合适。相反，如果光栅周期接近波长、几何轮廓可以分层描述，RCWA 的效率和精度组合几乎是最优的。

从那以后，我每次做一维光栅仿真，都强制先设 N=41，算一次能量守恒，再扫波长；深槽结构一律走 S 矩阵，三角轮廓必做自适应分层。这个习惯让很多“玄学”变成了有明确误差界的数据。希望帮到你。这套完整实现，正好在资源包里整理成了带注释的脚本和算例，值得拿源码对照着跑一遍。

本文还有配套的精品资源，点击获取