从麦克斯韦方程到散射矩阵：RCWA 严格耦合波分析详解

一束光照到周期光栅上，哪些方向会出现衍射光？每个方向分走多少能量？把光栅换成一层纳米柱，透射光又会获得怎样的相位？

光栅方程可以给出允许的传播方向，却不能独自回答各级衍射效率和相位。要回答这些问题，需要在实际材料分布中求解麦克斯韦方程。严格耦合波分析（Rigorous Coupled-Wave Analysis，RCWA）利用结构的横向周期性，把这个电磁边值问题转化成有限维的矩阵问题。

它的主线可以概括为：横向用傅里叶谐波表示，纵向用层内本征模传播，层间用电磁边界条件连接，再由散射振幅计算功率与相位。

图 1｜RCWA 的几何与端口。图中箭头用于说明方向和衍射级，不代表后文算例的定量角度或能量。本文六组图均由明确几何、公式或实际数值计算绘制，绘图与计算源码见文末。

本文面向学过电磁学、复数波动表示和线性代数的读者。第 1–5 节建立一维 TE 推导，第 6–9 节讨论稳定算法、TM 偏振和二维扩展，第 10–12 节给出实现与验证，第 13–15 节把这些工具连接到超构透镜设计。

1. RCWA 在求解什么问题？

1.1 “严格”指的是方程层面

RCWA 也常称为傅里叶模态法（Fourier Modal Method，FMM）。对于线性、时间不变的周期分层结构，它直接处理频域麦克斯韦方程，没有预先采用标量衍射、近轴传播或只保留两束耦合波的近似。早期经典工作已经把光栅衍射写成了可计算的矩阵边值问题。Moharam 与 Gaylord，1981

但“严格”不意味着电脑输出自动精确。实际计算仍然受到傅里叶截断、几何分层、材料数据、数值精度和边界建模的限制。更准确的说法是：在所选物理模型下，通过增加展开自由度和改善离散表示，可以检验解是否趋于稳定。

1.2 最适合的结构形态

典型模型包含：上方均匀半空间、若干周期图案层、下方均匀半空间。每层内部允许介电常数随横向位置变化，但沿厚度方向保持不变。

矩形条纹光栅可以直接用一层表示；圆柱纳米柱在固定高度范围内也可直接用二维周期图案表示。斜侧壁、锥形柱或曲面则需要沿厚度切片，每片用不同横截面近似。沿厚度分层和横向傅里叶截断是两种独立的误差来源。

RCWA 的标准模型在横向是无限周期延拓。把有限透镜放进大超胞，计算的仍然是周期重复的透镜阵列；如果希望用它逼近孤立器件，还要检查超胞尺寸和重复像之间的耦合。

2. 先把坐标、材料和相位约定说清楚

本文先讨论一维光栅：沿 $x$ 以周期 $\Lambda$ 重复，沿 $y$ 不变，层叠方向取 $+z$，入射光从上方向下传播。

采用时间因子 $e^{-i\omega t}$，于是无自由电荷和电流的频域麦克斯韦方程为

$$ \nabla\times\mathbf E=i\omega\mu_0\mathbf H, \qquad \nabla\times\mathbf H=-i\omega\varepsilon_0\varepsilon_r\mathbf E. $$

这里假设材料线性、局域、各向同性且非磁性，$\varepsilon_r$ 为相对介电常数。对采用本时间约定的普通无源材料，损耗对应 $\operatorname{Im}\varepsilon_r>0$。下面的标量推导还要求 $k_y=0$，也就是入射面为 $xz$ 平面。

记号	含义
$\lambda_0$、$k_0=2\pi/\lambda_0$	真空波长、真空波数
$n_i,n_t$	均匀入射区与出射区的折射率
$\theta_i$	入射波与 $+z$ 方向的夹角
$G=2\pi/\Lambda$	一维倒格矢基本量
$m=-M,\ldots,M$	保留的傅里叶级次，数量 $N=2M+1$
$Z_0=\sqrt{\mu_0/\varepsilon_0}$	真空波阻抗
TE	$\mathbf E$ 沿 $y$，非零分量为 $E_y,H_x,H_z$
TM	$\mathbf H$ 沿 $y$，非零分量为 $H_y,E_x,E_z$

注意三种不同对象：$m$ 是傅里叶谐波编号；层内本征模通常是多个谐波的叠加；只有在均匀外部区域，单个谐波才对应一个独立的平面波衍射级。

3. 周期性怎样变成矩阵？

3.1 Bloch 条件决定允许的横向波矢

平面波入射时，场满足准周期条件

$$ E_y(x+\Lambda,z)=e^{ik_{x,0}\Lambda}E_y(x,z), \qquad k_{x,0}=k_0n_i\sin\theta_i. $$

因此可以写成 Bloch 相位乘上周期函数，即

$$ E_y(x,z)=\sum_{m=-\infty}^{\infty}e_m(z)e^{ik_{x,m}x}, \qquad k_{x,m}=k_{x,0}+mG. $$

这并不是假定场在每个周期完全相同；斜入射时，相邻周期之间带有确定的相位差。

在折射率为 $n_j$ 的均匀、无损外部区域，

$$ k_{z,jm}=k_0q_{jm},\qquad q_{jm}=\sqrt{n_j^2-(k_{x,m}/k_0)^2}. $$

