企业官网建设流程全解析

第一次在COMSOL的特征频率列表里看到某个模式的虚部小到10的-7次方量级时，我的第一反应是求解器又在拿伪模式糊弄我。把网格加密三遍、周期性边界逐项核对、换了两台机器跑了一整夜之后，才确认那确实是光子晶体超表面里的连续束缚态（BIC）——一个Q值超过千万的暗模式。这类模式在能带图上表现为“离散能级嵌进连续谱却完全不辐射”，做超表面的人听到都会多看两眼：理论上Q可以无限大，意味着无限窄的共振线宽、超强局域场，是做滤波器、激光器、传感器的黄金结构。但BIC有一个绕不过去的毛病——偏振选择性太强，常见对称保护BIC只对特定偏振的入射光有效，换个偏振方向共振就没了。这篇文章记录的是我最近一次基于介质光子晶体超表面做极化无关BICs多极子分析与COMSOL仿真模拟的完整过程，内容包括双偏振BIC结构设计逻辑、COMSOL周期性超表面建模仿真细节、从近场数据做多极子展开定位辐射消失物理源头的方法，以及网格、求解器、后处理阶段一批实打实的踩坑记录。全文围绕“极化无关”和“多极子分析”两个核心词展开，适合正在做超表面仿真、想复现BIC现象的同行，也适合刚接触COMSOL特征频率计算的人。

1. 从“零线宽”说起：BICs的物理本质与极化无关思路

1.1 连续束缚态不是直觉里的“束缚”

连续束缚态最早是量子力学里的反常概念：在连续谱能带内部，居然存在一个不衰减的束缚态。放在光子晶体语境下，这个角色由本征模式来承担——模式频率落在辐射连续谱范围内，却因为没有通道把能量带到无穷远而“卡”在结构中。等价描述是远场辐射完全被禁止，共振线宽的极限为零，Q因子发散。

这个“禁戒”不是偶然的。光子晶体超表面中的BIC，从对称性角度大致分两类：一类是对称保护BIC，靠结构对称性硬生生堵住所有辐射通道；另一类是参数型BIC，靠结构参数调谐使不同辐射通道在远场干涉相消。前者更直观，后者更隐蔽，但也更有意思。

我这次仿真选择的是正方形晶格介质圆柱阵列，这是最常见、最容易复现的BIC承载结构。介质圆柱本身在Γ点附近就能支持对称保护BIC，而且圆柱结构保留了C4v旋转对称性，为后来做极化无关提供了非常自然的延伸空间。介质选硅，折射率3.48，与空气背景形成高折射率对比，能带间隙更大，BIC局域性更强，仿真和实验都容易对齐。

1.2 多极子分析在BIC问题里扮演什么角色

远场辐射本质上是多极子辐射的叠加。从源的角度看，单胞内部感应出来的电流/极化可以分解为电偶极（ED）、磁偶极（MD）、电四极（EQ）、磁四极（MQ）等分量的贡献，远场就是这些分量相干的叠加。BIC出现的物理路径只有两条：要么所有辐射分量同时消失，即模式在结构对称操作下具有特定的奇偶性，任何常规辐射通道都无法匹配这种对称性；要么源分量并不为零，但不同分量的辐射在远场干涉相消。

光看近场场图，你只能得出“这模式挺好看”的直观印象，却判断不出它为什么不辐射，更回答不了“把几何参数改一点，它还能不能保持不辐射”这个工程上更关心的问题。多极子分析的价值就在这里——把近场数据展开到多极子基上，把每个分量的幅值和相位都算出来，才能解释辐射消失的机制，也才能在参数空间里准确预测下一个BIC出现在哪里。

我在COMSOL里做这件事的思路，不是用复杂的第三方工具箱，而是直接从求解出的电场数据积分出偶极矩、四极矩。这个方法的好处是透明、可控，每一步物理含义都清楚，调参数时也能非常直观地看到是哪个多极子分量在变化。具体怎么做，后面第3部分详细拆解。

1.3 “极化无关”的双偏振设计逻辑

单根圆柱在Γ点得到的对称保护BIC，通常只对应一种特定的场偏振；换成正交偏振激励，耦合路径跟着改变，共振条件就没了。极化无关设计的目标，是让同一周期结构对x偏振和y偏振的入射波都提供近似等同的BIC特征。

工程上有三条常见路线：

