跟做超表面或者周期性光栅的人聊 COMSOL最常听到的一句话是透射谱算出来了但共振来源说不清楚。我自己的体会很深——谱线在某段波长掉下去一个深谷S21 从 0.9 掉到 0.15形态看着像教科书里的共振线型可真被人问一句“这到底是电偶极共振还是磁偶极共振跟前后的谷是不是同一个模式”光盯着谱线根本答不上来。这时候需要的是多极子展开。它把结构里感应出来的电流按电偶极、磁偶极、电四极、磁四极、环偶极这样一组空间对称性基函数展开每一项对应一种“辐射语言”透射谱上的谷、峰、不对称线型都能转化为某个或某几个多极矩的增强或互相干涉。COMSOL 里做周期性结构的光学仿真透射谱计算本身是成熟流程多极子分解则需要你在后处理中额外做一组体积分这也是很多人卡住的地方。这篇文章就把从建模到透射谱、再到多极子分解的完整链路讲一遍并给出可以直接抄作业的积分表达式和操作顺序。下文会以一个方形晶格的高折射率介质纳米圆柱阵列为例单元周期 P600 nm圆柱半径 R120 nm、高 H300 nm背景为空气材料折射率按 n≈3.5 的无损介质近似。这套参数在近红外波段会给出清晰可辨的透射谷足够我们分析共振归属。有基底、斜入射的情况我会在关键步骤处补充说明。1. 先搞清楚一件事透射谱里的谷到底是谁造出来的1.1 一个谷可以对应共振也可以对应两个模式的相消透射谱是平面波穿过周期性结构后零阶透射强度随波长的变化。听起来很直接但它实际上是整个单元里所有散射波相干叠加的结果。周期性结构处于亚波长周期时每个单元被同样的入射波激发但相邻单元之间还有多次散射最终远场的所有贡献互相干涉谱线上才出现峰和谷。一个谷可能来自一个共振模式被激发也可能来自两个模式的相消干涉。尤其是暗模式或类电磁感应透明的情况谱线上反而会出现一个透明的窗口窗口两侧才是峰。你很难从谱线形态直接判断是“哪种对称性的电流在振荡”——这就像一个乐队同时演奏时你听到整体音量变小却没法知道是哪把乐器在某个瞬间停止了。透射谱是结果模式藏在谱线后面。1.2 多极子展开为什么是回答这个问题的好工具多极子展开的本质是把一个体积内的感应电流分布写成一系列整体矩的叠加。这跟声音分解成基频和泛音是同一个思想。电偶极矩描述电荷重心做振荡磁偶极矩描述环形电流的整体旋转电四极矩描述两个反向偶极构成的形变环偶极矩则是一种首尾相接的电流环形成的闭合磁感线构型。每一项都有明确的辐射方向图和远场强度公式。对周期性结构来说多极子展开通常在单个颗粒内做体积分。每个粒子的感应电流以多极矩向外辐射单元之间的耦合效应则体现在多极矩随波长的变化里——因为外部入射场已经包含了周围单元散射场的贡献。所以多极矩并不是一个纯粹的“单颗粒属性”它是在阵列环境下提取出来的有效响应。1.3 展开到哪一阶才够用亚波长颗粒的尺寸远小于波长时高阶多极的远场贡献会迅速衰减通常算到电四极、磁四极和环偶极已经足够。如果结构接近波长或者里面有高折射率的强波导模式就需要考虑更高阶项。判断是否截断合理有一个很实用的验证把偶极项、四极项、环偶极项的辐射功率加起来和 COMSOL 算出的总散射/反射能量做对比如果差值很小说明展开截断够了。如果差很多不要急着加更高阶先排查原点和网格。这两个因素对高阶矩的影响往往比“项数不够”更致命。2. COMSOL 建周期性结构透射谱的关键边界条件优先于几何2.1 单胞几何和材料一个可以复现的起点几何不用做得太花哨。单元范围取 x∈[-300, 300] nmy∈[-300, 300] nm圆柱中心放在全局原点半径 120 nm高度沿 z 从 -150 nm 到 150 nm上下各留约 500 nm 空气层。把圆柱中心放到原点是一个很值得养成的习惯因为后面做多极子体积分时坐标原点就是多极展开的中心能少犯一类非常隐蔽的错误。材料用高折射率介质n≈3.5对应相对介电常数 εr12.25。先不要加损耗把流程跑通后再考虑色散和损耗。背景是空气background 折射率 n1。如果你想模拟有基底的情况最省事的办法是把空气背景换成等效折射率例如 n_bg1.45 的低折射率衬底但这会在透射谱里引入衬底本身的反射更严格的做法是把基底作为一层实体建模端口设置在端面处。2.2 周期性边界与端口“一个单元怎么代表无穷阵列”这件事别想当然侧面用的是 Floquet 周期性边界。垂直入射时相位差为零斜入射时需要按 kx k0 sinθ cosφ、ky k0 sinθ sinφ 把相位填进去。上下两个面不用周期性边界而是用 Periodic Port——这是算透射谱的关键。端口 1 设置入射平面波并打开 Excitation端口 2 作为透射输出端口COMSOL 会自动提取 S11 和 S21。最常见的一个错误是在周期边界里填了 kx、ky但端口设置里没有同步填或者填了两个不一致的参数。周期边界和端口必须共享同一个入射波矢量参数否则模式匹配会出错得到的 S 参数数值完全不靠谱。另一个常见问题是端口方向反了现象是 S21 恒等于 1、谱线没有谷或者 S11 和 S21 的角色对调。2.3 端口模式阶次亚波长单元最省心但别掉以轻心当周期 P 小于波长时所有高级衍射级次都变成倏逝波端口上能传播的模式只有零阶此时透射率用 |S21|² 没问题。一旦波长缩短到 P 以下比如周期 600 nm、波长扫到 550 nmCOMSOL 的 Periodic Port 会给出多个端口模式S21 对应的是零阶还是高阶需要你自己确认。处理方法是在端口设置里查看每个模式对应的切线波矢零阶模式的波矢切向分量为 (0,0)。手动限制模式数量也可以但一定要确保包含所有传播级次。透射率也应汇总所有传播的零阶模式功率而不是只看一个 S21。3. 扫透射谱从建模到 S21 提取的完整操作顺序3.1 参数化扫描用波长而不是频率光学里的人习惯看波长COMSOL 的频域接口习惯用频率。最顺手的做法是在全局参数里定义 λ然后让接口频率等于 c0/λ。参数化扫描时扫 λ比如先粗扫 800 nm 到 1300 nm、步长 10 nm。扫描完成后再围绕透射谷加密步长取 2 nm 甚至 0.5 nm。粗扫步长不建议一开始就设 1 nm因为周期性结构的共振对网格和端口设置很敏感先粗扫找出异常区域再细扫确认形态能节省大量时间。如果粗扫曲线看起来有锯齿或剧烈抖动不要直接跳进细扫先回去查 Floquet 相位和网格多半是设置问题而不是物理现象。3.2 透射率怎么取S21 的平方附带一条能量守恒检查一维绘图里直接画 abs(S21)^2 就是透射率。反射率用 abs(S11)^2。对于无损材料T R 应该接近 1。如果明显小于 1不是有吸收而是端口模式没取全或能量漏进了高阶衍射级次。有人喜欢直接用域内的坡印廷矢量积分算透射功率也可以但要注意积分面面积和周期单元的面积关系还得考虑 Floquet 端口本身的归一化。既然用 Periodic Port直接读 S 参数最省事Port 内部已经处理了功率归一化。只要基底和端口设置正确S21 就是可靠的。3.3 网格和求解器决定谱线可信度的两件事网格策略可以概括为颗粒内部细颗粒表面更细空气区域粗。圆柱域内最大网格尺寸按 R/6 左右控制20 nm 上下空气区域按波长除以 8 到 10 控制。如果共振位置的局域场增强明显再对这个区域做一次局部网格加密。求解上用频域直接求解器通常就够了PARDISO 是我习惯的选择。参数扫描本质上是一个波长一个波长地求解每个波长之间没有物理连续性所以不要指望“用上一个解做初值”能带来太明显的加速但它确实可以缩短部分波长点的迭代时间。网格收敛检验一定要做取几个关键波长点对比粗网格和细网格下的透射谷位置如果谷位偏移超过 1 nm网格还要加密。4. 多极子展开的核心位移电流上的体积分以及 COMSOL 实现4.1 先理解“多极矩”到底在描述什么入射光在颗粒中激励出振荡电流。远场观测者看不到电流分布的每一处细节只能看到这个电流的“整体矩”。电偶极矩是电荷重心整体振荡磁偶极矩来自电流沿某个方向形成的环流电四极矩则是两个反向偶极的对称形变环偶极矩对应一种更隐蔽的环形电流排列它在远场的行为和电偶极相近但物理机制不同。对高折射率介质颗粒来说磁偶极共振和电偶极共振是两种最常见的低阶响应。两个共振还会互相叠加产生明显的定向散射效应——也就是前向散射增强、反射减弱的所谓惠更斯行为。但如果你只看透射谱很难区分这个透明窗口到底来自“磁偶极共振”还是“电偶极与磁偶极干涉”这时多极子分解就是唯一能直接给出答案的工具。4.2 从电场到多极矩只需要两类积分在 COMSOL 中提取多极矩并不需要额外物理场接口。核心思路是对颗粒内部做体积分被积函数是感应电流密度 J。对于非磁介质感应电流来自极化电流Jx -i*omega*eps0*(epsilonr_mat - epsb)*Ex Jy -i*omega*eps0*(epsilonr_mat - epsb)*Ey Jz -i*omega*eps0*(epsilonr_mat - epsb)*Ez这里的 epsilonr_mat 是颗粒材料接口中的相对介电常数变量epsb 是背景相对介电常数。所谓“位移电流”在频域里就是 iωD 中与背景不同的部分。把 epsilonr_mat 这个材料变量直接拿来用空气区域会自动等于 epsb所以 J 在背景中自然为 0后面即使不小心选到空气区域积分也不会引入额外贡献。建立积分算子 intop1区域选颗粒域。下面是几组最关键的多极矩表达式采用文献中常见的“辐射可比”写法磁偶极矩和环偶极矩里已经并入光速因子最终数值可直接与电偶极矩在同一比例上比较px intop1(Jx)/(i*omega) py intop1(Jy)/(i*omega) pz intop1(Jz)/(i*omega) mx intop1(y*Jz - z*Jy)/(2*c0) my intop1(z*Jx - x*Jz)/(2*c0) mz intop1(x*Jy - y*Jx)/(2*c0) Qxx intop1(x*Jx - (x*Jxy*Jyz*Jz)/3)/(i*omega) Qxy intop1(x*Jy y*Jx)/(2*i*omega) Qxz intop1(x*Jz z*Jx)/(2*i*omega) Qyy intop1(y*Jy - (x*Jxy*Jyz*Jz)/3)/(i*omega) Qyz intop1(y*Jz z*Jy)/(2*i*omega) Qzz intop1(z*Jz - (x*Jxy*Jyz*Jz)/3)/(i*omega) Tx intop1((x*Jxy*Jyz*Jz)*x - 2*(x^2y^2z^2)*Jx)/(10*c0^2) Ty intop1((x*Jxy*Jyz*Jz)*y - 2*(x^2y^2z^2)*Jy)/(10*c0^2) Tz intop1((x*Jxy*Jyz*Jz)*z - 2*(x^2y^2z^2)*Jz)/(10*c0^2)强调一句如果你没有把几何中心放在全局原点那么表达式里所有 x、y、z 都要换成相对粒子中心的局部坐标例如 xc x - x0。原点错了高阶项会混入偶极项环偶极子和电四极子会突然出现没有物理意义的巨大数值。4.3 环偶极子容易漏也不能乱扣帽子环偶极子的公式看起来最长但它并不是一个可有可无的高阶项。在某些超表面结构里窄带透明或非对称线型恰恰来自环偶极子的增强。不过环偶极子与电偶极子的远场辐射方向图非常接近只比较辐射强度时容易混淆。我见过不少人直接拿电偶极子的增强去解释一个本来就是环偶极主导的共振结论完全反了。要区分两者得看两个矩的相对相位或者观察电流分布在结构内部是否呈现出首尾相连的环形形态。在 COMSOL 里画出粒子截面的电流箭头图比只盯着透射谱猜测有用得多。4.4 单位、相位和归一化怎么让各阶矩能放进同一张图不同多极矩的量纲不同直接画 abs(px)、abs(mx)、abs(Qxx) 没法公平比较。需要把它们按辐射功率或归一化权重折算到同一尺度。一个简单做法是画无量纲权重W_ED abs(px)^2 abs(py)^2 abs(pz)^2 W_MD abs(mx)^2 abs(my)^2 abs(mz)^2 W_EQ k0^2*(abs(Qxx)^2 abs(Qxy)^2 abs(Qxz)^2 abs(Qyy)^2 abs(Qyz)^2 abs(Qzz)^2) W_TD k0^2*(abs(Tx)^2 abs(Ty)^2 abs(Tz)^2)然后把 W_ED、W_MD、W_EQ、W_TD 分别除以四者之和得到 0 到 1 之间的权重曲线。这个权重图可以直观展示每个波长下是哪个多极矩占主导。注意这是“辐射权重近似”不是严格的多极辐射功率精确值但用来判断共振归属已经足够可靠。还要记住一点这些多极矩都是复数相位信息非常重要。当两个矩强度接近时透射谱上的峰谷往往取决于它们的相位差是相长还是相消。所以除了强度权重图建议同时把主导矩的相位随波长画出来。5. 把透射谱和展开系数放在同一张图里才算完成分析5.1 一图两轴S21 和多极矩权重摆在一起在 COMSOL 的 Results 里创建一个一维绘图组左侧轴放透射率 abs(S21)^2右侧轴放多极矩权重或 log10 之后的多极矩强度。不同版本里添加副坐标轴的方式有差别但逻辑一致先画一条曲线再右键添加副坐标轴把另一组数据绑定到副轴。这样做的好处非常明显。某一波长附近透射谷出现时如果权重曲线里 W_MD 同步冲上去而其他权重平缓那“这个谷来自磁偶极共振”就是板上钉钉的事如果透射谷本身不对称还把两个矩的峰重叠在一起说明这是模式干涉产物而不是单个模式的谱线。5.2 从谱线形状到物理结论的三种常见对应我总结出三种最常见的对应关系你可以用它做初步判断第一种透射谱出现对称的深谷同时某个多极矩在相同波长出现干净的峰这是单一模式主导的最典型信号。第二种透射谱出现非对称线型多极权重图上有两个矩形峰值部分重叠谷位往往落在两个峰之间这是两个模式相干干涉的结果。第三种透射谱整体平坦看不出明显谷但某个多极矩在某个频率明显增强这往往是暗模式或高 Q 模式对零阶透射几乎不贡献需要结合近场分布或非线性响应来确认。前两种在常见介质纳米阵列里最容易遇到第三种多出现在高 Q 超表面中。如果你发现 S21 谱和多极矩谱“完全对不上”不要急着怀疑公式先检查是不是出现了第三种情况。5.3 当展开系数和透射谱对不上时先查这三个地方第一个是原点。多极矩是以粒子中心为展开原点定义的原点偏了会导致高阶矩被污染尤其是电四极子和环偶极子会莫名其妙出现大值。第二个是积分区域。积分算子一定要选在颗粒域不要把整个空气域选进去——虽然理论上空气区域没有极化电流贡献但数值噪声和端口附近的无用场可能会污染积分结果。第三个是基底。如果模型里有实体基底基底本身也会贡献感应电流这时直接对颗粒域积分只能得到颗粒的矩但透射谱反映的是整个层结构的响应。处理办法是把基底折射率并入背景或者把基底也作为积分对象并单独解释。提示做任何多极子分析之前先算一个最简单的验证把多极矩总辐射功率和 COMSOL 的远场积分或能量守恒结果对比。对不上时问题多半出在原点和网格而不是物理本身。6. 我的固定流程和几个常被忽视的检查点6.1 我的拆解顺序粗扫、定位、细扫、分解、验证我个人习惯的固定顺序是第一步跑粗谱10 nm 步长标记所有透射谷或异常区域第二步以谷为中心做细扫确认线型和极值位置第三步在谷附近提取多极矩画出权重图第四步把主导矩的共振波长和谷位对比看差值是否在线宽范围内第五步做能量守恒和截断验证确认结果可信。这个顺序能省下一半时间。很多人一上来就把网格加密、波长步长设到 0.1 nm结果跑一次要几个小时最后发现模型边界条件有问题全部白算。先粗后细等于给整个流程加了一道安全阀。6.2 几个容易翻车但常被忽视的检查点端口模式顺序是最容易翻车的。高级衍射级次出现后COMSOL 会按波矢大小给端口模式编号S21 可能不是零阶模式。你需要在端口设置里查看每个模式的切向波矢确认自己取的是 (0,0) 那个零阶项。Floquet 相位参数的一致性也值得复查。周期边界和端口必须引用同一个 kx、ky 参数不少人在垂直入射时直接写成 0后来改成斜入射时只改了边界没改端口结果谱线形态完全错误。材料色散是第三个隐蔽点。如果你用固定折射率跑完再换成实际材料的色散数据共振峰位置可能移动几十纳米这通常不是物理问题只是模型没同步更新。最后说一个我自己的习惯每次跑完一组参数我都会顺手把透射谷位置、对应的主导多极矩、以及该矩的相对相位记录在一个表里。信息量最大的往往不是某个矩的强度而是两个矩之间的相对相位。后续如果要解释反常透射、干涉相消或者做参数优化这个表能帮你快速定位是哪些项在互相干涉。多极子分解一旦融入你的日常仿真流程透射谱就不再只是一条曲线而是一套能拆开来自我解释的语言了。