做超表面和纳米光子学仿真这几年，我拿到一个透射谱最怕听到一句话：“帮我解释一下这个峰是哪来的。”用 COMSOL 扫参数出透射谱只是第一步，要把透射峰和物理机制对上号，十有八九绕不开多极子展开。这篇文章围绕周期性结构里最常用的操作流程：模型搭建、Floquet 边界设置、透射谱计算、多极子分解和可视化，讲讲我在实际操作中踩过的坑和验证过的做法。适合正在做超表面、人工电磁材料、周期性天线阵列仿真的新手，也适合那种模型能跑、但解释不了结果为什么长这样的进阶玩家。

1. 为什么透射峰必须靠多极子展开归因

1.1 透射谱只能告诉你“在哪”，多极子才能告诉你“为什么”

周期结构里最常见的仿真输出就是透射率随波长变化的曲线。正常情况下，扫一个参数范围，能得到若干个透射谷或者透射峰，比如法诺共振、电偶极共振、磁偶极共振、晶格模式耦合，这些都可能在同一条谱线上留下痕迹。但只看谱线本身，根本分不清哪个峰来自电偶极子、哪个峰来自磁偶极子、哪个峰是四极子贡献、哪个峰是暗模式被某个对称破缺悄悄点亮。

我做过的模拟项目X里，一个单胞里有三到四个介质柱，透射谱上出现了一个很窄的谷。第一次我以为是最低阶电偶极子导致的，后来把各阶偶极矩算出来，才发现电偶极子贡献在这个波段几乎为零，真正活跃的是一个磁四极子和一个电四极子的叠加。如果没有多极子展开这步，我大概率会拿着场分布图硬编一个故事。场图确实能看到涡旋，但涡旋对应的是磁响应还是电响应的四极分量，肉眼很难下结论。

所以我把多极子展开定位成“谱线归因”的必备工具。它能把一个复杂的散射过程拆成有限几项，每一项对应一个明确的电荷/电流振荡模式。对周期性结构来说，这个拆解还能直接和透射谱、反射谱、吸收谱对上账。

1.2 多极子展开到底在展开什么

很多做光学设计的人第一次接触多极子，是在米氏散射那一章。一个球粒子受到平面波照射后，散射场可以按矢量球谐函数展开，每一组球谐系数对应一个电多极子或者磁多极子。零阶是电偶极子，接着是磁偶极子、电四极子、磁四极子，依次类推。

对任意形状的粒子，也可以用笛卡尔多极矩来近似。核心思想是：在远比波长小的区域里，感应电流的积分效应可以用几个积分量代替。第一个积分量是总电荷振荡形成的偶极矩，第二个积分量是环形电流形成的磁偶极矩，第三个积分量是电流的空间分布不对称形成的四极矩。粒子一但比波长大，高阶项就开始变得不可忽略。

这些多极矩一旦算出来，每种模式的辐射功率就有解析公式。偶极辐射功率和频率的四次方成正比，四极辐射功率和频率的六次方成正比。这个频率依赖关系直接影响谱线的峰位和宽度。低频段高阶多极子天然压不过低阶项，高频段四极子甚至更高的八极子就有机会冒出来。所以多极子展开不是简单的数学展开，它实际上在告诉我们每一个共振模式在频谱上的“势力范围”。

1.3 周期结构里的坑：单粒子公式不能直接套

周期性结构和孤立粒子最大的区别，是单元之间存在多次散射。孤立粒子散射出去的波直接进入远场，但周期结构里，一个粒子辐射的波会被相邻粒子再次散射，形成晶格集体模式。严格处理这个问题，需要引入晶格求和或者周期格林函数。

实际操作中，很多人直接把自由空间的偶极辐射公式套到周期结构上，这在小占空比、弱耦合单元的情况下误差可以接受，单元间距一密，误差立刻凸显。我通常把周期结构的单元设计成“类孤立粒子”状态，比如相邻单元间距不小于四分之一波长，且中间没有强相互作用的介质连接。在这种前提下，用自由空间公式算出的多极子功率依然有参考价值，但解释结果时必须说明这是单粒子近似。

如果单元非常密，或者晶格模式参与了共振，我会把透射谱计算作为主要锚点，把多极子展开作为辅助分析，而不是反过来。就像用放大镜看细节之前，得先用尺子确认整体尺寸没错。

2. 模型搭建与边界条件的关键设置

2.1 用电磁波频域接口搭最小可行模型

COMSOL 里做这类仿真，最常用的是电磁波，频域接口，也就是 ewfd。三维模型用全矢量求解，二维模型可以模拟无限长柱体的横截面，二维等效相对简单。我建议周期结构定性分析先用二维跑通流程，再上三维。二维模型里粒子是无限长柱体，多极子展开公式要做对应修改，但边界条件、端口设置、透射谱提取的思路完全一致。