对向 $+z$ 的出射波，选取传播时 $q_{jm}>0$、倏逝时 $\operatorname{Im}q_{jm}>0$ 的分支。$q_{jm}$ 为正实数的级次可以把能量带到远处；纯正虚数对应随距离衰减的倏逝波。等号 $q_{jm}=0$ 是级次开启或关闭的截止点。

由横向动量关系可直接得到

$$ n_j\sin\theta_{jm}=n_i\sin\theta_i+m\frac{\lambda_0}{\Lambda}. $$

这个式子回答“往哪里走”；RCWA 接下来求解“每个方向的复振幅是多少”。

3.2 介电常数的乘法变成卷积

每层内的周期介电常数展开为

$$ \varepsilon_r(x)=\sum_{p=-\infty}^{\infty}\varepsilon_p e^{ipGx}, \qquad \varepsilon_p=\frac{1}{\Lambda}\int_{-\Lambda/2}^{\Lambda/2} \varepsilon_r(x)e^{-ipGx}\,dx. $$

空间中的乘积 $\varepsilon_r E_y$ 对应傅里叶系数之间的卷积：

$$ (\varepsilon_rE_y)_m=\sum_n\varepsilon_{m-n}e_n. $$

截断后，定义

$$ \mathbf e=(e_{-M},\ldots,e_M)^\mathsf T, \qquad \mathcal E_{mn}=\varepsilon_{m-n}, \qquad K_x=\operatorname{diag}\left(\frac{k_{x,m}}{k_0}\right). $$

$\mathcal E$ 是 Toeplitz 卷积矩阵，对角线相同位置差的元素相等。材料均匀时，只有零阶介电系数非零，谐波彼此解耦；材料有横向变化时，非对角元素把不同谐波连接起来。这就是“耦合波”的矩阵来源。

实现时有一个常见细节：场只保留 $-M$ 到 $M$，但索引差 $m-n$ 会覆盖 $-2M$ 到 $2M$。构造卷积矩阵所需的材料傅里叶系数范围，比场的级次范围更宽。

3.3 矩形条纹的解析系数

设一个周期内，居中的高折射率条纹宽度为 $F\Lambda$，占空比为 $F$，条纹和间隙介电常数分别为 $\varepsilon_a,\varepsilon_b$。则

$$ \varepsilon_p=\varepsilon_b\delta_{p0} +(\varepsilon_a-\varepsilon_b)F\,\operatorname{sinc}(pF), \qquad \operatorname{sinc}(u)=\frac{\sin(\pi u)}{\pi u}. $$

这里的 sinc 与 NumPy 的 np.sinc 定义一致。条纹中心若平移到 $x_c$，非均匀部分还要乘 $e^{-ipGx_c}$。只改变单元坐标原点就可能改变某些衍射振幅的相位，因此比较两套程序时必须对齐几何和相位参考。

图 2｜左：不连续函数的傅里叶重建在界面附近出现 Gibbs 振荡，增加阶数主要使振荡区域变窄。右：介电卷积矩阵的非对角元素体现谐波耦合。图示采用占空比 0.45、介电常数 4 和 1。

4. 从麦克斯韦方程推导 TE 层内本征问题

4.1 把标量波动方程投影到傅里叶基底

在当前 TE 条件下，由麦克斯韦方程消去磁场，得到

$$ \frac{\partial^2E_y}{\partial x^2} +\frac{\partial^2E_y}{\partial z^2} +k_0^2\varepsilon_r(x)E_y=0. $$

对 $x$ 求导等价于给第 $m$ 个谐波乘 $ik_{x,m}$，材料乘法则变成卷积。逐项比较相同谐波的系数：

$$ \frac{d^2e_m}{dz^2}-k_{x,m}^2e_m +k_0^2\sum_n\varepsilon_{m-n}e_n=0. $$

收集成矩阵形式：

$$ \frac{d^2\mathbf e}{dz^2} +k_0^2\underbrace{(\mathcal E-K_x^2)}_{\mathcal A}\mathbf e=0. $$

原本同时依赖 $x,z$ 的偏微分方程，变成了仅依赖 $z$ 的常系数耦合常微分方程。

4.2 对角化后，每个模态独立传播

求解

$$ \mathcal A W=WQ^2, \qquad Q=\operatorname{diag}(q_1,\ldots,q_N). $$

$W$ 的每一列是一个层内模态在傅里叶基底中的权重；$k_0q_\nu$ 是该模态的纵向传播常数。

若当前层顶部为 $z=0$，则形式解为

$$ \mathbf e(z)=W\left[e^{ik_0Qz}\mathbf a^+ +e^{-ik_0Qz}\mathbf a^-\right]. $$

这一步说明：层内不必再沿 $z$ 用许多小步做数值积分，传播可以通过矩阵指数的对角元素直接计算。RCWA 因而常被称作半解析方法。

4.3 磁场不是额外猜出来的

由旋度关系

$$ H_x=-\frac{1}{i\omega\mu_0}\frac{\partial E_y}{\partial z}. $$

定义与电场同单位的切向磁场变量 $\mathbf v=-Z_0\mathbf H_x$，则

$$ \mathbf v(z)=WQ\left[e^{ik_0Qz}\mathbf a^+ -e^{-ik_0Qz}\mathbf a^-\right]. $$

记 $V=WQ$，在同一个参考面上有

$$ \begin{pmatrix}\mathbf e\\\mathbf v\end{pmatrix} = \begin{pmatrix}W&W\\V&-V\end{pmatrix} \begin{pmatrix}\mathbf a^+\\\mathbf a^-\end{pmatrix}. $$

