RCWA与PWE：周期结构光学计算的两种核心方法解析

简介一套面向光学与电磁场研究者的MATLAB实现围绕严格耦合波分析RCWA与平面波展开法两大数值方法展开可用于周期性介质衍射、光栅设计、光子晶体带隙分析等场景。压缩包共24个文件以13个m脚本为核心涵盖setupGrid、convmat、calcLayer、calcReflectionSide等关键函数配合8个dat数据文件和README说明便于直接运行与二次开发整体仅17KB。已有82人学习浏览适合微纳光学、计算电磁学方向的学生与工程师快速上手。代码按RCWA与TransferMatrixMethod等模块组织既包含二维网格设置、PML边界处理也给出器件反射/透射计算流程能够帮助读者理解傅里叶展开、耦合波方程组求解及平面波叠加的实现细节为新型光栅、滤波器和波导结构设计提供可复用的仿真基础。 做周期性光学结构计算这几年我慢慢发现一个规律很多人一上来就开FDTD仿真软件网格一剖光子晶体一画波长一扫一个晚上就过去了结果只是为了验证一个本可以用几行MATLAB代码算出来的趋势。严格耦合波分析RCWA和平面波展开法PWE就是这样一对被低估的兄弟它们都建立在平面波展开的数学框架上但解决的问题完全不同RCWA给你算光栅的各级衍射效率PWE给你算光子晶体的能带结构。今年我为了一个超表面设计项目把这两套代码完整重写了一遍踩了不少坑也把实现思路理顺了这篇就按我实际写代码的顺序来聊聊。1. 先搞明白RCWA和PWE分别回答什么物理问题我在知乎上看到不少人把RCWA和PWE混为一谈实际上这两套方法虽然共享平面波展开的数学工具但面对的物理问题有天壤之别。RCWA回答的是给定一个入射光这个周期结构把光衍射到哪些方向、每个方向分到多少能量它是一个散射问题、一个输入输出问题。PWE回答的则是一个无限大的周期介质结构电磁波可以在哪些频率上存在它是一个本征值问题、一个模式分析问题。打个比方。RCWA像你拿手电筒照一块光栅问光被反射到哪些方向、哪些方向变亮了PWE像你敲一块周期性打孔的板子问这块板子本身有哪些固有频率。二者都属于周期结构的电磁分析但在实际项目里往往是配合使用的先用PWE确定结构的带隙位置和模式分布判断哪些波长能通过再用RCWA精细计算具体光栅的衍射效率。MATLAB里这两套方法的核心代码都能控制在几十行比商用电磁仿真软件灵活得多也更容易嵌入参数优化循环。从物理机制上说两种方法的共同出发点都是Floquet定理在周期介质中电磁场可以展开成一系列平面波的叠加叠加的权重由材料分布和入射条件决定。RCWA把这个展开应用到一个有限厚度的分层结构中逐层求解模式再匹配边界PWE则把这个展开代入无界周期介质的麦克斯韦方程直接把问题变成一个矩阵特征值方程。理解了这一点后面所有代码都只是同一思想在不同边界条件下的具体实现。方法对比速查表维度RCWAPWE问题类型散射/衍射本征模式/能带适用结构有限厚度周期光栅、超表面无限周期光子晶体核心输入入射角、波长、结构分层晶格类型、材料折射率、波矢路径核心输出各级反射/透射衍射效率频率-波矢色散关系附加能力可处理损耗、多层复杂剖面可提取带隙、模式分布典型局限依赖层状近似难以处理有限尺寸和缺陷2. RCWA实现从介电常数傅里叶展开到衍射效率输出RCWA写代码的思路可以拆成四步设置结构和入射参数、计算衍射级次的横向波矢、在光栅层内做特征分解、层间边界匹配并提取效率。下面按这个顺序拆开讲。2.1 结构分层与衍射级次波矢RCWA处理有限厚度的周期结构时首先沿光传播方向设为z把结构切成若干层。最简的一维光栅只需要一层光栅介质加上覆盖层和衬底层。如果结构剖面是倾斜的或者多层堆叠的就沿z方向多切几层切片厚度建议不超过结构最小特征尺寸的十分之一。以TE偏振、正入射的一维方波光栅为例我用这一段代码做参数初始化% 一维光栅RCWA参数设置 lambda 0.633; % 入射波长 period 0.4; % 光栅周期 ff 0.5; % 占空比 d 0.3; % 光栅深度 n_cover 1.0; % 覆盖层折射率(空气) n_grat 2.0; % 光栅材料折射率 n_sub 1.5; % 衬底折射率 N 10; % 截断级次取 m -10 .. 10 共21阶 m (-N:N); % 正入射波矢 theta 0; phi 0; k0 2*pi/lambda; kx_inc n_cover * k0 * sin(theta) * cos(phi); % Floquet条件各级衍射波矢 kx kx_inc - m * (2*pi/period);关于kx kx_inc - m * (2*pi/period)这行初学者容易搞混符号其实符号只取决于Floquet展开的约定方向不同文献可能差一个负号。判断对错的方法是看后面计算的衍射角是否符合物理直觉对正入射m0对应偏向一侧的衍射级m0对应另一侧。另外注意当某个衍射级次的横向波矢绝对值超过k0时这个级次在z方向的波矢分量会变成纯虚数它成为倏逝波不携带远场能量只影响近场耦合。2.2 介电常数的傅里叶系数与Toeplitz矩阵RCWA的关键操作是把光栅层的介电常数展开成傅里叶级数。对占空比为ff的方波光栅介电常数ε(x)在n_cover²和n_grat²之间切换傅里叶系数可以解析写出% 方波光栅介电常数的傅里叶系数 eps_avg ff * n_grat^2 (1-ff) * n_cover^2; eps_diff n_grat^2 - n_cover^2; eps_coeff zeros(2*N1, 1); for idx 1:length(m) mm m(idx); if mm 0 eps_coeff(idx) eps_avg; else eps_coeff(idx) eps_diff * sin(pi*mm*ff) / (pi*mm); end end % 由傅里叶系数构建Toeplitz卷积矩阵 E toeplitz(eps_coeff(N1:end), fliplr(eps_coeff(1:N1)));这个Toeplitz矩阵是空间域卷积在谐波域的表现形式。它之所以长成这样是因为两个谐波级次之间的耦合只取决于它们的级次差也就是m-m。手写这个矩阵时最常犯的错误是取错系数排列方向我的经验是先用小N验证一下矩阵是否满足厄米性质当介质无损耗时E矩阵必须是对称Toeplitz不满足就翻转移位。2.3 光栅层内的特征值与模式匹配TE偏振下光栅层内的电磁场问题可以化成一个二阶常微分方程在归一化坐标下变成求矩阵Kx^2 - E的特征值。这段代码做了核心计算% 归一化横向波矢对角阵 Kx diag(kx) / k0; % TE偏振的特征矩阵 A Kx^2 - E; % 特征分解得到本征传播常数 [V, D] eig(A); gamma sqrt(diag(D)); % 归一化传播常数复数包含损耗/倏逝信息这里gamma的物理意义很直接实部对应传播模式的衰减/相位变化虚部对应倏逝模式的指数衰减。对无损耗介质A是对称实矩阵本征值应为实数或成对共轭如果算出来明显不对称大概率是前面的Toeplitz矩阵构错了。拿到本征模之后还需要做层间边界匹配。这一步的理论是电场和磁场的切向分量在界面上连续实际操作中建议用S矩阵算法而不是T矩阵算法。T矩阵在层厚较大或者材料损耗较强时会包含exp(gamma*d)这类指数增长项导致数值溢出S矩阵算法把整个多层结构从顶部到底部逐层连接避免了指数插值放大稳定性好很多。新手如果发现反射率超过1或者出现NaN十有八九是用了T矩阵连接。边界匹配的代码量较大核心思路是把每个层内的前向和后向模式振幅作为未知量组装一个线性方程组最终从总S矩阵中提取反射和透射衍射效率。效率的计算公式是所有衍射级次的振幅平方乘以相应方向余弦的比值。2.4 TM偏振和Li规则如果你只需要算TE偏振上面代码基本够用。但要算TM偏振就必须面对RCWA历史上最有名的一个坑Li的傅里叶分解规则也叫Li规则。简单说TM偏振时方程中出现了(1/ε) * (某个场分量)这种乘积项。如果直接把1/ε的傅里叶系数展开再做标准卷积Field分量在介电常数突变界面上的切向边界条件无法自动满足导致结果严重偏离物理。Li在1996年证明了对这种第一类间断函数的乘积正确的做法是在谐波域使用倒规则——先把ε的傅里叶系数矩阵求逆而不是对1/ε直接做卷积。对应的代码差别很微妙% TM偏振错误做法是对1/eps的系数直接toeplitz p a hrefhttps://download.csdn.net/download/xinkai1688/91926403 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p