• 第一条是双原子基元法，单胞里放两个形状相似但朝向相差90度的散射体，分别支起TE-like和TM-like模式，各带一套对称保护BIC，再通过子晶格耦合把两支模式的频率凑近。
• 第二条是几何微扰法，在保持C4v旋转对称性的圆柱阵列基础上微调几何，让模式对x与y偏振的响应等价。最省事的情况就是圆柱阵列本身——在Γ点同时存在TE-like和TM-like两支模式，它们在对称性上各自禁戒辐射，天然具备极化无关潜力，关键是调整高度和半径让两支模式的工作频率重合或接近。
• 第三条是引入各向异性介质，利用双折射材料的折射率椭球分别设计TE/TM两个通道的k空间模式分布。

我最后走的是第二条路线，因为结构最简单、参数最少、可复现度最高，而且与多极子分析的结合最顺。如果你要对比不同路线，可以看下面这个表。

2. COMSOL光子晶体超表面建模仿真的关键设置

2.1 几何与材料：模型不用复杂，但边界要用对

很多新手问“光子晶体超表面COMSOL仿真怎么建”，其实是卡在边界条件的逻辑上。模型本身很简单：一个正方形单胞，中间一根硅圆柱，周围是空气。具体参数我这次用的是晶格常数a=700 nm，圆柱半径r=150 nm，高度h=600 nm。硅折射率设为3.48，损耗先关闭——研究理想BIC时千万别开材料损耗，否则Q因子会被材料吸收压到几千甚至几百，完全看不出趋势。

COMSOL里具体操作步骤：

1. 组件几何：在二维工作平面画一个a乘a的正方形域，中间画一个圆形域，向下拉伸成硅圆柱，再向上、向下延伸空气层作为包围域。空气层的厚度至少要留出半个波长以上，后面才不会影响近场积分。
2. 物理场选择“电磁波，频域（ewfd）”。材料节点里设置硅的折射率或相对介电常数，空气保持默认。
3. 单胞的x、y边界设置为Floquet周期性条件，把布洛赫波矢的分量指定为全局参数kx和ky，注意单位一致。
4. 如果只算无限薄板的模式，z方向两端可以用周期性条件封死；但要做远场分析和Q值衰减特性，建议在z方向留出空气层后加完美匹配层（PML），避免边界反射制造的假模式。

2.2 特征频率研究与k空间扫描

BIC一定出现在k空间里某个孤立点，最常见的就是Γ点附近。研究类型选择“特征频率”，在全局参数里定义kx和ky，然后加上“参数化扫描”，扫从0到π/a方向的若干k点。

这里有个关键经验：COMSOL特征频率求解器默认输出模式数量有限，如果不设置基准频率，它经常只给最低阶的几个模式，而你要找的BIC可能在几百太赫兹的高频区。所以第一步先做一次粗扫描，比如搜索范围设在100 THz到400 THz，一次求出10个模式，用电场模分布区分哪些是硅圆柱内的驻波模式、哪些是空气腔里的杂散模式。锁定目标模式后，把搜索基准收窄到该模式附近，再精细扫k，数量设为2即可。

k空间扫出来的结果，BIC在数值上并不会严格给出“虚部为零”，因为数值离散存在误差，Q值只会随网格细化不断提高。判断BIC是否存在的标准，是虚部的绝对值能否随网格加密而持续下降，而不是直接等于零。

2.3 6.x版本的操作体验与批处理

我一直习惯在Linux服务器上跑大批参数扫描，模型量一多就受不了手动点界面。COMSOL 6.4对特征值求解器和并行效率的提升比较明显，同一套模型，多核并行时的内存占用比老版本稳定很多。如果你是跑参数扫描找BIC，建议直接上批处理，用COMSOL with MATLAB接口循环改参数、跑研究、存结果。比如下面这种循环：

model = mphopen('bic_supercell.mph'); kxScan = linspace(0, 0.02, 21); results = zeros(length(kxScan), 4); for i = 1:length(kxScan) model.param.set('kx', kxScan(i)); model.study('std1').run(); freq = mphglobal(model, 'freq'); Q = abs(real(freq(1))/(2*imag(freq(1)))); results(i,:) = [kxScan(i), real(freq(1)), imag(freq(1)), Q]; end

这段示意代码里的freq变量名要按模型里实际定义的全局表达式来取，但循环结构是通用的。如果你的运行环境没有MATLAB，COMSOL的Java API也能做类似的事，社区里还有封装好的脚本框架，可以用Python之类的外部语言间接控制，不过版本兼容要格外小心。

3. 多极子分析的完整实现路径

3.1 从复电场到多极矩：积分公式与坐标陷阱

多极子分析不需要额外买模块，靠COMSOL的派生值计算就能完成。核心思路是：先得到单胞内的复电场分布，然后构造极化电流密度，再按多极子定义做体积分。

假设时谐因子为e^{jωt}，单胞内极化电流密度可以写成：

