COMSOL纳米圆柱多极散射分析：建模、积分与结果解读

在COMSOL里把纳米圆柱的散射谱算出来很多人习惯性地截一张彩色的电场分布图然后就没有下文了。实际上那张频谱图上的每一个共振峰背后都对应着一种或几种辐射模式——如果你不做多极散射分析就说不清它到底是电偶极、磁偶极还是四极主导更没法接着回答“怎么让前向散射更强”“怎么把背向散射压下去”这类实际问题。这篇文章我就完整复盘一下在COMSOL里做纳米圆柱多极散射分析的流程包括理论背景、模型搭建、多极矩积分实现以及我在实际计算里踩过的一些坑。1. 多极展开从“一团光斑”到“模式清单”1.1 为什么球能算解析解圆柱却只能靠数值球形颗粒的散射问题在1908年就被Gustav Mie用分离变量法严格解出来了。入射平面波按球谐函数展开散射场也按球谐函数展开最后每一阶分波的系数都能用Bessel函数和Hankel函数表示。这就是经典的Mie理论。它最大的好处是结果天然按电磁多极子分类——每一个共振峰对应哪一阶电多极TM或磁多极TE模式一目了然。但到了圆柱形颗粒这里事情就没这么简单了。无限长的圆柱确实可以用二维柱坐标下的波函数展开得到一组“类Mie”的解析表达式可是实际研究里的“纳米圆柱”往往是有限高的三维结构比如纳米盘、纳米柱、圆盘天线。这些结构的边界条件不匹配任何可分离变量的坐标体系只能靠数值方法求解。问题在于数值方法算出来的结果通常是“总场”“总散射截面”它不告诉你到底哪种模式在共振。1.2 多极子展开的物理图像多极子展开的物理图像其实很直观。任意一个散射体在光场作用下产生的电荷位移和电流分布都可以看成是一系列基本辐射源的叠加电偶极子、磁偶极子、电四极子、磁四极子以及更高阶的项。总散射场就是这些多极子辐射场的相干叠加。更关键的是在远场球面上各阶多极子辐射场是相互正交的所以总散射功率刚好等于各阶多极子辐射功率的标量求和不需要考虑交叉干涉项。打个比方这就像你录了一支交响乐如果只给你一个总音量表你很难分辨当前是弦乐还是铜管在响。多极子展开就相当于把每个乐器单独拉出一条音轨虽然声音是同一支乐队发出来的但你可以分别检查每一轨的响度甚至做混音调整。在纳米光学中这就是“模式归属”能力。共振峰会出现在某个波长附近但到底是谁在共振必须把散射功率按多极子通道拆开才能说清楚。1.3 高折射率纳米圆柱里藏着磁共振为什么要单独强调纳米圆柱因为圆柱/圆盘结构在“磁响应”方面有特殊价值。金属纳米球主要支持电多极等离激元共振磁响应很弱但高折射率介质颗粒硅、锗、二氧化钛能在可见或近红外波段支持强磁偶极共振。这个磁共振来源于颗粒内部位移电流形成的环形分布等效于一个磁偶极子。纳米圆柱因为几何自由度比球多直径和高度可以独立调节磁偶极共振与电偶极共振的频谱位置更容易调谐重叠。当电偶极和磁偶极共振重叠时会产生所谓的Kerker效应在一个特定波长附近前向散射显著增强背向散射被抑制。这种定向散射特性在天线、光伏、传感领域都有应用价值。而要在仿真中确认这种效应是“电偶极加磁偶极”共同造成的就必须做多极分解。2. 模型搭建前必须拍板的三个选择2.1 二维模型还是三维模型先看圆柱的纵横比在动手建几何之前第一个要决定的事情是用二维还是三维模型。如果是无限长的纳米线且光源垂直入射可以用二维模型近似把问题分解成TM偏振电场平行柱轴和TE偏振磁场平行柱轴两个独立问题网格量小一个量级。但“纳米圆柱”这个说法在大多数文献里指的是有限高的三维圆柱或圆盘这时候必须用三维模型。一个粗略的判断标准圆柱高度如果大于两倍直径可以尝试二维线模型高度和直径在同一量级或者高度明显小于直径也就是纳米盘必须建三维。本文以下都以三维有限高圆柱为例二维模型的多极展开逻辑不同它应该用二维柱谐波展开而不是三维球谐波展开公式要单独查。2.2 材料库里的数据不总是可靠的COMSOL材料库确实带了不少光学材料数据但用之前一定要检查数据的来源、波长范围和虚部。金属纳米颗粒最常用的是Johnson and Christy在1972年发表的银、金介电常数数据它覆盖约0.1到6 eV能量范围在可见光波段基本够用。但要注意薄金属膜和纳米颗粒的实际介电常数会因为晶粒尺寸、氧化层、表面散射效应而偏离块体数据仿真结果和实验对不上时多半是材料参数的问题不是求解器的问题。高折射率介质方面硅的数据要特别小心。晶体硅在近红外约1.1微米以上是透明低吸收的折射率约3.5但在可见光波段吸收明显复折射率实部和虚部都随波长剧烈变化。如果你仿真的波长范围是400 nm到900 nm一定要用插值形式的波长相关折射率数据不能只填一个常数。我的习惯是优先用Palik或文献实验数据把它们整理成“波长-实部”“波长-虚部”两张插值表导入COMSOL而不是直接用材料库默认的单一折射率。2.3 总场与散射场多极分解必须用散射场公式在COMSOL的“电磁波频域”接口里默认求解的是总场。对于散射问题我强烈建议改为求解散射场。原因有两个。第一散射场公式下入射平面波是解析已知的数值求解只针对散射场这样计算域边界上的反射误差小很多第二多极矩积分需要用的量是散射体内部的“额外极化电流”它等于总电场乘以介电常数增量而散射场公式天然把背景场和散射场分开后处理时表达式清清楚楚。当然这里有一个非常容易踩的坑COMSOL中选用散射场公式后默认的后处理变量ewfd.Ex往往指的是因变量也就是散射场不是总场。你要算总场得写成ewfd.Ex ewfd.Ebx之类背景场变量的具体名称因版本而异常见的是ewfd.Ebx、ewfd.Eby、ewfd.Ebz也有直接用Ebg的。多极矩积分公式里那个电场必须是总电场如果误用散射场结果会完全不对而且看起来还很合理很难排查。3. 几何、PML与背景光的落地细节3.1 几何尺寸的建议实际建模时我通常把纳米圆柱放在坐标原点底面在z0高度为H直径为D。周围包一个球形空气域最外面再包一层球形完美匹配层PML。空气域的半径没有一个严格公式但经验值是从圆柱表面到PML内边界的距离至少要大于最大工作波长的二分之一。如果你扫描的波长范围是400 nm到900 nm那这个距离就按900 nm对应的半个波长来留也就是至少450 nm。然后PML厚度取最大波长的四分之一到三分之一也就是230 nm到300 nm左右。这里要注意一个常见误区空气域半径不是看颗粒尺寸而是看波长。颗粒直径往往只有100到200 nm但空气域要到600到800 nmPML还要再加两三百纳米。域开得小了PML会反射数值噪声导致散射谱上出现莫名其妙的毛刺。我常用的一个参考配置如下表工作波长范围圆柱直径/高度空气域半径颗粒中心到PML内边界PML厚度400–900 nm100–200 nmλ_max/2 颗粒半径λ_max/4约230 nm600–1200 nm200–300 nmλ_max/2 颗粒半径λ_max/4约300 nm可见光窄带如632.8 nm100 nm左右约450 nm约180 nm3.2 PML的设置要点选中PML域后在“完美匹配层”节点里把几何类型设为“球形”或“圆柱”这取决于你外层空气域的形状。球形空气域配球形PML是最省心的组合。COMSOL默认的PML参数层数、比例因子对大多数光学散射问题都够用但如果你发现远场方向图出现不自然的条纹或者散射功率在某些波长出现跳变第一反应应该是把PML的比例因子调一调或者干脆把PML厚度加倍做对比测试。PML的另一个要点是网格。COMSOL的PML内部有特殊的坐标缩放算法网格不能太粗但也不需要特别细。我的做法是给PML域单独指定“扫掠网格”或“映射网格”径向方向至少划5层以上切向加密到与相邻空气域网格一致。用自由四面体网格直接剖PML也能算但相比之下反射会大一点收敛性不如扫掠网格稳定。3.3 定义入射平面波在散射场公式下入射平面波用“背景场”节点定义。假设一束x方向偏振的光沿z轴正方向入射背景电场表达式可以写成Ex E0 * exp(-i*k0*n_bg*z) Ey 0 Ez 0在实际输入时E0设为1 V/m即可因为截面计算时入射强度会被归一化掉。k0是自由空间波数n_bg是背景介质折射率如果背景是空气n_bg就是1。这里有一个新手特别容易犯的错直接把背景波矢写成2*pi/lambda0忽略了背景介质折射率。如果颗粒是放在水或者油里背景折射率不再是1这一项错了后面的所有结果都会系统性偏移。另外偏振方向、传播方向和颗粒轴向的相对关系会显著影响多极模式的激发效率。圆柱绕轴旋转对称但线偏振平面波不具备这个旋转对称性所以不能把模型简化成一个二维轴对称问题。不过如果入射光沿轴向入射圆柱关于x-z平面和y-z平面各有一个镜像对称面可以借助理想导体PEC和理想磁导体PMC边界条件把模型缩小到原来的四分之一大幅节省内存。热词里有人搜“COMSOL 4内存”跑这种三维频域扫描确实吃内存用对称性是性价比极高的手段。4. 多极矩计算在COMSOL里直接做体积分4.1 从极化强度出发的多极矩定义多极矩的计算方法在文献里有好几个版本单位制和时谐约定不同公式长得也不太一样。我使用的是一组以极化强度 ( \mathbf{P} \varepsilon_0(\varepsilon_r-\varepsilon_b)\mathbf{E} ) 为积分量的公式在COMSOL里比较好落地。电偶极矩[ \mathbf{p} \int_V \varepsilon_0(\varepsilon_r-\varepsilon_b)\mathbf{E},dV ]磁偶极矩[ \mathbf{m} \frac{i\omega}{2c}\int_V \mathbf{r}\times\mathbf{P},dV ]电四极矩张量[ Q_{\alpha\beta} \frac12\int_V \left[ r_\alpha P_\beta r_\beta P_\alpha - \frac{2}{3}\delta_{\alpha\beta},\mathbf{r}\cdot\mathbf{P}\right] dV ]这些公式里最关键的一个物理量是 ( \varepsilon_r-\varepsilon_b )也就是颗粒材料相对介电常数与背景介质相对介电常数之差。注意这里用的是复介电常数所以电场越强、介电常数对比越大极化就越强多极矩也就越大。用这个差值而不是用颗粒自身的绝对介电常数可以确保积分只在颗粒体积内有非零贡献。4.2 COMSOL时谐约定对公式的影响为什么我特意说“我使用的是一组公式”因为COMSOL内部默认的时谐因子通常是 ( e^{i\omega t} )而许多光学文献习惯用 ( e^{-i\omega t} )。这两种约定下“入射波沿方向传播”的表达式会差一个符号极化电流密度 ( \mathbf{J}\partial\mathbf{P}/\partial t ) 也会差一个虚数符号。如果你的公式是从某个文献直接抄来的务必确认它的时谐约定否则计算出来的多极矩相位会错。好在纳米圆柱多极分解最关心的是各通道的散射功率而功率公式里用的是 ( |\mathbf{p}|^2 )、( |\mathbf{m}|^2 )虚数符号对模方没有影响所以只算功率的话这个坑不明显。但如果你后续要做多极子之间的干涉分析比如算Kerker条件里的相位关系或者要重建远场方向图符号就很重要。4.3 积分算子加全局计算的操作在COMSOL里实现多极矩积分我用的是“组件耦合”里的“积分”算子而不是每次手动选域。具体做法是在“定义”节点下添加一个“积分”算子比如命名为int_np然后选择圆柱所在的域。接下来在“全局计算”里直接写积分表达式它会自动遍历频域扫描的每一个波长点输出一个结果表。假设你已经在变量或表达式中定义了epsr_np颗粒相对介电常数epsr_bg背景相对介电常数Px eps0_const*(epsr_np-epsr_bg)*ewfd.Ex注意这里要用总场如果后处理变量默认是散射场要补充背景场项那么电偶极矩的三个分量可以写成px int_np(Px) py int_np(Py) pz int_np(Pz)磁偶极矩的分量会稍微长一点例如mx int_np(i*omega/(2*c0)*(y*Pz - z*Py)) my int_np(i*omega/(2*c0)*(z*Px - x*Pz)) mz int_np(i*omega/(2*c0)*(x*Py - y*Px))然后在全局计算里定义散射功率表达式。我带一个电偶极通道作为例子P_ED c0*Z0_const*k0^4/(12*pi) * (abs(px)^2abs(py)^2abs(pz)^2)Z0_const是真空波阻抗部分版本可能写作Z0你用之前先确认一下常数名。电四极通道的功率公式系数不一样用我之前写的张量分量求模方求和。具体的功率公式系数我建议以你选定的文献为准公式版本很多最重要的是保证自洽。4.4 用小球和Mie理论做自检多极分解的公式写进COMSOL后绝对不能直接拿去算圆柱。我每次新建一个多极分解模型第一件事是最先做一个半径30 nm的小球验证。因为小球的Mie理论有解析解各阶分波贡献是已知的。用同一个积分公式算出小球的电偶极和磁偶极散射截面再和Mie理论对比如果对得上说明你的公式、符号、单位都是自洽的如果对不上先不要碰圆柱把问题查清楚再往下走。这个自检步骤看着简单但价值极高。多极分解公式涉及单位制、时谐约定、背景介质修正、变量名称任何一个环节错了结果都不可信。而Mie理论解析解就是最靠谱的标尺。小球半径取很小是为了保证偶极通道占绝对主导高阶项可以忽略这样对比时只看前两项就够了。5. 网格、频域扫描与收敛性检查5.1 网格尺寸的经验标准网格划分直接决定纳米圆柱散射仿真的成败。我按三个区域分别控制网格密度颗粒内部、颗粒周围的近场区域、外层空气域和PML。颗粒内部网格要分辨电场在颗粒内的空间变化尤其是高折射率介质和金属。经验上最大单元尺寸取 ( \lambda_{\text{np}}/10 ) 到 ( \lambda_{\text{np}}/15 )其中 ( \lambda_{\text{np}} \lambda_0/\text{Re}(n_{\text{np}}) )。对于硅颗粒在近红外波段这个值大约是几十纳米。金属颗粒麻烦一点金和银在可见光波段的趋肤深度只有十几到几十纳米电场主要集中在表面需要在颗粒表面加若干层边界层网格第一层厚度要小于趋肤深度。颗粒周围的近场区域网格取 ( \lambda_0/8 ) 左右就够了。外层空气域和PML的网格可以更粗一些PML内部用扫掠网格保证径向方向有足够的层数。一个常见问题是颗粒边缘或底角处曲率很大电场会在这个位置出现剧烈变化要在几何建模时做适当圆角或者在圆柱边缘添加局部细化网格。5.2 频域扫描的设置技巧频域扫描在“研究”节点里用参数化扫描实现。把波长lambda0设成扫描参数一开始用5 nm到10 nm的步长粗扫找到共振峰位置后再用2 nm甚至1 nm的步长在峰附近加密。COMSOL的“参数化扫描”支持一个表达式列表可以直接写400[nm] 450[nm] 500[nm] 550[nm] ...也可以用范围表达式加列表的组合。但注意材料折射率插值表要以lambda0为插值自变量参数化扫描的变量名和插值表自变量必须一致否则材料数据不随波长变化结果会是一堆完全错误的值。由于每个波长点都需要重新组装并求解一次矩阵波长点越多计算时间越长。三维模型上百个波长点用直接求解器跑一次普通台式机可能要跑小半天。我的建议是先用粗扫描摸清峰值位置再局部加密扫描而不是一上来就全波段高分辨率扫。5.3 三步收敛性检查做完一批结果后我习惯做这样三件事。第一把颗粒内部网格尺寸减半重新算共振峰处的散射截面看变化是否小于1%。如果变化大说明网格还不够密。第二把PML厚度翻倍或者把空气域半径加大看远场有没有变化。PML反射会造成远场轻微的振荡或者散射功率的虚假偏移。第三也是最容易被忽视的把“远场球面积分得到的总散射功率”和“各多极子散射功率之和”做对比。理论上这两个量应该相等如果差别超过百分之几说明多极分解公式有问题或者远场积分区域不够大、PML有反射需要回头排查。这三步做完结果才算可信才能进入物理分析环节。直接拿未验证的结果去写论文或做实验对比很容易得出完全错误的结论。6. 结果解读与常见坑6.1 多极散射谱线的归属判断多极分解的真正乐趣在结果解读。把全局计算得到的电偶极、磁偶极、电四极散射功率随波长的变化画到同一张图上你会看到一个总散射谱上的每个峰实际上是若干通道的组合。比如说一个直径120 nm、高度120 nm的硅纳米圆柱在可见光到近红外波段可能有两个明显的散射峰。其中一个峰对应的电偶极散射功率最高另一个峰来自磁偶极第三个峰则可能包含电四极和磁四极的贡献。把这几个通道的功率画成堆叠图共振峰的物理来源就一目了然了。如果是实验样品设计你看图就可以知道应该怎么调几何参数想让磁共振红移就把直径或高度加大想让电偶极和磁偶极重叠就同时调这两个参数。6.2 与文献对比时的单位、定义检查单对比文献数据时多极分解的坑非常隐蔽因为在不同文章里“多极子散射截面”的定义可能不一样。有的文章用多极子辐射功率除以入射强度得到的是“散射截面贡献”有的文章用的是多极子激发系数的模方还带有一个球谐函数归一化因子还有的文章把背景介质中的修正因子 ( n_b ) 拿走或加进来。直接对比数值之前先把定义查清楚否则你会觉得自己的结果比文献大好几倍或者小好几倍。我的检查顺序是单位制SI还是高斯、时谐约定、背景介质折射率、多极矩定义极化流还是电流还是位移电流、公式中的Z0和k0是否已换成背景介质中的值。这条检查单过一遍大多数对不上的问题都能解决。6.3 多极矩原点漂移造成的高阶模式假阳性这是一个我踩过很深、而且不太容易自己发现的坑。多极矩的数值严重依赖积分坐标原点的选择。如果原点不在颗粒中心即使原本纯粹的偶极散射也能在四极通道里分解出虚假的大分量甚至高阶矩超过低阶矩导致你误以为存在某个高阶共振模式。物理上这是因为各阶多极子之间在非原点坐标下并不正交平动会混合不同阶的贡献。所以在建模时一定要把圆柱放在坐标原点多极矩积分算子选中的域精确对应该域。如果因为几何装配问题导致颗粒不在原点可以新建一个坐标系把积分表达式里的坐标变量改成相对于颗粒中心的局部坐标。网格不对称也会造成类似的效应尽量保证关于颗粒中心的网格对称分布。6.4 消光截面、散射截面和吸收截面要分清很多人在初学时会混淆消光、散射、吸收这三者。消光截面等于散射截面加吸收截面。实验上常见的透射光谱测量得到的基本是消光截面贡献而暗场散射光谱测到的是散射截面。在仿真中散射截面通过远场坡印廷矢量的积分得到吸收截面则通过对颗粒体积内的功率损耗密度积分得到消光截面是两者之和也可以通过光学参量直接积分得到。COMSOL里吸收功率密度在不同版本有不同的内置变量常见的是ewfd.Qh、ewfd.Qrh、ewfd.Qmh分别对应电阻损耗、介电损耗和磁损耗。对金属纳米颗粒电阻损耗主导对半导体纳米颗粒介电损耗主导。实际计算时可以先把几种损耗都积出来看哪种机制贡献大这对理解材料在共振中的角色也很有帮助。6.5 绘制远场方向图时原点一定要放在圆柱中心COMSOL的“远场”节点会以计算域中某个点作为相位参考点来绘制远场方向图默认通常是全局坐标系原点。如果你建模时没有把圆柱放在原点或者远场参考点选在了别的位置画出来的远场方向图会有额外的相位倾斜前后向散射比也可能被扭曲。修正方法很简单确保圆柱中心在原点或者在远场节点设置中手动指定参考坐标等于圆柱中心。我已经见过很多次有人算出的“前向散射增强”结果其实只是参考点位置不对导致的假象。另外绘制远场方向图时检查一下你画的是正z方向还是负z方向的前向。入射光沿z传播时前向散射方向是z后向是-z。COMSOL三维远场方向图默认用的是球坐标极角和方位角建议先画二维切片方向图确认前后向趋势再做三维全空间图。7. 我在这类仿真中沉淀的几个习惯最后分享几个我个人的操作习惯。第一本机跑三维频域扫描时内存不够先从对称性下手。前面说的PEC/PMC半模型或四分之一模型可以在不损失精度的情况下把自由度降到原来的四分之一以下时间也快很多。第二多极分解脚本固定下来后保留小球自检用例以后换版本、换电脑、换材料库之后都先跑一遍自检。第三所有材料数据单独放一个插值文件不要散落在材料节点里这样换波长范围、换背景介质时改起来很快。最近我在尝试把同一套多极分解流程用到带衬底的圆柱和周期性阵列上结果会牵扯到衬底反射、阵列耦合之类的额外效应多极子之间的正交性也会被破坏需要做更多验证。这也是纳米圆柱多极散射分析下一步比较有挑战的方向。如果你也在这条路上调试希望上面这些细节能帮你少走几步弯路。