这个矩阵把“模态振幅”转成可以直接应用边界条件的“切向场”。负号来自反向传播波的电磁场关系，不能因为只关心强度就省略。

5. 层间边界条件如何求出反射和透射？

在没有自由表面电流的水平界面上，切向 $E_y$ 和 $H_x$ 连续。设界面左、右两侧的正反向模态幅度分别为 $\mathbf a_L,\mathbf b_L$ 和 $\mathbf a_R,\mathbf b_R$，则

$$ W_L(\mathbf a_L+\mathbf b_L)=W_R(\mathbf a_R+\mathbf b_R), $$

$$ V_L(\mathbf a_L-\mathbf b_L)=V_R(\mathbf a_R-\mathbf b_R). $$

把输入端口 $(\mathbf a_L,\mathbf b_R)$ 与输出端口 $(\mathbf b_L,\mathbf a_R)$ 分开，可以写成线性方程：

$$ \begin{pmatrix}W_L&-W_R\\V_L&V_R\end{pmatrix} \begin{pmatrix}\mathbf b_L\\\mathbf a_R\end{pmatrix} = \begin{pmatrix}-W_L&W_R\\V_L&V_R\end{pmatrix} \begin{pmatrix}\mathbf a_L\\\mathbf b_R\end{pmatrix}. $$

求解这个矩阵方程，就得到界面的散射矩阵。外部均匀区域可以直接使用平面波基底，即 $W=I$；TE 情况下 $V=\operatorname{diag}(q_m)$。

从上方以零级单位振幅入射、下方没有入射时，输入向量只在 $m=0$ 位置取 1。整叠结构的解会返回所有保留级次的 $r_m,t_m$，其中既包含传播级，也包含倏逝级。

即使最终只想要零级透射，也通常不能只保留零级谐波。 近场中的高阶和倏逝分量会参与边界匹配，再通过耦合影响零级结果。

6. 为什么需要散射矩阵，而不是一路乘传输矩阵？

6.1 倏逝模带来的指数放大

设某个模态 $q=i\beta$，$\beta>0$。穿过厚度 $d$ 后，两个形式上的传播因子为

$$ e^{ik_0qd}=e^{-k_0\beta d}, \qquad e^{-ik_0qd}=e^{+k_0\beta d}. $$

后者并不意味着无源结构凭空放大能量，而是同一个倏逝解被从不合适的端面反向外推。在多层矩阵连乘时，同时携带极大和极小的数字会使条件数恶化，甚至溢出。

稳定 RCWA 的关键之一，就是重新组织未知量。经典文献讨论过增强透射矩阵等稳定实现；本文使用适合逐层组合的散射矩阵方法。它们都是对数值稳定性的处理，不应把某一种实现技巧和 RCWA 本身等同起来。Moharam 等，1995

6.2 把两个方向的幅度分别参考到两端

对于顶部为 $0$、底部为 $d$ 的一层，更适合数值计算的场表达是

$$ \mathbf e(z)=W\left[e^{ik_0Qz}\mathbf a_L +e^{ik_0Q(d-z)}\mathbf b_R\right]. $$

这里的 $\mathbf a_L$ 参考到顶部，$\mathbf b_R$ 参考到底部。对普通无源介质选好模态分支后，这两个传播因子都沿相应方向衰减。

定义

$$ \begin{pmatrix}\mathbf b_L\\\mathbf a_R\end{pmatrix} = \underbrace{\begin{pmatrix}S_{11}&S_{12}\\S_{21}&S_{22}\end{pmatrix}}_S \begin{pmatrix}\mathbf a_L\\\mathbf b_R\end{pmatrix}. $$

没有内部界面的均匀图案层，其传播散射矩阵为

$$ S_{\mathrm{prop}}= \begin{pmatrix}0&P\\P&0\end{pmatrix}, \qquad P=e^{ik_0Qd}. $$

图 3｜左：散射矩阵把两端输入映射到两端输出。右：指数增长与衰减因子说明朴素传输矩阵为什么可能病态；曲线为解析函数，不是结构的增益谱。

6.3 用 Redheffer 星积连接两段结构

设结构 $A$ 在前，结构 $B$ 在后，内部相连端口使用同一个模态基底和参考面。消去内部往返波后得到合成矩阵 $C$。为避免矩阵乘法次序混乱，先定义

$$ D=I-B_{11}A_{22},\qquad X=D^{-1}B_{11}A_{21},\qquad Y=D^{-1}B_{12}. $$

那么

$$ \begin{aligned} C_{11}&=A_{11}+A_{12}X, & C_{12}&=A_{12}Y,\\ C_{21}&=B_{21}(A_{21}+A_{22}X), & C_{22}&=B_{22}+B_{21}A_{22}Y. \end{aligned} $$

这就是一组与本文端口约定一致的 Redheffer 星积公式。$D$ 描述两段结构之间的多次反射；它接近奇异时，可能对应真实的强共振，也可能暴露数值问题，需要结合收敛和条件数判断。

公式写成逆矩阵便于理解，代码中应使用 solve(D, rhs) 求解线性方程，避免显式计算逆矩阵。散射矩阵能显著改善厚层和倏逝模造成的不稳定，但不能消除所有截止、模态简并和病态共振问题。

7. TM 偏振为什么更容易收敛缓慢？

7.1 差别来自材料边界条件