建模时我喜欢把单胞做成一个方形盒子，粒子放在中心。盒子周围四个侧面设置周期性条件，上下两个面设置入射和出射端口。材料参数要注意色散，尤其是介质柱在可见光波段折射率随波长变化很剧烈，直接用常数折射率会在共振位置产生明显偏差。

几何上用布尔操作把粒子域和背景域分开，这样后续积分多极矩时可以单独选中粒子域。背景域如果是空气，直接用默认真空属性就行。如果背景是衬底或者折射率匹配液，需要额外定义背景域的材料，并把多极矩积分域严格限定在粒子内部。一个重要提醒：背景介质中也会存在感应电流吗？如果是纯介质、无损、无增益，均匀背景中平面波不会产生净体电流，所以只对粒子域积分是可以的。但若粒子周围有金属膜或导电层，就必须包含这些区域的电流。

2.2 Floquet 周期边界和端口激励的设置要点

三个周期性侧面用 Floquet 周期边界条件。设置时最关键的是定义 k 向量。正入射时 k 向量的 x、y 分量为零，斜入射时要按入射角和方位角填入波矢分量。很多人斜入射算错，是因为 Floquet 边界上给定的 k 分量和端口激励角度没对齐。

我习惯把入射端口设置在最底部的面，采用“端口”边界条件，并把模式指定为平面波。出射端用另一个端口收集透射功率。端口大小和单胞一致，这样透射谱对应的是零级衍射波。如果波长低于晶格常数，会出现更高阶衍射模式，这时候端口设置会更复杂，要加多个端口对应多个衍射级。先确定工作波长是否在一个衍射级以下的区间，能省掉很多麻烦。

端口设置完成后，检查一下输入功率的方向。COMSOL 端口边界会返回 S 参数，但我建议不要只依赖 S 参数。直接用积分算子统计出射面功率，再和入射功率相除，得到的结果通常更直观，也更方便后续做能量守恒检查。

2.3 网格和求解器的坑

网格这块，最大的问题是粒子内部网格密度不够。表面共振模式对粒子边界非常敏感，边界层网格至少要留两到三层，而且要用三角形或四面体。波长分辨率的经验值是每波长至少六个二阶单元，但共振位置附近要更密。我做参数扫描时，会在共振峰附近加密网格重新验证一次，如果峰位移动超过百分之几，说明网格还没收敛。

求解器方面，频域扫描通常用直接求解器比较稳。周期性边界条件耦合了不同边界上的节点，迭代求解器经常收敛很慢。内存允许的话，直接选 MUMPS，它对这类问题表现稳定。

模型搭建阶段最容易犯的错，是忽略单元的镜像对称设置。如果入射光是圆偏振，或者结构本身有旋转对称性，可以利用对称边界把模型缩小到四分之一甚至八分之一。但要注意对称边界上的场对称类型必须和入射场一致，设置反了，计算结果会完全错误。我建议第一次做周期性结构时先别用对称缩减，跑全单胞，确认没问题后再优化。省下来那点算力不值得冒出错的风险。

3. 多极子展开的数据处理流程

3.1 从电场到感应电流的换算

COMSOL 频域求解完成后，后处理里能拿到复电场分量 Ex、Ey、Ez。多极子展开需要的是感应电流密度 J，而不是总场本身。均匀无源介质中，感应电流密度和电场的关系是：

J = i * omega * epsilon0 * (epsilonr - 1) * E

这里 epsilonr 是粒子的相对介电常数，epsilon0 是真空介电常数，omega 是角频率。这个表达式里已经默认了谐波因子是 exp(i * omega * t)。如果你定义变量时发现实部和虚部跟理论预期对不上，大概率是时谐约定反了，把 i 改成 -i 再试。

注意一个细节：这个公式里的 epsilonr 是复数。有损耗的材料，epsilonr 的虚部不为零，必须用完整复数值，不能只取实部。否则多极子积累出来的功率会偏小。

对于金属粒子，不能用简单折射率，最好直接用光学常数。金属中诱导电流主要在表层，网格要加密到趋肤深度能满足。多极子积分对表面电流分布很敏感，电流积分域哪怕只缺一小块边缘，高阶多极矩都会出现假值。

3.2 用集成算子实现多极矩的实时计算

在后处理前，先在组件定义里添加积分算子。例如创建一个 intop1，作用范围选择粒子域。之后定义自用变量，就能在参数扫描结束后直接输出多极矩。

我用过的常用变量如下：

