简介严格耦合波分析RCWA是分析周期性结构光学性质的常用方法这份资源围绕一维光栅衍射效率计算展开包含一个m格式脚本文件和一个caj格式研究文献。脚本可设置光栅周期、填充因子、材料折射率、入射波长、偏振与入射角等参数求解各衍射级次效率并输出结果文献详细展示了相位光栅衍射特性的分析流程可与脚本相互印证。资源共两个文件压缩包大小约4.84MB适合光学工程、光栅设计与超材料研究方向的科研人员以及具备一定波动光学基础的高年级学生。已有1939人浏览学习。通过该资源可快速掌握RCWA数值实现步骤理解不同衍射级次效率随波长、入射角及结构参数的变化规律为光谱仪、激光器和光通信器件中的光栅选型与性能优化提供理论依据和可修改的代码起点。1. RCWA计算一维光栅衍射效率先想清楚要解的到底是什么做衍射光学时一维光栅的衍射效率往往是第一个要定量算准的量。无论是设计分束光栅、抗反射表面浮雕结构还是给光谱仪选型拿到周期、占空比、膜厚和入射条件之后最直接的诉求就是知道每个衍射级次分走了多少能量。标量衍射理论在这个问题上只能给远场角度关系碰到亚波长周期或者折射率差比较大的结构时效率误差经常超过工程可接受范围。RCWA严格耦合波分析把光栅层里的电磁场按周期展开成一组平面波模态在傅里叶空间里数值求解麦克斯韦方程因此能同时覆盖偏振、入射角、层厚和材料吸收是一维光栅衍射效率计算的默认工具。这篇笔记会从物理模型讲到可运行代码再把收敛性、偏振处理、复折射率符号这些容易翻车的点单独拿出来谈。2. 从光栅方程到傅里叶模态RCWA一维光栅的物理模型2.1 衍射级次从哪来光栅方程的一次数值检查一维光栅的本质是在一个方向上引入周期性的折射率调制。光波入射到这个周期界面时切向动量会按空间频率发生平移每个平移量对应一个衍射级次。入射条件加上周期约束后衍射级次的切向波矢满足k_x,m k_0 n_i sinθ_i − m·(2π/Λ)其中 k_02π/λ 是真空波数n_i 是入射介质折射率θ_i 是入射角m 是整数级次Λ 是光栅周期。这个式子比常见的光栅方程多写了一个 m 的符号约定不同文献里 m 可能取正号也可能取负号实际影响只是级次编号方向相反效率数值不变。判断某个级次是否真实传播要看它在入射介质和出射介质里的纵向波矢是否为实数。纵向波矢为k_z,m sqrt( n²·k_0² − k_x,m² )如果根号内为正级次以平面波形式传播为负则成为倏逝波能量只在界面附近存在不进入远场。算衍射效率时倏逝级次虽然不直接贡献效率但在边界匹配中会参与场分布这也是 RCWA 必须保留足够高阶展开的原因之一。以空气入射、633nm 波长、1000nm 周期、正入射为例光栅方程给出的传播级次只有 m0、±1 三个。m2 时 k_x±4π/Λ对应纵向波矢已经是虚数。这个例子里周期略大于半波长但还远没到能够支撑更高阶传播的程度。实际设计时先查传播级次数再决定要保留多少展开阶数是避免后面收敛性翻车的第一步。2.2 周期介电常数怎么进麦克斯韦方程傅里叶系数与截断一维光栅的折射率分布是 x 方向的周期函数RCWA 的做法是把这个周期函数展开成傅里叶级数。设光栅脊的介电常数为 ε_r沟槽为 ε_g占空比为 f一个周期内脊宽为 f·Λ则介电常数的第 m 阶傅里叶系数为ε_m (ε_r − ε_g) / (j·2π·m) · ( 1 − exp(−j·2π·m·f) )m≠0零阶系数是占空比加权平均 ε_0 ε_r·f ε_g·(1−f)。注意这个表达式里包含了相位因子相位零点取在一个周期内脊的起点位置如果相位零点平移非零阶系数会整体乘一个相位项但最终效率结果不受影响。把介电常数展开后光栅层内的波动方程变成一组耦合方程耦合强度由 ε_{p−q} 这个 Toeplitz 矩阵决定。数值实现时必须截断无限级数把级次限制在 −M 到 M 之间共 2M1 项。截断阶数直接决定最终方程组的规模每层场的未知量、边界匹配矩阵的边长都随 N2M1 线性增长整体求解复杂度约与 N³ 成正比。M 取 5 时计算量很小M 取 30 时单次扫描已经能明显感到耗时因此实际工程里通常先做 M 收敛性扫描再固定一个相对小的值去做参数扫描。对于一维光栅如果占空比均匀、截面是简单矩形手写傅里叶系数并不难。一旦出现梯形侧壁或非矩形截面就需要把光栅层切成多个薄片每片用当地的占空比作为等效矩形处理。这种分层近似会让 RCWA 的计算量随层数线性增加但换来了能处理真实工艺形貌的能力。2.3 TE/TM偏振必须分家同一个光栅两种算法规格一维光栅在 y 方向是均匀的因此可以把偏振分为两类电场沿 y 方向的是 TE 偏振磁场沿 y 方向的是 TM 偏振。TE 情形下电场只有 E_y 一个分量波动方程是标量亥姆霍兹方程数值上最省事TM 情形下磁场只有 H_y 分量但电场在界面处存在法向不连续微分方程里会出现介电常数倒数 1/ε数值展开规则完全不同。这带来的直接后果是同一个一维光栅TE 和 TM 的效率曲线可能差异巨大。举个常见例子介质光栅占空比 0.5、周期 0.8 微米、厚度 0.6 微米正入射时 TE 可能在某个波长出现明显的 0 级透射凹陷TM 却不明显。原因是电场方向决定它主要感受的是介电常数分布还是其倒数分布对应的有效折射率变化量不同。RCWA 实现里 TE 和 TM 必须分别建立矩阵不能用一个偏振的代码直接套另一个。工程上常说的“偏振无关光栅”做设计时恰恰要同时跑两套 RCWA然后优化参数让两个偏振效率差足够小。只算一个偏振就下结论是新手最容易犯的建模错误。后续第 4 章会详细讲 TM 偏振在介电常数展开上的特殊处理这里先记住结论一维光栅 RCWA 从物理模型开始就要区分偏振。3. Python实现TE偏振RCWA从参数到衍射效率的最小可运行代码3.1 先摆好参数和光栅结构一份能直接跑的输入清单写代码之前先把一维光栅的结构参数整理成一张表。我习惯用纳米作为长度单位介电常数直接用折射率平方表示所有角度都转成弧度后再传入三角函数。参数含义示例值wl入射光真空波长633.0 nmperiod光栅周期1000.0 nmduty占空比脊宽/周期0.5h光栅层厚度300.0 nmtheta_deg入射角0 度eps_inc入射介质介电常数1.0eps_sub基底介质介电常数1.0eps_ridge脊材料介电常数2.25eps_groove沟槽介电常数1.0M正负保留级次数10这些参数对应一个悬浮在空气中的矩形介质光栅折射率 1.5 对空气。算完这个基准情况后把 eps_sub 改成玻璃的 2.25就能得到一个贴合实际的反射型一维光栅。入射角建议先从 0 度开始这样可以先避开掠射角附近的数值问题等结果稳定后再扫角度。输出衍射效率之前还要做一次传播级次预检确认当前周期和波长下哪几个级次是传播波。import numpy as np wl 633.0 period 1000.0 M 10 k0 2.0 * np.pi / wl K0 2.0 * np.pi / period m np.arange(-M, M 1) kx 0.0 - m * K0 # 正入射theta0 eps_inc 1.0 kz_inc np.sqrt(eps_inc * k0**2 - kx**2 0j) propagating np.abs(kz_inc.imag) 1e-12 print(传播级次 m:, m[propagating]) print(对应 kz:, kz_inc[propagating])这段代码用正入射条件直接生成全部截断级次的切向波矢然后通过纵向波矢虚部判断是否传播。虚部接近零表示实数传播波虚部绝对值较大则是倏逝波。参数说明里要留意0j是为了强制 kz 数组为复数类型避免后续运算里实部虚部被自动截断如果去掉numpy 的 sqrt 遇到负实数会报 nan这是新手最常见的小坑。3.2 介电常数傅里叶矩阵与光栅模特征值求解核心计算分两步先构造介电常数的 Toeplitz 卷积矩阵再求解光栅层内的特征值问题。TE 偏振下光栅层内每一阶场的耦合由 ε_{p−q} 矩阵描述横向算子写作A k₀²·[ε_conv] − diag(k_x²)这个 A 的特征值就是光栅层内允许的纵向传播常数平方。注意特征值可能为负对应的传播常数是纯虚数代表倏逝布洛赫模。def eps_hat(k, eps_ridge, eps_groove, duty): 方波光栅介电常数的第 k 阶傅里叶系数 if abs(k) 1e-12: return eps_ridge * duty eps_groove * (1.0 - duty) return (eps_ridge - eps_groove) / (1j * 2.0 * np.pi * k) * ( 1.0 - np.exp(-1j * 2.0 * np.pi * k * duty) ) N 2 * M 1 eps_conv np.zeros((N, N), dtypecomplex) idx np.arange(-M, M 1) for i in range(N): for j in range(N): eps_conv[i, j] eps_hat(idx[i] - idx[j], 2.25, 1.0, 0.5) kx_arr kx.reshape(N, 1) A k0**2 * eps_conv - kx_arr * kx_arr.T eigvals, W np.linalg.eig(A) gamma np.sqrt(eigvals) print(gamma 实部前 5 个:, np.real(gamma)[:5]) print(gamma 虚部前 5 个:, np.imag(gamma)[:5])上面这段以 3.1 节的参数代入eps_conv[i,j]存的是 ε_{i−j}也就是第 i 个衍射级与第 j 个衍射级之间的耦合系数。这个双重循环在小 M 下很快M 超过 20 后建议改用 scipy.linalg.toeplitz 一次性构造效果等价但速度快得多。np.linalg.eig返回的特征向量列顺序不固定每次运行可能有变化但特征值对应的物理意义不变边界匹配时只要使用同一组特征向量即可。gamma np.sqrt(eigvals)默认取主分支也就是实部非负的那个根对纯虚特征值来说正负分支代表向上和向下衰减的模这里暂时用主分支后续边界条件会通过方向判定校正。3.3 边界条件拼装与衍射效率提取把顶部入射区、光栅层、基底三个区域的场在光栅上下两个界面处匹配得到一组线性方程。未知量是反射系数 r、透射系数 t、光栅层正向模幅度 c 和反向模幅度 d方程组尺寸为 4N×4N。对一维 TE 光栅这组方程可以直接用 numpy.linalg.solve 求解。kz_inc_vec np.sqrt(eps_inc * k0**2 - kx.astype(complex)**2) kz_sub_vec np.sqrt(eps_sub * k0**2 - kx.astype(complex)**2) # 两块对角矩阵 diag_ki np.diag(kz_inc_vec) diag_ks np.diag(kz_sub_vec) diag_g np.diag(gamma) # 按 [r, c, d, t] 排列未知量 mat np.zeros((4*N, 4*N), dtypecomplex) rhs np.zeros(4*N, dtypecomplex) inc_vec np.zeros(N, dtypecomplex) inc_vec[M] 1.0 # m0 对应下标 M # z0 界面电场连续 mat[0:N, 0:N] np.eye(N) mat[0:N, N:2*N] -W mat[0:N, 2*N:3*N] -W rhs[0:N] -inc_vec # z0 界面磁场连续对 TE 等价为电场导数连续 mat[N:2*N, 0:N] -1j * diag_ki mat[N:2*N, N:2*N] 1j * W diag_g mat[N:2*N, 2*N:3*N] -1j * W diag_g rhs[N:2*N] -1j * kz_inc_vec * inc_vec # z-h 界面电场连续 P_minus np.diag(np.exp(-1j * gamma * h)) P_plus np.diag(np.exp(1j * gamma * h)) mat[2*N:3*N, N:2*N] W P_minus mat[2*N:3*N, 2*N:3*N] W P_plus mat[2*N:3*N, 3*N:4*N] -np.eye(N) # z-h 界面磁场连续 mat[3*N:4*N, N:2*N] -1j * W (diag_g P_minus) mat[3*N:4*N, 2*N:3*N] 1j * W (diag_g P_plus) mat[3*N:4*N, 3*N:4*N] 1j * diag_ks sol np.linalg.solve(mat, rhs) r_coef sol[0:N] t_coef sol[3*N:4*N] R_eff np.abs(r_coef)**2 * np.real(kz_inc_vec) / np.real(kz_inc_vec[M]) T_eff np.abs(t_coef)**2 * np.real(kz_sub_vec) / np.real(kz_inc_vec[M]) print(反射效率:, np.real(R_eff)) print(透射效率:, np.real(T_eff)) print(总效率:, np.sum(np.real(R_eff)) np.sum(np.real(T_eff)))边界条件拼装的核心逻辑是电场在界面处连续TE 偏振下磁场连续等价于电场对 z 的导数连续。矩阵每一行块对应一个物理条件列块对应一类未知量求解时把入射波放在 rhs 右侧。衍射效率计算不能直接用幅度平方还要乘上纵向波矢的比值因为不同级次的传播方向与能量流密度不同绝对值平方只是电场振幅比不是能量比。注意kz_inc_vec[M]是正入射 0 级波的纵向波矢做归一化基准。如果入射角不为零这里的归一化要换成入射级 m0 对应那个纵向波矢公式不变只是下标对应的切向波矢变了。代码跑完后最好立刻检查总效率是否接近 1误差在 1e-6 量级说明边界条件拼装正确误差在百分之几则说明符号或分支设置有误。4. RCWA避坑实录收敛、偏振、复折射率与掠射角的常见翻车4.1 截断阶数不够衍射效率随M剧烈振荡现象同一个光栅M 从 3 改到 5算出 0 级透射效率从 0.62 跳到 0.47M 再改成 8效率又变回 0.58。效率曲线不是随 M 缓慢变化而是在小范围 M 内剧烈振荡。原因保留的傅里叶阶数太少无法完整表示光栅层内的高次倏逝场。这些倏逝布洛赫模虽然在远场不传播能量但它们承担了界面处电磁场边界条件的匹配截断越少边界误差越大。尤其当周期接近波长、光栅厚度较厚时高次模的贡献不会快速衰减。解决先用小 M 跑一遍然后逐步增加 M比如 M5、8、12、16、20观察目标效率的变化量。变化量小于 1e-3 就可以认为收敛。工程经验是 M 至少要比最大传播级次多 5~8 个。对 633nm 入射、周期 1000nm 的情形传播级次是 0 和 ±1M 取 10 通常已经足够对于金属光栅或高折射率差M 往往要取到 20 以上。我现在的习惯是把 M 收敛性检查做成固定流程每个设计点最少跑一次不省这一步。4.2 TM偏振直接复用TE代码大折射率差时结果严重错误现象把 TE 代码原封不动用于 TM 偏振入射角约 30 度、折射率比 1.0 对 3.5 时算出的透射效率超过 1而且不满足能量守恒。原因TE 和 TM 的横向本征方程不同。TE 偏振的波动方程里只出现介电常数 ε 本身TM 偏振的波动方程在部分项里出现的是 1/ε因为电位移矢量的法向分量需要穿过介电常数突变界面连续。若直接对 ε 做傅里叶展开再求逆得到的等效介电常数在中高折射率差下会严重偏离物理实际这就是傅里叶因式分解问题文献里常称为 Li 规则。解决TM 偏振单独实现横向算子中用到的是介电常数倒数矩阵而不是介电常数矩阵本身。具体做法构造 ε_conv 之后再构造 (1/ε)_convTM 的耦合矩阵由这两个矩阵按规范组合而成组合方式与 TE 完全不同。对于折射率差较大的材料这是决定计算结果正确性的关键步骤。如果只是想快速核验 TM 结果最稳妥的方式是用成熟仿真软件跑一个对照算例而不是在自写 TE 代码上打补丁。4.3 复折射率虚部符号设反总效率竟然超过1现象计算金属铝光栅时反射效率加透射效率之和为 1.08明显超过能量守恒上限。原因复数折射率的虚部符号约定与时间因子绑定。若全程使用 e^{jωt} 约定吸收型材料的折射率应写为 n − jκ若用 e^{−jωt} 约定则应写为 n jκ。两套约定对应不同的光场相位演化方向复介电常数的虚部符号也必须随之取反。代码里如果从文献直接抄了折射率却没校准时间因子就会把吸收误算成增益导致能量越算越多。另一个常见来源是 gamma 的平方根分支吸收介质中的纵向波矢虚部必须保证场向传播方向衰减分支选反等于把吸收材料当成增益介质。解决在代码里加一道能量守恒自检计算 RT 之和。若是无吸收介质RT 应严格等于 1若是吸收材料RT 应小于 1 且等于 1 减去吸收率。一旦出现超过 1 的情况优先检查折射率虚部符号和 gamma 分支两个位置。我通常写一段测试代码先跑一块无光栅的均匀介质膜这时 RCWA 应退化为标准薄膜干涉结果以此验证符号约定没问题后再算光栅。4.4 掠射与临界角附近条件数过大奇异点伪装成衍射峰现象入射角从 80 度扫到 89 度某个级次的衍射效率在 86 度附近突然出现一个极窄的峰值峰值宽度不到 0.1 度但两侧效率数值正常。原因入射角接近掠射时某个衍射级次的纵向波矢接近零对应所谓瑞利异常。此时该级次从传播模式转为倏逝模式转换点附近边界匹配矩阵的条件数急剧增大数值解对舍入误差非常敏感容易在异常点附近出现假峰。另一个原因是边界匹配矩阵中同时存在传播模和倏逝模数量级相差悬殊导致线性方程组病态。解决先用条件数检查确认np.linalg.cond(mat) 超过 1e10 就要警惕。对接近掠射的入射角角度步长要加密并观察结果是否平滑过渡若在物理上不可能出现尖锐峰的波段看到窄峰优先怀疑数值伪峰而不是设计亮点。对金属光栅或高吸收结构适当增加 M 或改用归一化矩阵可以让病态问题缓解。实际工程扫描中我一般把入射角限制在 85 度以内超过这个范围的结果只做定性参考。5. 结果验证与工程化使用让RCWA真正进入设计闭环5.1 收敛性检查每个算例都花十几秒做M扫描拿到一组设计参数后不要直接信任 M10 的输出。先写一个循环让 M 从 6 跑到 20记录目标级次效率的变化。判断标准不是绝对收敛到某个固定值而是在你的精度需求下波动足够小。比如分束光栅要求效率比误差在 1% 以内那 M 扫描中波动小于 0.002 就可以接受。这一步十几秒就能跑完却能挡掉大量由于截断不足导致的返工。5.2 能量守恒与库符号自检每次修改代码或换材料库后先跑一个无光栅对照算例再跑一个金属光栅算例分别验证 RT1 和 RT1 到底成立不成立。这两个测试能快速抓住 4.3 节里说的符号错误也能发现边界条件矩阵拼装是否多了负号。我遇到过几次本来物理模型没问题、只是材料库折射率虚部带错符号的情况没有这道自检流程根本查不出来。5.3 批量扫参与优化注意点一维光栅 RCWA 单次计算在毫秒级适合直接做网格扫描。但要注意扫参时不同参数点对 M 的需求可能不同折射率差变大、占空比接近 0 或 1、光栅厚度增大都会让收敛变慢。固定 M 做全波段扫描前先在最小和最大折射率差两个极端点各做一次 M 收敛验证确认同一个 M 在整个扫描范围内都可靠。我现在的习惯是把 5.1 和 5.2 写成两个独立函数放进每次调用的公共模块里新项目启动时先跑一遍后面所有参数扫描都基于同一套验证过的代码。这个流程帮我挡掉多数自写 RCWA 的坑也希望帮到你。本文还有配套的精品资源点击获取