RCWA 算一维光栅衍射效率：从原理到代码的完整指南

简介这份资源围绕严格耦合波分析RCWA方法计算一维光栅衍射效率展开面向光学、光子学方向的学生与科研人员以及需要评估光栅性能的工程师。内容涵盖RCWA的基本原理、衍射效率的定义与物理意义并延伸到相位光栅衍射特性的研究适合具备一定波动光学与数值计算基础、希望快速上手光栅仿真的读者。资源包共2个文件包含1个m脚本与1个caj文献压缩包约4.84MB前者应为MATLAB实现的RCWA计算程序可用于设置光栅周期、填充因子、折射率等参数并求解各衍射级效率后者则提供相位光栅衍射特性的理论分析与研究参考。目前已有1945人学习下载。通过脚本与文献配合读者可理解RCWA相对FDTD等方法在计算效率上的优势掌握衍射效率随波长、入射角及结构参数变化的规律为光栅设计与光学器件优化提供可复用的计算思路与参考依据。1. 为什么 RCWA 算一维光栅比 FDTD 更值得先搞懂很多人第一次接触光栅衍射效率计算直觉是上 FDTD 或 FEM——时域有限差分、有限元通用、能算任意结构听起来更“稳”。但真到一维周期光栅这种场景你会发现一个反直觉的结论RCWA严格耦合波分析才是这类问题的“专用解”它不是近似方法而是把 Maxwell 方程在周期结构里做严格展开后的解析-数值混合解。换句话说FDTD 是拿网格去“逼近”场RCWA 是拿物理上天然存在的“衍射级次”去“展开”场后者在一维光栅里几乎不浪费自由度。这件事值得搞懂是因为它决定了你算衍射效率时到底在算什么。FDTD 里你调的是网格、PML、时间步RCWA 里你调的是谐波数衍射级次截断、层数切片、本征模求解。前者是数值逼近后者是物理展开。一维光栅的衍射效率——尤其是透射/反射各级次的能量分配——本质上是周期结构对入射波的模式耦合问题RCWA 的数学骨架天生就是干这个的。所以如果你要算的是光栅、超表面、光栅耦合器这类周期结构先把 RCWA 的底层机制吃透比盲目堆 FDTD 网格划算得多。2. RCWA 的数学骨架把 Maxwell 方程拆成“级次 本征模”2.1 从周期性和 Floquet 条件说起一维光栅在 x 方向周期为 Λy 方向均匀z 方向分层。入射平面波打在光栅上因为结构在 x 方向周期场在 x 方向必须满足Floquet-Bloch 条件任意 x 处的场等于 xΛ 处的场乘一个相位因子 exp(j·k_x0·Λ)。这直接推出一个结论——场可以展开成一组离散的“空间谐波”E(x, z) Σ_m E_m(z) · exp(j·(k_x0 m·K)·x)其中 K 2π/Λ 是光栅矢量m 是整数代表第 m 级衍射级次。k_x0 是入射波在 x 方向的波矢分量。这个展开不是人为选的是周期性边界条件逼出来的。衍射级次不是“算出来的”是结构周期性天然定义的——这是 RCWA 和 FDTD 最根本的分野。2.2 分层切片把连续结构变成“多层均匀光栅”一维光栅在 z 方向往往不是均匀的——比如矩形光栅有脊和槽梯形光栅有斜边。RCWA 的做法是沿 z 方向把光栅切成若干薄层每一层内认为结构在 z 方向不变只在 x 方向周期变化。层数足够多时阶梯近似收敛到真实轮廓。每一层里介电常数 ε(x) 也是周期的同样展开成傅里叶级数ε(x) Σ_p ε_p · exp(j·p·K·x)这里 ε_p 是介电常数的第 p 阶傅里叶系数。矩形光栅的 ε_p 有解析式梯形光栅则要数值积分。这一步是 RCWA 精度的第一个来源ε 的傅里叶展开阶数要和场的谐波数匹配否则会出现收敛问题后面避坑章会细说。2.3 本征模求解每层里场怎么传播把场的展开式和 ε 的展开式代进 Maxwell 方程利用各谐波的正交性可以把偏微分方程化成一组耦合的常微分方程。对 TE 偏振电场沿 y 方向耦合方程写成矩阵形式d²E_m/dz² -k0² · Σ_n A_mn · E_nA 是耦合矩阵元素由 ε 的傅里叶系数和波矢决定。这是一个矩阵本征值问题解出 A 的本征值和本征向量就得到每层里各个级次的传播常数和模式分布。层内场写成这些本征模的线性组合系数由边界条件确定。对 TM 偏振耦合矩阵形式不同涉及 ε 的逆的傅里叶展开这也是 RCWA 实现里最容易翻车的地方之一——TE 和 TM 的矩阵构造不对称不能简单套用。2.4 S 矩阵为什么不用传输矩阵直接连得到每层的本征模后要把所有层“串”起来。最朴素的做法是写传输矩阵T 矩阵逐层相乘。但 T 矩阵在层厚较大或级次较多时会出现数值指数溢出——因为里面含 exp(j·kz·d) 这种项kz 有正有负负虚部会指数增长。工程上通用的解法是S 矩阵散射矩阵把每层看成一个小散射体定义入射到该层的各级次波与出射波之间的线性关系然后用 Redheffer 星积把各层 S 矩阵级联。S 矩阵的构造保证数值稳定不会出现指数爆炸。这是 RCWA 能算厚光栅、多层堆叠的关键也是很多开源实现里代码量最大的部分。3. 从公式到代码一维 RCWA 的核心实现细节3.1 傅里叶系数矩阵怎么建以一维矩形光栅为例周期 Λ脊宽 w脊折射率 n1槽折射率 n2。介电常数 ε(x) 是方波傅里叶系数有解析式import numpy as np def eps_fourier_rect(n1, n2, w, Lambda, P): 计算矩形光栅介电常数的傅里叶系数 eps_p, p -P ... P n1: 脊折射率, n2: 槽折射率 w: 脊宽, Lambda: 周期, P: 最大谐波阶数 eps1, eps2 n1**2, n2**2 coeffs {} for p in range(-P, P 1): if p 0: # 第 0 阶 占空比加权平均 coeffs[p] eps1 * (w / Lambda) eps2 * (1 - w / Lambda) else: # 方波傅里叶系数注意 sinc 形式 coeffs[p] (eps1 - eps2) * (w / Lambda) * np.sinc(p * w / Lambda) \ * np.exp(-1j * p * np.pi * w / Lambda) return coeffs这段代码里有两个容易忽略的点。第一p0 阶是占空比加权平均不是简单平均因为 ε(x) 在脊和槽里各占不同宽度。第二非零阶的相位因子 exp(-j·p·π·w/Λ) 取决于你把坐标原点放在哪——原点选在脊中心还是脊边缘相位不同但最终衍射效率不变。这是常见的“结果对不上”来源两个人代码不同其实只是原点定义不同。对 TM 偏振还需要 ε 的逆的傅里叶系数。注意ε 的逆的傅里叶系数不等于 ε 傅里叶系数的逆必须重新对 1/ε(x) 做傅里叶展开。这是 TM 实现里第一个大坑。3.2 耦合矩阵的构造与截断拿到 ε_p 后构造耦合矩阵。对 TE 偏振矩阵元素是def build_A_TE(eps_coeffs, kx0, K, k0, P): 构造 TE 偏振的耦合矩阵 A eps_coeffs: dict, p - eps_p kx0: 入射波 x 方向波矢分量 K: 光栅矢量 2*pi/Lambda k0: 真空波矢 2*pi/lambda P: 最大谐波阶数 M 2 * P 1 # 总级次数 A np.zeros((M, M), dtypecomplex) for i in range(M): m i - P kxm kx0 m * K for j in range(M): n j - P # 对角项kxm^2 / k0^2 if i j: A[i, j] kxm**2 / k0**2 # 耦合项eps_{m-n} A[i, j] - eps_coeffs.get(m - n, 0) return A矩阵维度是 (2P1)×(2P1)。P 的选择直接决定精度和速度P 太小高阶衍射级次被截掉能量不守恒P 太大矩阵本征值求解变慢而且可能出现伪模。经验上对浅光栅 P10~20 够用对深光栅或金属光栅 P 要加到 30 以上。判断标准是能量守恒所有级次的透射反射效率之和应该等于 1无损耗时偏差超过 1e-3 就说明 P 不够。3.3 本征值求解与层内场表达对每层解 A 的本征值问题def solve_layer_modes(A, k0, d): 解 A 的本征值得到层内模式 A: 耦合矩阵 k0: 真空波矢 d: 层厚 返回: 本征值 q, 本征向量 W, 以及传播矩阵 eigvals, W np.linalg.eig(A) # 本征值 q 对应 kz^2 / k0^2 q np.sqrt(eigvals.astype(complex)) # 层内传播矩阵exp(j*k0*q*d) V W np.diag(np.exp(1j * k0 * q * d)) return q, W, V这里 q 是归一化传播常数k0·q 是实际 kz。注意 sqrt 的符号选择要保证透射层里场是衰减或向 z 传播的反射层里是衰减或向 -z 传播的。符号选错S 矩阵会给出物理上荒谬的结果比如反射效率大于 1。3.4 S 矩阵级联把各层“缝”起来每层解出本征模后构造该层的 S 矩阵然后用星积级联。核心逻辑def star_product(S1, S2): Redheffer 星积把两个 S 矩阵级联 S1: 上层 S 矩阵, S2: 下层 S 矩阵 每个 S 矩阵含四个块: S11, S12, S21, S22 I np.eye(S1[S11].shape[0]) # 中间矩阵避免直接求逆 D I - S2[S11] S1[S22] D_inv np.linalg.inv(D) S_new {} S_new[S11] S1[S11] S1[S12] D_inv S2[S11] S1[S21] S_new[S12] S1[S12] D_inv S2[S12] S_new[S21] S2[S21] S2[S22] D_inv S1[S21] S2[S11] S_new[S22] S2[S22] D_inv S1[S22] return S_new星积的物理含义是把两层之间的多次反射全部“折叠”进新的 S 矩阵。D 矩阵求逆是数值稳定性的关键如果 D 接近奇异说明两层之间存在共振耦合需要检查层厚或级次截断。级联完所有层后加上入射层和透射层的半无限边界条件就能解出各级次的反射和透射系数。3.5 衍射效率从系数到能量最后一步把各级次的场系数转成衍射效率。对反射第 m 级η_R,m Re(kz_R,m / kz_inc) · |R_m|²对透射第 m 级η_T,m Re(kz_T,m / kz_inc) · |T_m|²kz 的比值因子不能漏它来自能流密度Poynting 矢量的 z 分量。很多人只算 |R_m|² 就当成效率结果能量不守恒——这是新手最常见的错误之一。另外当 kz 是纯虚数时该级次是倏逝波不传播能量效率应为 0代码里要用 Re() 自动处理。4. 避坑与排查RCWA 算一维光栅最容易翻车的五个地方4.1 能量不守恒先查级次截断再查 kz 因子现象所有级次效率之和明显偏离 1比如 0.85 或 1.15。原因最常见的是谐波数 P 太小高阶衍射级次被截掉能量“漏”了。其次是效率公式里漏了 kz 比值因子或者 kz 的实部符号取错。解决先把 P 翻倍看能量是否收敛到 1。如果收敛说明是截断问题如果不收敛检查 kz 计算和效率公式。对金属光栅还要注意 ε 是复数kz 的 Re() 要正确取。4.2 TM 偏振结果异常ε 逆的傅里叶展开错了现象TE 偏振结果正常TM 偏振效率完全不对甚至出现负值。原因TM 偏振的耦合矩阵涉及 1/ε(x) 的傅里叶展开。很多人直接拿 ε_p 的倒数当 1/ε 的傅里叶系数这是错的。1/ε(x) 的傅里叶系数必须对 1/ε(x) 这个函数重新做傅里叶积分。解决单独写一个函数算 1/ε 的傅里叶系数。对矩形光栅1/ε(x) 也是方波只是脊和槽的值变成 1/n1² 和 1/n2²用同样的方波傅里叶公式即可。4.3 厚光栅数值溢出T 矩阵换成 S 矩阵现象光栅厚度增加后结果出现 NaN 或无穷大。原因用了传输矩阵逐层相乘里面含 exp(j·kz·d) 项kz 的虚部在厚层里指数增长浮点数溢出。解决改用 S 矩阵级联。S 矩阵的构造天然避免指数增长因为每层的 S 矩阵元素是有界的。如果已经用了 S 矩阵还溢出检查星积里的 D 矩阵求逆是否遇到奇异。4.4 收敛震荡阶梯近似层数不够现象增加谐波数 P 时效率不单调收敛而是来回震荡。原因对梯形或弧形光栅z 方向切片层数不够阶梯近似引入的误差和傅里叶截断误差耦合导致震荡。另外矩形光栅如果 ε 的傅里叶展开阶数和场的谐波数不匹配也会震荡。解决先固定 P增加 z 方向层数直到结果稳定再固定层数增加 P 直到能量守恒。两个方向分开收敛不要同时调。对矩形光栅ε 的傅里叶阶数建议至少取到 2P避免吉布斯现象。4.5 结果对不上文献先对齐坐标原点和符号约定现象代码逻辑没问题但衍射效率和文献差一个相位或符号。原因坐标原点定义不同脊中心 vs 脊边缘、时间约定不同exp(jωt) vs exp(-jωt)、傅里叶展开符号不同都会导致中间量不同但最终效率应该一致。如果效率不一致多半是 kz 符号或效率公式的 Re() 取错。解决先算一个已知解析解的简单案例——比如零深度光栅就是平面界面看反射率是否符合 Fresnel 公式。对齐这个基准后再算真实光栅。5. 设计取舍为什么 RCWA 在一维光栅里“刚刚好”RCWA 的定位很微妙它比解析方法比如严格耦合波之前的各种近似严格又比 FDTD/FEM 在周期结构里高效。这个“刚刚好”来自三个设计取舍。第一展开基的选择。RCWA 用 Floquet 模式展开这是周期结构里 Maxwell 方程的本征函数。FDTD 用网格上的局部基通用但自由度浪费在非周期方向。一维光栅里y 方向均匀z 方向分层只有 x 方向周期——RCWA 的自由度几乎全用在“刀刃”上。第二本征模求解 vs 时间步进。RCWA 解的是频域本征值问题一次求解给出所有级次的稳态响应。FDTD 要时间步进到稳态对高 Q 结构可能很慢。但 RCWA 的代价是每层都要解本征值层数多时矩阵大。所以它适合层数适中、周期方向自由度可控的场景。第三S 矩阵 vs T 矩阵。S 矩阵数值稳定但构造复杂T 矩阵简单但会溢出。这个取舍在实现里体现为代码量S 矩阵级联的代码通常是 T 矩阵的三四倍但它是 RCWA 能算厚光栅的前提。我自己的习惯是拿到一个一维光栅问题先用 RCWA 算一遍把 P 和层数收敛性跑清楚得到一个可信基准。如果结构有非周期缺陷或需要宽谱扫描再上 FDTD 做补充。从那以后我每次写 RCWA 代码都强制先跑零深度光栅的 Fresnel 基准再跑能量守恒检查最后才看衍射效率曲线——这三步走完基本不会翻车。希望帮到你。本文还有配套的精品资源点击获取