对于竖直的条纹侧壁，TE 电场 $E_y$ 沿界面切向，因而连续。TM 的 $E_x$ 是界面法向分量，会在介电常数跳变处发生跳变；连续的是法向电位移 $D_x=\varepsilon_0\varepsilon_rE_x$。

无限傅里叶级数下，等价变形可以给出相同结果；但有限截断之后，先相乘再截断，与先截断再相乘，并不总能互换。尤其当两个不连续函数的乘积恰好连续时，直接卷积容易给出很慢的收敛。这是傅里叶因子分解规则要处理的问题。Li，1996

定义第二个卷积矩阵

$$ \mathcal B=\left[\!\left[ 1/\varepsilon_r\right]\!\right], \qquad \mathcal E=\left[\!\left[\varepsilon_r\right]\!\right]. $$

双括号表示“先对函数求傅里叶系数，再组装有限卷积矩阵”。一般有

$$ \boxed{\mathcal B\ne\mathcal E^{-1}.} $$

对侧壁的切向电场，可以按直接乘积规则使用 $\mathbf d_t/\varepsilon_0=\mathcal E\mathbf e_t$；对法向电场，更合适的截断关系是

$$ \mathbf e_n=\mathcal B\frac{\mathbf d_n}{\varepsilon_0}, \qquad \frac{\mathbf d_n}{\varepsilon_0}=\mathcal B^{-1}\mathbf e_n. $$

因此不能在所有位置统一选 $\mathcal E$，也不能把所有 $\mathcal E$ 都替换为 $\mathcal B^{-1}$。应先判断分量相对材料界面是法向还是切向。降低介电重建曲线的振铃，也不能代替正确的因子分解。

7.2 一维 TM 的可实现矩阵形式

令 $u=Z_0H_y$。TM 的连续方程可写为

$$ \frac{\partial}{\partial z}\left(\frac{1}{\varepsilon_r}\frac{\partial u}{\partial z}\right) +\frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_r}\frac{\partial u}{\partial x}\right) +k_0^2u=0. $$

对于层内沿 $z$ 不变的一维条纹，采用与侧壁边界条件相容的因子分解后得到

$$ \mathcal B\frac{d^2\mathbf u}{dz^2} +k_0^2\left(I-K_x\mathcal E^{-1}K_x\right)\mathbf u=0. $$

于是本征问题变为

$$ \mathcal A_{\mathrm{TM}}W=WQ^2, \qquad \mathcal A_{\mathrm{TM}}=\mathcal B^{-1} \left(I-K_x\mathcal E^{-1}K_x\right). $$

对应连续的第二个切向场变量取 $\mathbf v=\mathbf E_x$，由 $\partial_z u=ik_0\varepsilon_rE_x$ 得到

$$ V=\mathcal B WQ. $$

后面的界面匹配与散射矩阵级联形式保持不变，只是 $W,V$ 的构造不同。在均匀外部区域，TM 的导纳因子为 $q_m/\varepsilon_r$，与 TE 的 $q_m$ 区别明显。

一个有用的自检是把材料设为均匀：此时 $\mathcal E=\varepsilon_r I$、$\mathcal B=I/\varepsilon_r$，TE 和 TM 都恢复 $Q^2=\varepsilon_r I-K_x^2$。

8. 从复振幅到衍射效率与吸收

反射率和透射率是功率比，并不总是等于复振幅的模平方。

对本文限定的均匀、无损、非磁性外部介质，且仅有上方零级入射，统一定义

$$ y_{jm}=\begin{cases} q_{jm},&\mathrm{TE},\\ q_{jm}/\varepsilon_j,&\mathrm{TM}. \end{cases} $$

其中 TE 的幅度变量是 $E_y$，TM 的幅度变量是 $Z_0H_y$。若入射幅度为 1，则每级效率为

$$ R_m=\frac{\operatorname{Re}y_{im}}{\operatorname{Re}y_{i0}}|r_m|^2, \qquad T_m=\frac{\operatorname{Re}y_{tm}}{\operatorname{Re}y_{i0}}|t_m|^2. $$

求和得到

$$ R=\sum_{m\in\mathcal P_i}R_m,\qquad T=\sum_{m\in\mathcal P_t}T_m,\qquad A=1-R-T, $$

$\mathcal P_i,\mathcal P_t$ 是两侧的传播级集合。无损外部中的单个纯倏逝波没有向远处输运的净法向功率，但它的场振幅可以非零，也会影响结构内部耦合。

例如正入射 TE 的零级透射，在两侧折射率不同时有 $T_0=(n_t/n_i)|t_0|^2$。若忘记这个因子，就可能误以为无源系统的“透射率”超过 1。

更一般地，应从时间平均坡印廷矢量出发：

$$ P_z=\frac12\operatorname{Re}\int_{\mathrm{cell}} (\mathbf E\times\mathbf H^*)\cdot\hat{\mathbf z}\,dA. $$

外部介质有损、各向异性，或存在两侧相干入射时，功率定义需要重新检查，不能机械套用上面的独立级次公式。S⁴ 等成熟实现直接提供层内功率通量计算接口。S⁴ 功率接口文档

另外，$A$ 很小只能说明离散系统近似满足功率平衡，不能证明 $T_0$、相位或近场已经收敛。第 11 节会给出一个直接例子。

9. 二维周期结构：纳米柱需要完整矢量 RCWA

圆柱、方柱或旋转纳米鳍通常在 $x,y$ 两个方向都有周期。此时使用二维倒格矢

