简介针对周期性介质结构中的光波传播与衍射问题这份MATLAB实现包提供了严格耦合波分析RCWA与平面波展开法的完整代码框架适用于光学工程、光栅设计、光子晶体等方向的研究者与学生。包内共24个文件以13个m脚本为主体涵盖RCWA主程序、网格生成、PML吸收边界、反射/透射计算等模块另有8个dat数据文件、1个txt说明、1个asv备份及1个md文档便于快速理解代码结构与运行流程。压缩包仅17KB小巧精炼。目前已有85人学习下载是快速上手两种数值方法的实用参考。通过源码可直观掌握傅里叶级数展开、耦合波方程组构建及平面波叠加的实现细节并可直接改造用于自定义周期性结构分析显著缩短理论到仿真的转化成本。1. 严格耦合波分析和平面波展开法的MATLAB实现先从两条倒空间思路说起严格耦合波分析RCWA和平面波展开法PWE的MATLAB实现是我在处理周期性光栅和光子晶体时最先选择自己动手写的一对频域工具。它们不依赖商业电磁仿真包也不需要调FDTD网格只要会用矩阵分解就能让衍射效率和能带图从你自己的数据里长出来。RCWA 擅长回答“一个周期光栅把入射光分到每个衍射级几次”PWE 擅长回答“这族周期结构在哪些频率存在传播模式”。适合刚接触光子模拟的人建立坐标系也适合熟手在收敛性问题上快速拆解算法。很多人以为这两个方法已经很老但在超构表面、硅光栅耦合器和光子晶体波导的设计里标准做法仍然是它们。MATLAB 里实现两者的最大价值不是性能而是把截断阶数、介电常数傅里叶系数这些容易当成“黑匣子”的参数全部暴露在眼前。后面各章会从原理一直走到可运行的代码和参数边界你可以把模型从单层光栅一路扩展到二维光子晶体。2. 先把倒格子语言备好严格耦合波分析和平面波展开法共同的数学根基2.1 光栅方程和布洛赫条件两种方法都在倒空间里截断场周期性介质里的电磁场不能再写成单一平面波。RCWA 处理的问题是入射波打在周期为 d 的界面后散射场的横向波数被周期性调制产生一系列衍射级满足光栅方程 (k_{x,m}k_0 \sin\theta - 2\pi m/d)。这个关系决定了每一级的横向波数也决定了哪些级次能真实传播、哪些是倏逝波。PWE 处理的则是无限周期介质中的本征模式场写成布洛赫展开 (\mathbf{H}(\mathbf{r})\sum_{\mathbf{G}}\mathbf{H}_{\mathbf{G}} e^{i(\mathbf{k}\mathbf{G})\cdot\mathbf{r}})其中 (\mathbf{G}) 是倒格矢(\mathbf{k}) 被限制在第一布里渊区。两种方法的共同点在于都把连续空间中的电磁场问题转换成倒空间中的矩阵问题。矩阵的大小完全取决于你保留多少个衍射级或倒格矢这就是截断阶数。从代码框架看RCWA 和 PWE 的步骤高度相似先构建周期结构介电常数的傅里叶系数再用这些系数组装某个矩阵最后做特征值分解或求解线性系统。理解了这一层后面所有参数调整都有了方向。2.2 方法选型算衍射效率用RCWA算能带用PWE别用反维度RCWAPWE适用结构有限厚度的周期光栅、多层膜无限周期介质的光子晶体输入入射波长、入射角、偏振、光栅几何晶格类型、晶格常数、介质分布主要输出各级反射/透射衍射效率、近场分布能带结构、带隙、模式频率边界处理入射/出射半空间和层间匹配周期性边界无外部入射数值代价矩阵规模随谐波数线性增长矩阵规模随倒格矢数平方增长选型最大的错误是拿 PWE 去推衍射效率。PWE 给你的只是本征频率没有外部入射波的耦合强度算不出“百分之几的光透过了第 m 级”。反过来拿 RCWA 去算无限周期光子晶体的能带也存在问题你需要引入复频率扫描或者加周期边界代价明显高于 PWE。我一般的判断标准是问题里有没有对结构的“外部入射”有选 RCWA没有只是问模式在哪些频率存在选 PWE。2.3 截断阶数与收敛性决定你矩阵规模的核心参数RCWA 的一维截断阶数 (N) 决定矩阵宽度为 (2N1)。PWE 在二维方格子中取倒格矢截断半径 (N_p) 时基函数数是 ((2N_p1)^2)。它们的共同规律是截断越大矩阵越稠密计算越慢结果越逼近真解但到某个值之后再增大只是把小数点后第三位推稳而已。收敛性的通用判据是做一组截断扫描。比如 RCWA 中固定波长和角度比较 (N8,12,16,20) 时的零级反射效率变化小于 0.1% 就可以认为稳定。PWE 中固定高对称点看第一布里渊区边界频率随 (N_p) 的变化通常 (N_p6\sim10) 对大多数介质柱结构足够。注意高折射率对比和尖锐边缘会让傅里叶系数衰减变慢需要更高截断。TM 偏振比 TE 收敛慢这是由介电常数傅里叶展开的性质决定的后面避坑章会展开。3. 严格耦合波分析的MATLAB实现一维光栅的最小可运行框架3.1 定义光栅参数和入射条件波长、周期、占空比和偏振一个能立即跑起来的 RCWA 模型参数越少越好。下面这套定义对应一维矩形光栅只在 x 方向周期变化入射面是 x-z 平面层内 y 方向均匀。这也是严格耦合波分析最经典的教学场景。lambda 1.55e-6; % 入射自由空间波长单位 m d 0.8e-6; % 光栅周期单位 m f 0.5; % 占空比一个周期内高介电材料占比 eps_g 12.25; % 光栅材料的相对介电常数比如硅 eps_b 1.0; % 背景材料空气 theta 30 * pi / 180; % 入射角相对表面法线 h 0.4e-6; % 光栅层厚度单位 m N 12; % 保留 -N..N共 2*N1 个衍射级 m (-N:N).; % 级次索引列向量这里的关键参数是周期和波长的比值。当 (d/\lambda) 较大时可传播的衍射级数多矩阵即使不涨也需要更多级次来描述能量分配。占空比 (f) 直接影响介电常数傅里叶系数定义为高介电材料在一个周期内占的比例。如果你模拟的是槽结构把背景和光栅材料对调即可。厚度 (h) 在这个单层模型里决定层内相位积累深沟槽光栅的 RCWA 对层内数值稳定性更敏感后面会说明。(N) 的初始值取 12 是经验做法先跑通再扫描收敛。3.2 组装介电常数Toeplitz矩阵与状态矩阵ARCWA 的核心是把空间周期的介电常数 (\varepsilon(x)) 展开成傅里叶级数。对于矩形占空比结构这个级数有解析形式零阶项是平均介电常数高阶项正比于 (\mathrm{sinc}) 函数。deltaE eps_g - eps_b; eps_m eps_b deltaE * f * sinc(m * f); % 介电常数傅里叶谐波 E toeplitz(eps_m); % 卷积运算对应的 Toeplitz 矩阵 kx sin(theta) - m * lambda / d; % 归一化横向波数 Kx diag(kx); A_TE Kx^2 - E; % TE 偏振状态矩阵代码中sinc(m*f)使用的是 MATLAB 的归一化 sinc 函数(\mathrm{sinc}(x)\sin(\pi x)/(\pi x))。当 (m0) 时它等于 1所以eps_m的第一个元素就是 ( \varepsilon_b \Delta\varepsilon f)对应空间平均介电常数。toeplitz(eps_m)构造的矩阵第 (i,j) 个元素是eps_m(i-j)恰好把傅里叶级数在倒空间的卷积变成了矩阵乘法。这是整个 RCWA 最容易写错的地方很多翻车都是从傅里叶系数正负号或索引错位开始。横向波数kx直接来自光栅方程。注意这里用了 (\sin\theta - m\lambda/d) 的符号约定它决定了倒空间索引的方向无论用正号还是负号最终衍射效率相同但和后续边界匹配的相位定义必须保持一致。A_TE Kx^2 - E来自 TE 偏振的波动方程电场只有 y 分量对 z 的二阶导数由横向波数平方和介电常数共同决定。TM 偏振需要换成逆介电常数矩阵见避坑章。3.3 求解层内模式特征值分解与传播常数分支选择状态矩阵A_TE的特征值就是层内归一化传播常数的平方。求解特征值后还要处理一个物理分支问题每个特征值开平方得到两个符号必须选择沿正 z 方向传播或衰减的那一支。[U, D] eig(A_TE); q2 diag(D); q sqrt(q2); flip (real(q) 0) | (abs(real(q)) 1e-12 imag(q) 0); q(flip) -q(flip); q q(:); % 确保列向量这段代码看起来简单但它是最容易引入玄学错误的地方。sqrt对负实数默认返回正虚部对复数默认返回实部非负的主分支但不保证满足沿正 z 方向的辐射条件。我一般把传播模式要求为实部大于 0倏逝模式要求为虚部大于 0再统一翻转符号。特征向量U不需要翻转因为 (q) 和 (-q) 对应同一个特征向量列。特征向量矩阵U的每一列代表层内一个“耦合波模式”它的行对应不同衍射级。后面的边界匹配就是把这些模式按边界条件叠加而不是直接对每个级次单独求解。3.4 边界匹配与衍射效率用分块线性系统解出反射和透射系数单层光栅的边界条件可以写成 4 组未知量反射系数 (R)、透射系数 (T)、层内正向模式幅度 (c_)、反向模式幅度 (c_-)。把电场和磁场的切向连续性组装成一个大分块矩阵一次\求解就能得到所有系数。k0 2*pi/lambda; kIz sqrt(eps_b*k0^2 - (kx*k0).^2); % 入射介质 z 向波数 kIIz sqrt(eps_b*k0^2 - (kx*k0).^2); % 出射介质这里是悬浮膜对称结构 n 2*N1; I0 zeros(n,1); I0(N1) 1; % 零级入射振幅为 1 P exp(1i*k0*q*h); % 正向传播相位 Pm exp(-1i*k0*q*h); % 反向传播相位 M [eye(n), zeros(n), -U, -U; -diag(kIz), zeros(n), -k0*U*diag(q), k0*U*diag(q); zeros(n), -eye(n), U*diag(P), U*diag(Pm); zeros(n), -diag(kIIz), k0*U*diag(q)*diag(P), -k0*U*diag(q)*diag(Pm)]; rhs [-I0; -diag(kIz)*I0; zeros(2*n,1)]; sol M \ rhs; R sol(1:n); % 反射系数对应 m -N..N T sol(n1:2*n); % 透射系数组装逻辑来自四条关系式在 (z0) 面电场和磁场切向连续在 (zh) 面电场和磁场切向连续。代码第一行对应电场连续第二行对应磁场连续第三、四行是光栅层内部的场传播关系。矩阵中-k0*U*diag(q)前面的k0来自对 z 方向求导后归一化不能丢掉丢掉它会得到看起来很合理但能量不守恒的结果。衍射效率的计算需要按能流归一化。电场系数平方只是振幅比真实衍射效率要乘上对应级次在 z 方向的波数实部比kIz0 kIz(N1); ref_eff abs(R).^2 .* real(kIz) ./ real(kIz0); tra_eff abs(T).^2 .* real(kIIz) ./ real(kIz0); sum(ref_eff) sum(tra_eff) % 悬浮膜无吸收时应该等于 1这个归一化对反射和透射都成立因为能流的 z 分量正比于 (\mathrm{Re}(k_z)|E|^2)。上下介质不同时只要分别用对应的入射和出射波数公式仍然成立。如果总和不等于 1先别怀疑算法按第 5 章的分支选择去查代码。提示这个线性系统适合单层和较厚的均匀层。层数很多或者有深度光栅时直接拼接矩阵会出现数值病态工程上应该改用散射矩阵S matrix递推。4. 平面波展开法的MATLAB实现二维光子晶体能带的完整流程4.1 从晶格常数到倒格矢生成波矢基组平面波展开法要先把实空间晶格翻译成倒空间基组。以二维正方晶格为例原胞为边长 (a) 的正方形倒格矢也是正方形坐标间隔 (2\pi/a)。先生成波矢基组后面的矩阵组装都在这个基组里完成。a 1; % 晶格常数后面归一化会用到 R 0.2*a; % 圆孔半径 eps_b 12.0; % 背景介质介电常数 eps_a 1.0; % 孔内介质空气 Np 8; % 倒格矢截断半径以 2*pi/a 为单位 f pi*R^2/a^2; % 孔占空比 g 2*pi/a; [gx, gy] meshgrid((-Np:Np)*g, (-Np:Np)*g); G [gx(:).; gy(:).]; % 2 x nG 的倒格矢矩阵这里 (N_p8) 表示每维保留 (2N_p117) 个倒格矢总基函数数为 289。对于圆孔光子晶体这个数值已经能把前几条带算到工程可用精度。如果你做的是三角晶格或复杂原胞网格生成方式要换成晶格对应的倒格矢公式但后续矩阵组装逻辑完全一致。G的列顺序由meshgrid决定第一列是倒空间原点这正好对应后面的 Toeplitz 索引。4.2 计算结构因子和傅里叶系数圆孔的解析形式圆孔的傅里叶系数有解析解不需要用实空间网格做 FFT。关键在于把“孔内为空气、背景为高介电”的介电函数写成背景加修正项。修正项的傅里叶系数就是圆孔结构因子。gn sqrt(sum(G.^2, 1)); S zeros(size(gn)); S(gn 1e-12) f; % 零倒格矢处的结构因子等于填充率 nz gn 1e-12; S(nz) 2 * f * besselj(1, gn(nz)*R) ./ (gn(nz)*R); epsC eps_b * double(gn 1e-12) - (eps_b - eps_a) * S; invC (1/eps_b) * double(gn 1e-12) - (1/eps_b - 1/eps_a) * S;代码里besselj(1, x)是一阶贝塞尔函数这个公式来自圆盘在倒空间的填充因子(S(\mathbf{G})f\frac{2J_1(GR)}{GR})。当 (G0) 时结构因子等于填充率 (f)于是epsC(0)正好是空间平均介电常数 (\varepsilon_a f \varepsilon_b(1-f))。invC不是简单地对epsC取倒数而是对倒数介电函数 (1/\varepsilon(\mathbf{r})) 做同样的展开。这一步对 PWE 的 TM 模式是必须的跳过它就会出现本征频率偏移。4.3 组装本征方程并扫描高对称路径TM 极化电场沿 z 方向的本征方程是 [ \sum_{\mathbf{G}}\varepsilon^{-1}(\mathbf{G}-\mathbf{G})(\mathbf{k}\mathbf{G})\cdot(\mathbf{k}\mathbf{G})E_{\mathbf{G}} \frac{\omega^2}{c^2}E_{\mathbf{G}}. ] 写成矩阵形式时介电常数卷积矩阵和点积矩阵按元素相乘而不是普通矩阵乘法。Etoep toeplitz(epsC); invEtoep toeplitz(invC); function [freq, evecs] pwe_eig(kx, ky, G, invEtoep) Kgx kx G(1,:).; Kgy ky G(2,:).; Amat Kgx * Kgx. Kgy * Kgy.; % (kG_i) 点乘 (kG_j) H invEtoep .* Amat; H 0.5 * (H H); % 显式埃尔米特化压制数值噪声 [evecs, evals] eig(H); freq sqrt(max(real(diag(evals)), 0)); freq sort(freq); freq freq(freq 1e-8); % 丢弃数值零模 end函数中invEtoep .* Amat是关键Amat(i,j)是两个波矢点乘invEtoep(i,j)是倒数介电常数的傅里叶系数两者逐元素相乘才等价于对 (\mathbf{G}) 求和。显式埃尔米特化对中心对称结构没有物理副作用却能把eig返回的极小虚部清掉。频率计算取sqrt后要过滤接近 0 的模式因为平面波基底存在零频率平凡解它们不代表真实传播模式。4.4 画出归一化能带频率排序与实部提取能带图需要沿着布里渊区边界的高对称路径扫描波矢 (\mathbf{k})。正方晶格常用路径是 (\Gamma \to X \to M \to \Gamma)坐标以 (2\pi/a) 为单位分别是 ((0,0)\to(0.5,0)\to(0.5,0.5)\to(0,0))。npts 30; path1 [linspace(0, 0.5, npts); zeros(1, npts)]; path2 [0.5*ones(1, npts); linspace(0, 0.5, npts)]; path3 [linspace(0.5, 0, npts); linspace(0.5, 0, npts)]; klist [path1, path2, path3]; xaxis [linspace(0,1,npts), linspace(1,2,npts), linspace(2,3,npts)]; Nb 6; bands zeros(Nb, length(klist)); for i 1:length(klist) freq pwe_eig(klist(1,i), klist(2,i), G, invEtoep); bands(:,i) freq(1:Nb) * a / (2*pi); % 归一化频率 a/λ end plot(xaxis, bands, LineWidth, 1.2); set(gca, XTick, [0 1 2 3], XTickLabel, {Gamma,X,M,Gamma});归一化频率是 (\omega a/(2\pi c))也就是 (a/\lambda)。扫描时每个 (\mathbf{k}) 点都重新组装Amat但invEtoep只依赖结构可以在循环外算好。频率排序后直接取前 (N_b) 条带就能看到带隙位置。如果某些带在带隙边缘出现小突起优先怀疑倒格矢截断不够而不是物理现象。5. 严格耦合波分析和平面波展开法会遇到的5个坑从现象到修复5.1 现象衍射效率总和不等于 1单层悬浮膜、无吸收介质的情况下反射加透射效率应该严格等于 1。我第一次跑通时总效率只有 0.82看上去像是能量凭空消失实际原因在传播常数的分支选择。原因主要有两个一是对q取sqrt后没有统一翻转符号导致层内相位因子exp(1i*k0*q*h)指向了错误的衰减方向二是kIz和kIIz对倏逝级次取了错误的虚部符号让边界匹配方程把能量“吸走”了。解决方式是先做自检取一个没有光栅或占空比趋近 0/1 的极限这时多层反射率应该趋近无光栅薄膜的解析解。再对零吸收结构打印效率总和如果在 (1\pm10^{-6}) 内分支选择基本没问题。如果差值很大逐项检查q的实部和虚部符号特别要关注real(q)0的翻转逻辑。5.2 现象TM 偏振收敛很慢甚至振荡TE 偏振用介电常数的 Toeplitz 矩阵E收敛还算正常。TM 偏振如果照搬同样流程衍射效率会随截断阶数振荡很多级次的能量跳来跳去看起来像数值噪声。原因是 TM 极化的波动方程里出现的是 (1/\varepsilon) 算符而 (1/\varepsilon) 的傅里叶系数不等于介电常数傅里叶系数的倒数。直接对每阶系数取倒数等于改变了材料的空间分布当然不收敛。解决方式是在组装状态矩阵时使用逆介电常数的 Toeplitz 矩阵而不是1/E的元素级逆。也就是说先对 (1/\varepsilon(\mathbf{r})) 做傅里叶展开再构造 Toeplitz。这样 TM 的收敛速度会明显改善。更严格的还是遵循傅里叶变换的卷积规则不同极化方式处理方式不同TE 用 (\varepsilon) 的 ToeplitzTM 用 (\varepsilon^{-1}) 的 Toeplitz。5.3 现象PWE 能带图出现虚假平带用 PWE 算能带时本征值谱里经常出现一条完全平坦的带频率不随波矢变化。它不是物理上的缺陷态而是矩阵组装时引入了额外的零空间。原因通常是逆介电常数矩阵的零倒格矢项和结构因子处理不一致。比如圆孔结构里epsC(0)算错或者toeplitz(invC)的顺序没有和G的列顺序对齐导致卷积矩阵不对称平面波基底不能正确解耦。解决方式是先把倒空间原点确认清楚G的第一列必须是 ((0,0))epsC的第一个元素必须对应原点。然后再看本征值取实部时有没有被sqrt引入负值。过滤频率时阈值 (10^{-8}) 足够去掉伪零模但如果H没有埃尔米特化会留下幅值小于 (10^{-4}) 的虚部分量推荐代码里显式做0.5*(HH)。5.4 现象RCWA 薄层线性系统病态单层线性系统在光栅层很薄时能跑通但层数变多或层厚差异巨大时M \ rhs开始不稳定效率总和剧烈波动。我的经验是不要试图扩展这个线性系统到多层那不是小修小补而是换算法。原因是正向相位因子P和反向相位因子Pm在薄层里值非常接近边界匹配矩阵的列几近线性相关条件数爆炸。散射矩阵方法把每一层看成独立的散射单元通过递归合并避免直接求大线性系统数值稳定性好很多。常见做法是保留单层代码作为验证工具但正式的多层光栅全部改用 S 矩阵更新。如果只是自己研究单层光栅线性系统就够了如果要搭滤波片或层叠光栅尽早切换到 S 矩阵不要在薄层病态问题上花太多时间。5.5 现象二维 RCWA 结果和三维仿真对不上严格耦合波分析这个名字容易让人误以为它是全矢量三维方法但前面给出的实现只适用于“结构在一维方向周期、另一个横向方向均匀”的问题也就是入射面内衍射。如果光栅在 x 和 y 两个方向都有周期或者入射波倾斜到包含 y 分量衍射级次需要用二维倒格子索引代码里的横向波数会变成两个分量。原因很简单二维光栅的 TE/TM 偏振不再能用一个标量状态矩阵描述必须写出完整的 3×3 电磁场分量耦合。很多人在这个问题上翻车是因为拿一维代码去算二维光栅结果效率复数或者角度不对。解决方式是检查问题维度结构只在 x 方向周期、入射面在 x-z 平面才适用上面的标量公式否则需要升级到三维 RCWA矩阵规模随两个方向的截断阶数乘积增长参数调整原则仍然一致但代码复杂度不可同日而语。6. 从能跑通到可发表验证方法、一个加速技巧和我的自检习惯6.1 用光栅方程和薄膜极限做交叉验证我每次写完新的 RCWA 或 PWE 代码第一件事不是看目标曲线而是先跑两个廉价验证。RCWA 验证入射角和衍射级次位置固定波长和周期把零级透射峰对应的衍射角与解析光栅方程比对。PWE 验证极限结构把圆孔半径调到很小能带应该逼近均匀介质的色散线 (\omega c|\mathbf{k}|/\sqrt{\varepsilon_b})。这两个验证都不需要额外工具代码里打印几行数值就能完成。如果极限对不上不要继续往下调参数问题一定在傅里叶系数或矩阵组装中。我见过太多同学在复杂结构上调整半天最后发现是 sinc 函数多了一个 π。6.2 加速参数扫描复用 Toeplitz 矩阵做波长扫描或入射角扫描时RCWA 的结构矩阵E只依赖光栅几何不依赖入射波长和角度。PWE 的invEtoep只依赖介质分布不依赖波矢。把这两项从循环里提出来可以省掉大量重复构造时间。更深一层的技巧是在 PWE 中对同一批倒格矢先算好 (G_x,G_y) 的坐标矩阵每个 k 点只需要做一次实对称矩阵特征分解。对 MEMS 光栅或温调滤波器这类参数扫描场景这种复用能把单次扫描从几分钟压到几十秒。特征分解这一步无法避免但至少不要重复计算不变的部分。6.3 我的自检习惯和最后的建议我现在的习惯是在每个脚本开头写一个小节“自检开关”一个极限结构测试、一个功率守恒断言。RCWA 里断言abs(sum(ref_eff)sum(tra_eff)-1) 1e-6PWE 里断言最低带频率在均匀介质极限下与解析值误差小于 1%。有了这两个断言后面改动参数时数值上稍微不干净就会立刻报警。严格耦合波分析和平面波展开法其实是一对很好的“数值显微镜”它们让你看到截断阶数如何在倒空间里决定精度。先跑通单层光栅的 TE 模型再扩展 TM先算最基础的圆孔正方晶格能带再加超胞。每个扩展都在已有的矩阵框架上加一个维度而不是推翻重来。希望这篇实践笔记能帮你少走几个我走过的弯路也希望你在跑通第一条能带曲线时能感受到把矩阵变成物理图像的那种踏实感。本文还有配套的精品资源点击获取