• 极化电流J = -jω(ε_r - 1)ε0 E
• 电偶极矩p = ∫ (ε_r - 1)ε0 E dV
• 磁偶极矩m = ½∫ r × J dV
• 电四极矩张量Q_ij = (j/(2ω)) ∫ [x_i J_j + x_j J_i - (2/3)δ_ij x_k J_k] dV

注意，这些公式里的坐标原点必须放在单胞几何中心。我看过不少初学者直接把坐标原点留在CAD导入时的角点上，结果高阶多极矩全部被平移带来的虚假分量污染，偶极矩大得离谱，整个分析失去意义。每次改模型后，都要检查一下坐标系原点是否还在单胞中心。

还有一个细节：如果研究的是磁谐振超表面，材料参数里需要加入磁导率，多极子展开还要补上磁流密度项。光学频段绝大多数介质超表面都可以忽略磁响应，所以上述公式够用。

3.2 近场积分提取多极子系数：COMSOL具体操作

在COMSOL中要做的是：

1. 在“定义”节点下创建一个“积分”算子，作用域选整个单胞域。
2. 在“派生值”的“全局计算”里，定义如下表达式：intop((epsilon_r_const-1)*epsilon_0_const*ewfd.Ex)，得到电偶极矩的x分量。
3. 类似地计算y分量、z分量；磁偶极矩表达式里用坐标分量x、y、z和电流密度分量ewfd.Jx等组合。
4. 电四极矩的表达需要逐项写出张量分量，虽然麻烦，但计算量很小。
5. 把这些全局表达式做成“表格”，每次求解完一键输出。

这套方法有两个好处：一是完全在COMSOL内闭环，不需要导出场数据，二是可以对每个特征频率和每个k点自动输出多极子分量，配合参数扫描非常方便。

3.3 远场拟合路径：什么时候该用

近场积分提取的是“源端”的多极子强度，但如果想精确分析远场干涉相消，特别是研究参数型BIC时不同通道的相位关系，需要走第二条路径——远场拟合。

在模型外边界（避免把PML包进去）激活“远场计算”，COMSOL会自动计算给定边界上的远场复振幅。得到远场数据后，用矢量球谐函数做最小二乘拟合，拟合系数就是各个多极子辐射系数。这个路径的优点是同时拿到幅值和相位，可以直接观察两支谐振通道之间的干涉；缺点是计算量比近场积分大，而且拟合阶数要取得适当，取少了精度不够，取多了反而引入数值噪声。

我的习惯是近场积分做主判据，远场拟合做辅助核对。对同一个BIC，若近场显示电偶极矩趋近于零、远场拟合也显示ED辐射系数为零，那分析就稳了。

3.4 多极子分量与对称性速查

这张速查表在做极化无关BIC分析时非常有用。比如你在Γ点看到一个TE-like模式，多极子展开结果只剩磁四极分量，说明它的核心机制是EQ/MQ高频通道的残余辐射被某种对称性压制；此时你只需要检查该模式在90度旋转操作下的对称性，就能预测它对不同入射偏振的响应。

4. 参数扫描寻找双偏振BIC：从能带到Q值奇点

4.1 搜索策略：在参数空间里找两条“消失线”

极化无关BIC的本质，是两支不同偏振的模式在同一个动量点上同时满足“零辐射”条件。实际操作中，这两支模式的频率通常并不重合，需要靠参数扫来拉近。我采取的策略是固定kx=ky=0，扫描圆柱半径r和高度h的二维参数空间，记录两支模式各自的频率和Q值。

BIC在参数空间里往往表现为Q值等高线上的一条脊——沿这条脊走，Q持续升高并向无穷发散。两支模式的BIC脊各有各的位置，如果它们在某一点交汇，那里就是双偏振BIC的候选点。

具体执行时，我用COMSOL的参数化扫描配合下面的流程：

1. 固定网格为全局较粗级别，扫描r和h，找到Q值超过1e4的候选区域。
2. 在候选区域附近加密网格，特别是圆柱顶部和底部界面的局部网格，重新计算Q值。
3. 对候选点做三次逐步加密网格的验证，观察log Q是否持续上升，如果出现平台说明并非真正BIC。
4. 微调两个参数，使两支模式频率差的绝对值最小化。

4.2 粗扫、细扫与局部细化：一个高效流程

第一遍扫描千万不要用很细的网格。BIC的Q值是出了名的对网格敏感，但频率实部没那么敏感，所以粗网格足以筛选候选区。粗扫的建议网格是每个波长至少8到10个单元。找到候选区后，第二遍做局部细化，把圆柱圆柱上下表面用边界层网格处理，边界层厚度取内密外疏的过渡，防止等距网格造成的空腔伪影。第三遍则是将整个单胞网格加密一倍，对比候选点Q值的增长趋势。