$$ \mathbf G_{mn}=m\mathbf b_1+n\mathbf b_2, \qquad \mathbf k_{\parallel,mn}=\mathbf k_{\parallel,0}+\mathbf G_{mn}, $$

并将每个场分量展开为

$$ F(x,y,z)=\sum_{m,n}F_{mn}(z) e^{i\mathbf k_{\parallel,mn}\cdot\mathbf r_\parallel}. $$

一般不能把两种偏振分别当成互不耦合的一维标量问题。消去纵向分量后，可用切向状态

$$ \boldsymbol\Psi= (\mathbf E_x,\mathbf E_y,Z_0\mathbf H_x,Z_0\mathbf H_y)^\mathsf T, \qquad \frac{d\boldsymbol\Psi}{dz}=ik_0\mathcal M\boldsymbol\Psi $$

构造一阶本征系统。若保留 $N_G$ 个二维倒格矢，状态维数为 $4N_G$，通常可进一步整理为涉及 $2N_G$ 阶矩阵的二阶问题。界面上同时匹配四个切向场分量。

更难的一步是因子分解：曲线边界的法向会随位置改变，不能再全局地把某个笛卡尔分量认定为法向。实际实现常借助法向量场或局部偏振基底组织材料矩阵。S⁴ 的文档列出了几种不同构造方式，并说明它们在适用性和精度方面的差别。S⁴ 傅里叶模态构造说明

因此，本文的一维代码适合帮助理解条纹光栅；不能把它的占空比参数改名为“柱直径”，就当作二维纳米柱求解器。

10. 一份能够逐步实现的算法流程

对于给定波长、入射角和偏振，可以按下面的顺序实现：

1. 确定谐波集合。 建立级次到数组下标的映射，并计算所有横向波矢。
2. 构造材料矩阵。 对每层求 $\varepsilon_r$ 和必要的 $1/\varepsilon_r$ 傅里叶系数，再按索引差组装矩阵。
3. 求层内本征模。 得到 $Q,W,V$，按时间约定与辐射条件选取模态分支。
4. 构造界面散射矩阵。 使用切向场连续条件求解线性方程。
5. 构造层传播矩阵并级联。 按入射区、各层、出射区的顺序计算星积。
6. 施加入射振幅。 读取反射与透射复系数，并按功率通量归一化。
7. 重算目标量并检查收敛。 增加谐波、加密几何切片，并比较振幅、相位或其他真正关心的量。

本文代码中最核心的层内步骤如下。E 表示 $\mathcal E$，B 表示 $\mathcal B$，K 表示 $K_x$：

if pol == "TE":
    A = E - K @ K
else:
    A = solve(B, I - K @ solve(E, K))

eigenvalues, W = eig(A)
q = outgoing_sqrt(eigenvalues)
V = W * q[None, :]
if pol == "TM":
    V = B @ V

完整脚本还包含界面求解、星积、功率提取、验证和绘图：下载 rcwa_demo.py。运行环境需要 Python、NumPy、SciPy；绘图另需 Matplotlib。

python rcwa_demo.py --output-dir results --figures-dir figures

也可以把它作为模块使用。长度统一用微米：

from rcwa_demo import rcwa

result = rcwa(
    wavelength=0.633,
    period=0.8,
    layers=[(0.35, 0.45, 2.0**2, 1.0)],
    M=35, n_in=1.0, n_out=1.45,
    theta_deg=0.0, pol="TE",
)
print(result["R"], result["T"], result["A"])

layers 中每个元组依次为厚度、占空比、条纹相对介电常数、间隙相对介电常数，按从上到下排列。复介电常数允许描述层内吸收。

这份教学实现限定为一维条纹、普通各向同性介质以及无损外部区域；它会拒绝精确的模态截止点，未处理二维图案、一般各向异性和近简并模态的高级数值问题。生产计算应使用经过充分验证的求解器；例如 S⁴ 明确采用了傅里叶模态展开与散射矩阵连接层结构。Liu 与 Fan，2012

11. 实际算例：总透射很高，零级却很低

11.1 参数与衍射级

以下是便于教学的无色散介质算例，不对应某种材料的实验光学常数：

参数	数值
真空波长	633 nm
周期	800 nm
条纹厚度	350 nm
占空比	0.45
条纹 / 间隙折射率	2.0 / 1.0
入射区 / 出射区折射率	1.0 / 1.45
入射方式	正入射，分别计算 TE 和 TM

根据光栅方程，反射和透射侧都只有 $m=-1,0,+1$ 能传播；其他保留级次是倏逝波。比如透射侧一级衍射满足 $|\sin\theta|=633/(800\times1.45)$，而二级已经超过 1。

保留 $M=35$、即 $N=71$ 个谐波，运行本文脚本得到：

量	TE	TM
总反射率 $R$	0.05748508	0.03652440
总透射率 $T$	0.94251492	0.96347560
零级透射率 $T_0$	0.00726361	0.01081522
$-1$ 级透射率	0.46762565	0.47633019
$+1$ 级透射率	0.46762565	0.47633019
与 $M=50$ 的总透射率差	$4.17\times10^{-6}$	$2.15\times10^{-6}$

由于几何和入射具有镜像对称性，两侧一级效率相同。这个例子也说明：“总透射接近 1”并不等于“光大部分沿零级透过”。 光栅把主要能量分配到了两个一级方向。

表中小数用于复算核对，不表示这些位数都已达到连续问题的真实精度。$M=50$ 是更高截断的比较基准，也不是严格误差上界。