名称: Jx_eff 表达式: i*ewfd.omega*epsilon0_const*epsr_eff*ewfd.Ex 单位: A/m^2

其中 epsr_eff 是材料相对介电常数。为了跟随材料定义，可以直接用 COMSOL 内置的材料属性，也可以手动定义色散函数，看你的材料库设置。

然后定义多极矩积分：

px = intop1(Jx_eff)/(-i*ewfd.omega) py = intop1(Jy_eff)/(-i*ewfd.omega) pz = intop1(Jz_eff)/(-i*ewfd.omega)

这里的负号取决于时谐约定，我在第三部分已经提醒过。多极矩的单位是库仑·米。

磁偶极矩的表达式稍复杂一点，需要叉乘位置矢量和电流：

mx = intop1( (y*Jz_eff - z*Jy_eff) ) / (2*clight_const) my = intop1( (z*Jx_eff - x*Jz_eff) ) / (2*clight_const) mz = intop1( (x*Jy_eff - y*Jx_eff) ) / (2*clight_const)

坐标原点必须放在粒子对称中心。如果模型几何原点不在粒子中心，计算前要把积分表达式的坐标偏移量减掉，否则会引入虚假的高阶分量。

电四极矩的表达式更繁琐，以 Qxx 为例：

Qxx = intop1( 2*x*Jx_eff - (2/3)*(x*Jx_eff + y*Jy_eff + z*Jz_eff) ) / (-i*ewfd.omega)

并需要类似 Qxy、Qyy、Qyz 等其他分量。实际使用中，我一般先把 px、py、pz、mx、my、mz、Qxx、Qxy、Qyy 这组算出来，检查量级，再看哪些分量占主导。

3.3 用全局计算批量导出数据

参数扫描结束后，我喜欢用“派生值”—“全局计算”直接输出所有多极矩变量随波长变化的表格。设置表达式时用 abs(px) 或者 abs(px)^2，前者是偶极矩振幅，后者正比于辐射功率。

如果 COMSOL 版本支持的话，直接在全局计算里写：

abs(px)^2 + abs(py)^2 + abs(pz)^2

能立刻得到电偶极子贡献的强度因子。同样的方式处理磁偶极子。四极子需要把六个独立分量的模平方求和，再加交叉项。在计算功率前，把所有矩的平方项导出，后续在 Origin 或者 Python 里再渲染成谱线。

这一步最值得注意的问题是：整数算子的积分范围必须固定。参数扫描时建议使用同一个网格，不要在每个频率点重新剖分网格，否则积分值会发生小幅跳变，虽然不影响定性判断，但在精细结构上会多出噪点。

4. 透射谱计算与多极子结果的对账

4.1 透射率的两种计算方式

透射率计算有两套并列方案。第一套是直接用端口 S 参数，读 S21 的幅值平方。第二套是在出射面用积分算子统计能流，除以入射功率。两者结果理论上一致，但在斜入射、高损耗、或者周期边界没对齐的情况下，S 参数会出现额外的归一化问题。

我的习惯是两套都算。参数扫描时把两种透射率都输出，如果两条曲线偏离超过百分之一，说明模型设置有问题，或者端口激励没有和 Floquet 边界匹配一致。

用能流法时，积分表达式如下：

P_out = intop_out( nx*ewfd.Poavx + ny*ewfd.Poavy + nz*ewfd.Poavz )

其中 intop_out 定义在出射端口面上，n 是外法向分量。COMSOL 的 Poavx 字段表示时均能流的 x 分量。如果出射面法向指向 z 方向，表达式可以简化成 intop_out(ewfd.Poavz)——前提是坐标轴方向已经对齐。

入射功率直接按解析公式算：

P_in = |E0|^2 * A_cell * cos(theta) / (2 * eta0)

这里 A_cell 是单胞面积，theta 是入射角，eta0 是自由空间波阻抗。正入射时 cos(theta) 等于 1。这比从端口读参考功率要透明得多，因为你清楚每一步除以了什么东西，不容易出现归一化错误。

4.2 用能量守恒抓隐形错误

透射率 T、反射率 R、吸收率 A 三者之和应该严格等于 1。周期性结构里如果有金属或吸收介质，吸收率不为零，这正好可以用来看模型自洽性。

反射率用同样方法在入射面统计反向能流得到。吸收率可以直接用全局变量ewfd.Qrh或者用体积分统计损耗功率密度。三者之间的偏差如果超过 2%，先去查网格，再去查边界条件。我遇到过一次偏差很大的情况，最后发现是出射端口的边界条件选成了散射边界，和周期端口不兼容，导致部分能量被边界反射回来。