这套流程下来，单个候选点需要跑的次数在10到20次之间，但每次计算都很快，总体效率远高于一上来就全细网格扫全参数空间。

4.3 极化无关的判读方法

判读极化无关是否实现，有几种方法。

最直接的是：对x偏振和y偏振的入射波分别求透射谱或反射谱，比较两个偏振下共振峰的Q值。如果两条Q曲线在目标参数点附近接近重合，说明结构对偏振不敏感。COMSOL里可以用频域研究加两个极化方向的背靠背扫描实现。

另一方法是看远场偏振分布的涡旋奇点。BIC在远场表现为偏振椭圆度分布里的拓扑涡旋，旋涡中心就是辐射为零的动量位置。两支模式如果分别在两个动量位置出现涡旋，调参让两个涡旋中心重合，就可以判定极化无关。COMSOL的远场绘图支持直接看偏振各分量的相位分布。

我这次仿真的结果，在r=148 nm、h=610 nm附近，两支模式频率差只差0.7%，Q值同时超过10的6次方，算是实现了本文说的极化无关BIC候选点。这个结果是在无损耗配置下得到的，实际样品因为材料吸收，Q值会显著下降，但工作机制不变。

5. 常见问题与排查技巧实录

5.1 特征频率求解器只给伪模式

症状是扫描k点时频率曲线不连续，模式场型突然出现孤岛震荡，或者一堆模式挤在一起分不清。排查重点有三个：

• 网格在Floquet周期边界两侧是否完全对称。周期边界要求两侧贴合，网格不对称会造成非物理的散射，直接破坏模式对称性。
• 基频设置是否覆盖目标频率。搜索区间太宽时，求解器会输出大量无意义的高阶模式，建议基准频率设置在目标模式附近。
• 布洛赫参数的用法是否一致。周期性条件里填写的波矢，要和你扫描的实际物理k值保持一致，常见错误是归一化系数错了一个2π。

5.2 Q值虚部像过山车

如果同一个模式在一次计算里Q只有5000，改变网格后跳到10的6次方，不用奇怪——这是BIC仿真的常态。Q因子对网格的敏感度远高于频率实部。解决方法是给圆柱上下表面加边界层网格，并逐步加密做收敛性分析。

如果你加了PML，还要检查PML厚度。PML太薄会把泄漏模式的虚部拉大，造成Q值被低估。我的经验是PML厚度至少取中心频率波长的二分之一，太薄时结果不靠谱。

5.3 多极子积分数值异常

如果算出的电偶极矩明显偏大，先检查坐标系原点是否在单胞中心。坐标系原点偏移会造成高阶矩“串扰”到低阶矩。另一个常见问题是积分域只选了硅圆柱区，漏掉了空气域中泄露出来的极化电流。多极子积分应该覆盖整个单胞，包括空气区域——高频模式在介质边界外会形成渐逝场，这部分也有辐射贡献。

5.4 如何防止“假BIC”翻车

曾经有人在群里晒过一组Q值高达10的9次方的高Q结果，后来一查是“数值BIC”——网格误差造成的伪收敛，换个剖分网格就没有了。要避免这种现象，有一个很实用的验证手段：在结构上引入一个极小的对称破缺，比如把圆柱截面改成椭圆，长短轴相差1%。如果原来的BIC是真正的对称保护BIC，它的Q值会瞬间跌到1e3以下；如果跌幅很小，很可能只是准BIC或数值假象。参数型BIC对微扰的响应更复杂，这时要多极子分析结果和远场相位分布一起看。

下面把常见问题整理成速查表。

5.5 材料损耗对BIC的影响

理想BIC仿真理论上Q无限大，但任何实验材料都有损耗。硅在近红外波段吸收虽然小，但也不是零。做仿真时要记住两步走：第一步关掉损耗，找理想结构参数，看BIC机制是否成立；第二步打开材料损耗，再算一次Q，得到实际可实现的器件指标。很多刚接触BIC的人只跑理想情况，结果信心满满去加工样品，测出来Q只有几百，误以为设计失败。其实那只是因为材料损耗压了上限，机制本身没有问题。

我个人在实际操作中最大的体会是：BIC仿真的结果一定不要轻信默认网格。同一个模型，粗网格下Q可能只有3000，细化圆柱上表面网格后直接能到10的7次方，实部却几乎不变。先粗扫找候选区、再局部细化验证收敛，是省时间又不误判的关键。最后想提醒一句，如果你要把这套方法延伸到更复杂的结构，比如双层超表面或非周期微扰结构，多极子积分和远场拟合的代码是可以复用的，只是积分区域和对称性分析需要重新核对。先在小模型上把流程跑顺，再逐步扩大参数范围，这条路最稳。