11.2 守恒与收敛要同时看

图 4｜同一结构的真实计算结果。左：总透射率随谐波数变化。右：相对 $M=50$ 结果的差，与能量残差作比较。曲线不必单调收敛；接近机器精度的能量残差也不能替代截断收敛。

即使只保留 $M=0$，离散系统也可以很好地满足 $R+T=1$，但此时根本没有表示真实存在的 $\pm1$ 衍射级。它的 TE 零级透射率约为 0.9512，与足够多谐波下约 0.00726 的零级透射相差极大。

这比单看总透射更能揭示错误：一个“很守恒”的计算，也可能把能量分到了错误的通道。

11.3 代码做了哪些验证？

脚本包含以下数值检查，详细结果保存在 validation.json：

• 均匀薄膜解析解： TE/TM、3 个入射角、3 个厚度，共 18 组，与独立 Fresnel 薄膜公式比较复反射和复透射振幅，最大绝对差约 $1.75\times10^{-14}$。
• 裸界面极限： 去掉所有图案层，空气到折射率 1.45 介质的正入射反射率为 0.03373594。
• 分层一致性： 把一层分为两层相同材料、各一半厚度，反射复振幅差小于 $10^{-10}$。
• 镜像与功率检查： 正入射对称光栅的正负级次效率一致，无损算例满足 $|1-R-T|<10^{-9}$。
• 有损与厚层检查： 正虚部介电常数产生正吸收；将图案层加厚到 50 μm 后，计算仍得到有限且守恒的结果。

这些检查能排除许多符号、参考面和级联错误，但不是与另一套完整 RCWA/FEM/FDTD 求解器的独立交叉验证。需要用于研究结果时，仍应增加与成熟求解器在同一模型、同一归一化下的比较。

12. 怎样判断 RCWA 结果值得相信？

12.1 对真正关心的量做收敛

若任务是求某一级效率，就监测该级效率；若要建立超构透镜相位库，就同时监测复透射振幅与相位。两次结果的相位差可以写成

$$ \Delta\phi=\arg\left(t^{(N_2)}[t^{(N_1)}]^*\right), $$

这样可以避免直接相减跨越 $-\pi/\pi$ 时造成的假跳变。当 $|t|$ 接近零，相位对微小误差极其敏感，此时应优先比较复数本身和透射功率，而不是硬给出相位精度。

还应分别加密横向谐波、纵向切片以及用 FFT 构造材料系数时的几何采样网格。金属、狭缝、尖角、高折射率反差和界面近场通常更难收敛。远场效率已稳定，并不意味着界面附近的点值场已经稳定。

12.2 截止与共振附近需要额外采样

当某级 $q_m\to0$ 时，系统接近 Rayleigh 截止，某些模态表示和导纳矩阵会退化。应检查截止两侧，避免把恰好落在奇点上的数值失败解释成物理发散。

窄共振也需要足够细的波长网格。RCWA 是频域算法，每个频点单独求解；扫得很稀时可能直接跳过窄线宽特征。增加傅里叶阶数与加密频率采样解决的是不同问题。

12.3 计算量为什么会增长得很快？

对直接使用稠密线性代数的实现，本征分解和矩阵求解通常具有近似立方时间复杂度、平方存储复杂度。对一维结构，矩阵维度随 $N$ 增长；对二维矩形截断，倒格矢数量为

$$ N_G=(2M_x+1)(2M_y+1). $$

如果同时把两个方向的截断尺度约翻倍，$N_G$ 约变为原来的 4 倍，单次稠密操作的时间可能接近 64 倍。实际比例还受对称性、求解器实现和硬件影响。

实践中可以缓存相同层的模态，复用同一结构对不同输入偏振的散射矩阵；但改变波长、材料或入射横向波矢之后，通常需要更新本征问题。

计算任务	RCWA 的特点	建模时重点检查
周期光栅、薄膜与周期超表面	傅里叶基底直接适配横向周期性	截断、因子分解和截止点
多波长宽带响应	每个频点单独求解，材料色散容易逐点代入	扫频采样、窄共振
有限器件与任意三维曲面	通常需超胞或较多切片	周期重复影响、几何逼近
FDTD/FEM 交叉验证	可比较同一周期单元或有限器件模型	边界、材料、源、参考面和功率定义必须一致

方法选择取决于问题结构，而不是某个算法名称是否听起来更“严格”。

13. 用 RCWA 建立超构透镜的单元响应库

超构透镜由许多不同几何尺寸或朝向的亚波长单元构成。常见设计路线是：先模拟一个单元被无限重复排列时的响应，再把这个响应作为实际透镜局部的近似。这个步骤称为局部周期近似（Locally Periodic Approximation，LPA）。

RCWA 很适合计算这种周期单元的复透射系数。已有超构透镜研究就使用 RCWA 单元库提取相位场，再以衍射传播估算聚焦表现。单元库与透镜传播的研究实例

13.1 应保存复数响应，而不只是相位

对单元几何参数 $\mathbf p$，记录

$$ t_0(\mathbf p,\lambda,\theta,\mathrm{pol}) =|t_0|e^{i\phi}. $$

除了复透射系数，数据库还应包含周期、厚度、材料色散、入射角、偏振、参考面、谐波集合和收敛信息。必要时保存反射、吸收和其他衍射级。

如果结构具有偏振转换，响应应写成 Jones 矩阵，例如固定线偏振基底下