能量守恒检查不只是为了确认透射率准不准，它还能帮你判断多极子展开的近似范围。用自由空间公式算出的全部多极子辐射功率，如果想要和透射谱“消失的能量”对账，必须考虑晶格耦合的修正。所以我不指望多极子总功率严格等于 1-T-R-A，只要求趋势一致，峰位对齐。如果多极子总功率的峰值波长和透射谷完全错开，那一定是处理公式或者积分范围出了问题。

4.3 用多极子功率和总散射功率做交叉验证

更严格的做法是，在 COMSOL 里再算一次总散射功率。用散射场在包围粒子的一个虚拟球面上做能流积分，得到总散射功率 P_sca。然后把电偶极、磁偶极、电四极这几项相加，和 P_sca 对比。

如果结构占空比很小，单元间耦合弱，理论上多极子各项相加应该接近总散射功率。如果差得远，通常有两个原因：一是高阶项没有取够，比如八极子甚至更高的贡献不可忽略；二是笛卡尔多极矩公式里的积分次幂不完整，比如缺少了某个分量。

我常用的对比点不是某个单一波长，而是整个扫描区间的积分。把每个波长的 P_sca 对所有波长积分得到总散射能量，再把多极子各项对所有波长积分，两者比值能反映展开的完整度。完整度低于百分之七十，我会提高展开阶数或者改用更严格的晶格多极子方法。

5. 可视化分析与常见问题排查

5.1 多极子谱线和占比堆叠图的画法

计算完成后，可视化是最容易出彩也最容易忽视的环节。不做多极子分解之前，你的透射谱只有一条二维曲线。做完展开后，你有至少五到六条曲线：电偶极、磁偶极、电四极、磁四极、多极子总和。把它们放在同一张图里，取对数纵坐标，往往能直接看到几个共振峰分别由哪个通道贡献。

堆叠图更直观。把每个波长下各多极子贡献的功率归一化，按百分比堆成面积图，一眼就能看出某一个窄带共振中哪个模式占绝对主导。我在做工作汇报时经常用堆叠图，因为非专业观众也能看懂“颜色变了就是模式切换了”。

画图建议用 Origin 或者 Python，导出 CSV 后在外部处理。不太建议直接在 COMSOL 里做太精细的排版，COMSOL 的画图引擎更适合检查和初稿。

5.2 把多极矩和远场方向图、电场分布叠在一起看

多极子展开虽然给出了每个模式的强度，但没有给出相位信息。强度决定共振峰的高度，相位决定干涉方向图。为了看清某个透射谷是电磁偶极干涉还是两个电四极子的对称劈裂，要把多极矩的复数振幅导出来，再算它们之间的相位差。

比如电偶极子和磁偶极子同时激发时，若它们的散射场在远场某方向干涉相消，就会在方向图上出现一个凹陷，同时背向散射变弱。这个现象用单纯谱线看不出来，但把多极矩的相位整理成矢量图，或者直接画出散射方向图，就非常直观。

方向上，COMSOL 的远场节点能直接输出远场辐射图。选择在模型周围定义虚拟球面，用ewfd.Efar之类的变量计算远场，然后画极坐标图。远场结果里能看到多极子干涉的特征。缺点是周期性结构在非垂直入射时，远场方向图会被栅瓣干扰，解读时要先确定所有高阶衍射级的位置，别把栅瓣当成物理模式。

5.3 常见问题速查

5.4 周期性结构高频失效与替代方案

如果工作波长小于晶格常数，单个单元周期内会出现多个衍射模式。零级透射波不再代表全部透射能量，多极子展开的单粒子近似也会失效，因为高阶 Floquet 模式会在单元之间传输能量，不能简单等效成孤立辐射源。

此时我通常改用两种方案。一是在 COMSOL 中显式建立多个端口，对应每个衍射级，把透射谱扩展成每个衍射级的谱线集合。二是自定义脚本实现晶格多极子求和，将每个单元的电流密度作为源，与周期格林函数卷积。后一种方案的代码量明显增加，但能处理强耦合和暗模式。

做周期性结构仿真这么多年，我的经验是：先把“单粒子近似、弱耦合”的模型跑熟，得到清晰的物理图像，再逐步增加复杂度和耦合效应。一上来就追求全严格处理，往往会被数值噪声掩盖真正的物理规律。多极子展开的好处恰恰在于它给了一个清晰的、可解释的中间层，让复杂问题不至于完全黑箱。

如果你只打算在项目里解决一个窄带透射峰的来源问题，前面的做法已经完全够用。如果还需要解释陌生成分或高阶混杂模式，就得结合场分布、远场方向图和多极矩相位做联合诊断。这一步做完，基本就能把一个纯仿真数据转化为有说服力的物理结论了。