$$ \begin{pmatrix}E_{x,\mathrm{out}}\\E_{y,\mathrm{out}}\end{pmatrix} = \begin{pmatrix}t_{xx}&t_{xy}\\t_{yx}&t_{yy}\end{pmatrix} \begin{pmatrix}E_{x,\mathrm{in}}\\E_{y,\mathrm{in}}\end{pmatrix}. $$

旋转纳米鳍的几何相位设计，需要同时关注目标偏振通道的转换效率；仅有正确的相位关系，还不足以保证高聚焦效率。

13.2 周期“亚波长”要相对于两侧介质判断

对正入射、正方晶格、无损均匀两侧介质，如果希望两侧只有零级可传播，一个简单条件是

$$ \Lambda<\frac{\lambda_0}{\max(n_i,n_t)}. $$

只满足 $\Lambda<\lambda_0$，不能保证高折射率基底内没有高阶传播通道。斜入射和一般晶格应直接用 $|\mathbf k_\parallel+\mathbf G|<n_jk_0$ 检查所有候选级次。

13.3 一个真实的一维响应扫描

为了展示“参数—相位—效率”的关系，本文继续使用一维 TE 程序，另取周期 280 nm、厚度 600 nm、条纹折射率 2.4、间隙 1、基底 1.45，在 633 nm 正入射下扫描占空比 0.08–0.90，保留 $M=18$。

图 5｜由附带 RCWA 程序实际计算。展开相位的扫描范围约为 6.66 rad，零级透射率约在 0.697–0.995 之间。这里的幅相对应固定顶部与底部参考面，材料折射率为教学常数。该扫描展示单元库的组织方法，不是二维圆柱超构透镜的设计结果，也未逐个参数完成高阶收敛认证。

扫描结果见 unit_cell_scan.json。相位覆盖超过 $2\pi$ 是有用的设计条件，但还要检查相位覆盖区域是否同时具备足够高的透射、制造可行性和对尺寸误差的容忍度。

比较不同单元时，应统一入射与输出参考面。单独给某些单元增加一段均匀介质传播，会改变它们的相位，进而污染相位匹配。还要区分软件输出的模态系数、切向电场系数和功率归一化系数；它们不应未经转换就混用。

14. 从单元库到有限孔径超构透镜

14.1 推导目标聚焦相位

考虑平面波正入射，透镜位于 $z=0$，希望在均匀出射介质中聚焦到 $z=f$。半径 $r$ 处到焦点的路径为 $\sqrt{f^2+r^2}$。为了让不同位置的场到达焦点时同相，需要透镜提供补偿相位

$$ \phi_{\mathrm{target}}(r,\lambda_0) =-\frac{2\pi n_t}{\lambda_0} \left(\sqrt{f^2+r^2}-f\right)+\phi_0 \pmod{2\pi}. $$

负号与本文 $e^{-i\omega t}$、沿正方向传播使用 $e^{+ikz}$ 的约定一致。$\phi_0$ 是空间统一的相位偏置，不改变单波长聚焦位置，却可能帮助找到整体透射更高的单元组合。

14.2 相位匹配要包含振幅与制造约束

最直接的映射，是给每个位置选择相位最接近目标的单元。更实用的做法是最小化复响应距离，例如

$$ \mathbf p^*(r)=\arg\min_{\mathbf p\in\mathcal C} \left|t_0(\mathbf p)-a_{\mathrm{target}}(r) e^{i\phi_{\mathrm{target}}(r)}\right|^2, $$

其中 $\mathcal C$ 只保留满足最小线宽、间距、纵横比和透射阈值的候选；$a_{\mathrm{target}}$ 与库中的 $t_0$ 使用相同振幅归一化。也可以把圆周相位误差、功率透射损失和尺寸敏感性分开加权。

不要直接用未包裹的普通相位差作为代价。例如 $1^\circ$ 与 $359^\circ$ 的物理相位差只有 $2^\circ$，不是 $358^\circ$。

图 6｜左：空气中波长 633 nm、焦距 10 μm 的目标相位截面，横向范围为 ±5 μm；这是解析相位需求，不是已经优化完成的透镜。右：RCWA 单元库到有限器件验证的设计流程。

14.3 局部周期近似什么时候需要警惕？

周期单元仿真默认该单元的邻居与自己相同；真实透镜相邻单元常常不同。因此，LPA 忽略了环境变化带来的部分非局域耦合。

对上述目标相位，边缘附近相邻单元的连续相位变化量可以估算为

$$ |\Delta\phi|\approx\left|\frac{d\phi}{dr}\right|\Lambda =\frac{2\pi}{\lambda_0}\Lambda\,\mathrm{NA}_{\mathrm{local}}, \qquad \mathrm{NA}_{\mathrm{local}}=\frac{n_t r}{\sqrt{f^2+r^2}}. $$

它提示：高 NA、较大周期或更短波长会使局部相位变化更快。但这不是通用的误差上界，更不能给所有材料和单元设定一个固定的 NA 失效阈值。共振强度、场束缚、相邻几何差别和偏振都很重要。全波研究也展示了高 NA 设计中局部近似可能明显偏离实际器件响应。大面积超表面的全波计算实例

可以先选取中心、边缘、相位回绕附近和几何变化显著的位置做邻域验证；再用包含若干不同单元的超胞评估局部偏转效率。超胞 RCWA 仍有周期重复边界，需要检查超胞选择是否影响结论。关键设计最终应结合有限孔径全波仿真或实验验证。

14.4 单元透射率不是聚焦效率

在 LPA 下，可以先构造透镜出口的近似复场

$$ E_{\mathrm{out}}(x,y)\approx t_0[\mathbf p(x,y)]E_{\mathrm{in}}(x,y), $$

再采用角谱传播等方法估计焦场。高 NA 和显著偏振效应下，需要检查标量传播是否足够，必要时使用矢量传播或全波场。

应分别报告总透射功率、目标偏振通道效率、焦距、焦斑宽度、旁瓣和聚焦效率。比如定义

$$ \eta_{\mathrm{focus}}= \frac{\displaystyle\int_{\Omega_f}S_z(x,y,z_f)\,dA} {P_{\mathrm{inc,aperture}}}. $$

必须明确焦点积分区域 $\Omega_f$、焦平面位置，以及分母是孔径上的入射功率还是总透射功率。两种分母衡量不同含义，不能在比较结果时混用。

15. 再向前一步：多波长与消色差设计

单波长相位覆盖 $2\pi$，并不意味着不同波长能聚焦在同一位置。对固定焦距目标，在无色散出射介质的简化条件下，令

$$ \Delta L(r)=\sqrt{f^2+r^2}-f, \qquad \phi(r,\omega)=-\frac{n_t\omega}{c}\Delta L(r)+\phi_0(\omega). $$

按本文的时间约定，传递函数相位对角频率的导数对应群时延，所以

$$ \tau_g(r)=\frac{\partial\phi}{\partial\omega} =-\frac{n_t}{c}\Delta L(r)+\phi_0'(\omega), \qquad \mathrm{GDD}(r)=\frac{\partial^2\phi}{\partial\omega^2}. $$

相对中心，边缘需要的延迟较小，以补偿更长的后续传播路径；上式中的负相对值不表示器件必须具有负的绝对传播时间。空间统一的 $\phi_0'(\omega)$ 提供整体延迟自由度。

若出射介质有色散，第一项应改写为 $-k_t(\omega)\Delta L$，群时延中出现 $dk_t/d\omega=n_{g,t}/c$，二阶导数也不能忽略。真实单元库同样应输入随频率变化的复材料参数，不能用一个固定折射率代替整个波段。

因此，消色差设计需要在一个固定结构上，同时满足多个频率的复振幅、相位以及必要的色散要求。连续宽带目标还要求检查频率采样是否足够密，而不仅是少数离散颜色的焦点重合。已有消色差超构透镜工作明确利用相位、群时延和更高阶色散控制来拓展工作带宽。Chen 等，2018

计算群时延时，要先沿频率对相位合理展开，再对 $\omega$ 求导；不能把对波长的导数直接当作群时延。若先得到 $d\phi/d\lambda_0$，则

$$ \frac{d\phi}{d\omega}=-\frac{\lambda_0^2}{2\pi c}\frac{d\phi}{d\lambda_0}. $$

共振附近和透射接近零的位置，数值微分尤其敏感，应加密采样并观察复透射轨迹。

一个可执行的多波长设计流程是：建立同一组候选几何在所有目标频率的复响应库；联合选择满足各频率要求的单元；固定整片透镜的布局；最后在各频率上检查焦距、焦斑与效率。对每个波长分别重新选择一片不同布局，属于多次单波长设计，不能据此证明同一片透镜消色差。

RCWA 在这条流程中提供的是周期结构的电磁响应及其色散。将它用于超构透镜时，还需要把单元响应、局部近似、有限孔径传播和器件验证连接起来，才能从“每个单元有合适相位”走到“整片透镜达到目标性能”。

参考文献与可复算附件

1. M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” JOSA 71, 811–818 (1981). DOI。经典耦合波矩阵建模。
2. M. G. Moharam, E. B. Grann, D. A. Pommet and T. K. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” JOSA A 12, 1068–1076 (1995). DOI。二元光栅的稳定实现。
3. M. G. Moharam, D. A. Pommet, E. B. Grann and T. K. Gaylord, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach,” JOSA A 12, 1077–1086 (1995). DOI。倏逝模与数值稳定性。
4. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” JOSA A 13, 1870–1876 (1996). DOI。不连续函数的傅里叶因子分解规则。
5. V. Liu and S. Fan, “S⁴: A free electromagnetic solver for layered periodic structures,” Computer Physics Communications 183, 2233–2244 (2012). DOI；官方文档。RCWA/FMM 与散射矩阵的开源实现。
6. “Nanoscale precision brings experimental metalens efficiencies on par with theoretical promises,” Communications Physics (2024). 原文。RCWA 单元库、出口场和聚焦效率的衔接。
7. “Fast multi-source nanophotonic simulations using augmented partial factorization,” Nature Computational Science (2022). 原文。有限大面积超表面计算及局部周期近似的局限。
8. W. T. Chen et al., “A broadband achromatic metalens for focusing and imaging in the visible,” Nature Nanotechnology 13, 220–226 (2018). 原文。相位、群时延和群时延色散的联合控制。

本文推导和代码统一使用前述时间、坐标与端口约定；阅读不同文献时，应先转换这些约定，再比较矩阵中的符号。本文图示、条纹算例和一维单元扫描为原创教学材料，未复现某篇论文的器件性能。

• 计算与六组绘图源码 rcwa_demo.py
• 数值验证、级次功率与收敛数据 validation.json
• 一维周期单元复透射扫描 unit_cell_scan.json