把文献里那套THz超构表面的BIC结构用COMSOL完整复现一遍是件看着容易、做起来容易卡壳的事。我第一版照着论文尺寸画完模型、扫完透射谱曲线平得像条直线BIC连影子都没有后来改成先做本征模式分析、再回头找透射特征几分钟就定位了问题。这篇文章就把整套思路、参数配置和踩坑过程写清楚给已经会COMSOL基本建模、但想复现光学/太赫兹超构表面文献里“高Q共振、准BIC”现象的同学做参考。1. 动手之前先把BIC超构表面这个模型拆明白1.1 先搞清楚BIC现象在仿真里长什么样连续谱束缚态Bound States in the ContinuumBIC这个名字听起来很绕但物理图像其实不难。常见的“困住光”的办法是法布里-珀罗腔靠两端反射镜把光来回反射Q值由镜子反射率决定。BIC是另一条完全不同的路结构是开放周期阵列光却能一直待在结构里不向外辐射不是因为有什么墙挡住它而是这个模式的对称性跟自由空间辐射通道“正交”电磁波在远场上相消干涉向外辐射的通道等于被关闭了理论上Q趋于无穷大。放在COMSOL里这个特性能直接体现在本征频率上我们算出来的特征频率通常写成 ( f \text{Re}(f) i,\text{Im}(f) ) 的形式虚部代表模式损耗。普通泄漏模式虚部很大BIC模式的虚部则趋近于零。所以复现BIC最直接的目标就是找到那个虚部极小的本征模式再通过对称性破缺把它变成能观测到的准BICquasi-BIC共振。THz波段在这个领域很常用因为超构表面的结构特征是百微米量级加工相对容易而且太赫兹源和探测器现在也比较普及实验上可以测到透射谱中的尖锐Fano峰。仿真上要复现的无非是这几件事本征模式频率、Q因子、透射谱线型以及Q值随结构不对称度的变化趋势。1.2 复现到什么程度算“复现成功”我给自己定了三条验收标准也建议你这样做不然很容易陷入“算完不知道对不对”的状态第一本征频率与文献报道的共振频率对得上一般偏差不超过几个百分点偏差主要来源是材料折射率的取值差异。第二Q因子的数量级和变化趋势对得上尤其是准BIC的Q值与破缺参数之间的幂律关系这个比单个点更可靠。第三透射谱上能看到经典的Fano不对称线型共振位置跟本征频率匹配线宽对应的Q值与从本征模式提取的Q值在同一量级。注意一点如果结构完全对称BIC点你在透射谱上是看不到任何共振特征的因为模式压根不往外辐射远场观测不到。必须引入一个对称性破缺参数把它变成准BIC才可能在透射谱上看到一个窄的Fano线型。这点新手特别容易踩坑我后面还会再强调。1.3 为什么选COMSOL来做这件事做超构表面仿真备选工具其实不少CST、Lumerical FDTD、COMSOL等。FDTD在时域里处理宽频透射谱很方便但遇到高Q值时要求极长的仿真时间因为脉冲能量在结构里迟迟不衰减CST在微波射频领域很顺手但做参数化扫描和模式分析时不如COMSOL直观。COMSOL的优势在于多物理场耦合和“本征模式频域扫描”双轨制。尤其算高Q模式时直接用频域扫频会疯狂吃网格和时间而本征模式求解器算尖锐共振是天然优势——它不依赖扫频分辨率直接解出复频率速度非常快。另外COMSOL的参数化扫描、辅助扫描和结果后处理流程很成熟适合做Q值随几何参数变化的趋势研究。所以我最终选择了COMSOL。2. 几何与材料搭出THz超构表面的最小晶胞2.1 单元结构选型和参数参考超构表面的BIC方案很多常见的有全介质椭圆柱阵列、圆柱阵列、开口环、以及介质板上打孔等。我这里用一个全介质硅椭圆柱阵列作为示例这也是太赫兹波段文献里非常经典的一种配置每个周期单元是一个高阻硅椭圆柱周期固定椭圆有一个旋转角参数。当椭圆长轴和短轴分别平行于晶格x/y方向时某些模式受对称性保护呈现BIC形态把椭圆整体旋转一个角度θ对称性降低BIC就变成准BICQ值随θ增大而降低。一些参考尺寸如下这些是依据常见文献量级取值的你也可以按自己手上的论文调整周期 (P_x P_y 200,\mu\text{m})椭圆长半轴 (a 50,\mu\text{m})椭圆短半轴 (b 40,\mu\text{m})柱高 (h 100,\mu\text{m})椭圆旋转角度 θ从 0° 扫到 10°这个尺寸下结构的米氏共振大约落在 0.8 THz 附近具体频率和材料折射率、柱高关系密切。圆柱或双圆孔的方案也可以照葫芦画瓢操作流程完全一样只是几何草图绘制时换成对应的形状。2.2 在COMSOL里搭建模型的基本操作顺序版本上我用的是COMSOL 6.x界面和5.x差别不大。流程如下新建模型时选择三维模型向导搜索“电磁波频域”emw这个物理接口就够用了。研究类型先选“特征频率”因为第一步要定位本征模式等你把模式摸清楚了再在同一模型中添加“频域”研究算透射谱。几何单位建议直接设成微米这样画草图时不用每次换算单位。COMSOL默认内置单位是米但几何里设置“单位”为微米后所有尺寸标注都会按微米显示求解时自动转成国际单位不影响结果。草图画法先在工作平面里画一个 200 μm × 200 μm 的矩形作为空气域底面再画一个以原点为中心的椭圆长半轴50、短半轴40椭圆旋转角度由全局参数θ控制。拉伸操作把两个域同时向上拉伸其中硅柱拉伸100 μm空气域建议上表面留出至少200 μm的空气层也就是空气域总高度至少300 μm后续算本征模式时还要加PML层所以空气层高度宁多勿少。材料分配上别偷懒。空气用内置Air材料相对介电常数1电导率0高阻硅需要手动配一个材料相对介电常数按折射率平方来设后面详细说。在COMSOL里材料节点里的域选择要确认选中的是圆柱体而不是空气域选反了共振频率会完全对不上。2.3 材料参数和单位陷阱高阻硅在太赫兹波段折射率大约在3.41到3.45之间取决于具体电阻率和工艺。很多文献直接取 (n3.45)也就是相对介电常数 (\varepsilon_r \approx 11.9)。仿真初期我建议先按无损材料算也就是设实部为11.9虚部为零。这样得到的Q值纯粹是辐射Q方便跟理论趋势对比。这里有个单位上的坑COMSOL的电磁波频域接口默认频率单位是Hz你写特征频率搜索值时如果想用THz直接写“0.8[THz]”是可以的COMSOL认识这个单位如果不放心就统一写 8e11 [Hz]永远不出错。更关键的是折射率与介电常数的换算千万不能用 (n3.45) 就直接填 (\varepsilon_r3.45)那样共振频率会偏得离谱。正确的是 (\varepsilon_r n^2)。3. 本征模式计算直接找到BIC并提取Q因子3.1 为什么要先算本征模式而不是直接扫透射谱这话我建议你刻在屏幕边上高Q结构不要一上来就做频域扫频。原因很现实THz超构表面单元尺寸是微米级但波长在亚毫米级整个模型如果要包含足够厚的空气层和PML自由度轻松几十万上百万。如果扫频范围要覆盖整个感兴趣波段步长还得小于共振线宽高Q共振线宽可能只有零点几GHz计算量会非常吓人。本征模式求解器则不同。它不关心你扫多少频点而是直接在一个搜索区间里解出若干个模式的本征频率和场分布一次求解就能定位所有候选模式。对于BIC这种虚部极小的模式本征模式计算的效率优势是碾压级的。我个人的习惯是先做本征模式把结构的所有模式都摸清楚尤其找到虚部最小的那个模式然后用它的实部频率作为频域扫描的中心只做窄带扫频。3.2 边界条件、周期条件和PML设置本征模式计算时模型四个侧面要设置周期条件。在COMSOL的物理场树里找到“周期条件”节点分别对x方向两个面、y方向两个面创建“周期”边界。周期条件类型选“Floquet周期”k矢量设为 ((0, 0, 0))对应第一布里渊区的Γ点。保护型BIC一般就在Γ点出现所以k矢量归零没问题。顶部和底部的开口边界需要慎重。最简单的做法是散射边界条件但高Q模式对外界数值反射非常敏感散射边界条件不一定能干净地吸收泄漏场虚部可能被抬高Q值偏低。更稳妥的是在顶部和底部再画一层空气和PML层。PML层材料很特殊COMSOL里直接在几何画好域后在“电磁波频域”物理接口下面添加“完美匹配层”特征选中对应的PML域即可。PML厚度建议至少半波长以上太薄吸收不干净也会污染高Q模式的虚部。算1 THz附近的模式时空气中波长300 μmPML厚度给200 μm起步比较踏实。3.3 从特征频率结果里把BIC选出来研究设置中特征频率搜索是一个很关键的地方。我把“所需模式数”设为20左右“搜索基准频率”设在0.8 THz附近。这个基准值怎么定可以先跑一次宽松搜索比如模式数设成10基准频率0.5 THz看看在0.5到1.2 THz区间的模式分布然后调整基准重新搜索。求解完成后结果里会给出每个特征频率的实部和虚部。在“全局”表格里能看到这样的数据实部单位Hz虚部也是Hz。如果某个模式的虚部是负值说明模式随时间衰减负得越少Q越高。Q值计算公式是[ Q \frac{\mathrm{Re}(f)}{2,\mathrm{Im}(f)} ]这里的 Im(f) 要取绝对值。比如实部 0.85e12 Hz虚部 -4e7 HzQ就是 0.85e12 / (2×4e7) ≈ 10625。但光看数值还不够一定要看场图确认模式形态。在结果节点里新建三维绘图组画电场模或磁场模分布。特征频率模式下生成的默认图通常会显示某个模式的场在“数据集”下拉框里可以切换不同模式。BIC模式的场应该集中在硅柱内部或者柱体周围很小的区域向外辐射的传播场分量几乎看不见普通辐射模式的场则会向空气层深处延伸有明显的行波特征。3.4 网格收敛高Q模式最容易翻车的地方高Q模式对网格质量极其敏感这句话怎么强调都不为过。频率和场分布算准了但Q值可能因为网格精度差了一个量级这在做网格无关性验证时经常遇到。一套比较稳的操作是先给整个域设置最大单元尺寸为 (\lambda_{\text{air}}/8)这里 (\lambda_{\text{air}}) 是搜索频率对应的空气波长。比如0.8 THz对应波长375 μm那空气域最大单元约47 μm。硅柱内部波长更短 (\lambda_{\text{eff}} 375/3.45 \approx 109,\mu\text{m})最大单元按 (\lambda_{\text{eff}}/10) 来加密到10 μm左右。然后做两套网格一套标准网格一套整体加密一倍对比Q值变化。如果Q值变化超过10%基本判定网格还没收敛需要继续加密。这里有个实践技巧Q值对结构内部场的网格更敏感而不仅仅是全局加密。所以可以给硅柱一个单独的“更精细”网格尺寸同时保持空气域网格中等这样Q值能收敛得快一些内存也不会爆炸。4. 对称性破缺把BIC变成准BIC并控制Q值4.1 为什么BIC对对称性这么敏感我这里说的对称性破缺指的就是破坏模式与辐射通道之间的正交条件。一个理想的BIC模式在远场的辐射分量为零本质原因是结构具有某种对称性比如沿x轴、y轴的镜面对称模式在这些对称操作下的变换方式与平面波辐射通道不同彼此不耦合。一旦你把椭圆旋转一个角度或者两个圆柱的半径变得不一样结构对称性降低原来被“禁掉”的辐射通道重新打开模式开始漏光但漏得不多——只要破缺程度小Q值依然很高这就是准BIC。Q值与破缺参数之间通常满足一次幂平方反比关系(Q \propto 1/\gamma^2)其中γ是表征不对称度的参数。这个关系是后面跟文献对比的核心也是一篇文章数据好不好看的关键。所以在COMSOL里我们要做的就是定义好这个破缺参数例如椭圆旋转角θ然后扫描θ从0到几度观察Q值如何下降。4.2 引入不对称参数的两种常见做法第一种是我示例里用的参数化旋转。全局参数θ控制椭圆的旋转角度θ0°时椭圆长短轴分别平行于x/y轴θ从0.5°、1°、2°一直扫到10°。几何上椭圆旋转角度可以通过“旋转”特征来设置角度直接引用全局参数θ。第二种做法是双圆柱或双椭圆结构。两个同尺寸圆柱对称地放在晶胞内当它们完全相同且位置关于中心对称时体系存在保护型BIC拉开直径差或位置偏移后BIC解除。这种方法在文献里也很常见参数化时控制的是半径比或者偏移量。我这篇文章里以旋转角θ为例其他结构殊途同归。4.3 参数化扫描与模式追踪在“研究”节点下展开“扫描”添加全局参数θ扫描列表里填一组角度0°、0.5°、1°、2°、4°、6°、8°、10°。这样一次求解就能得到所有角度的本征模式。但问题来了扫描θ时COMSOL会为每个θ独立求解不会自动告诉你“这个模式是上一个角度那个模式的延续”。从BIC到准BIC的过程中模式频率变化不大但计算时可能混入其他杂散模式。怎么追踪我的土办法比较有效每个θ计算完后把特征频率按实部排序挑出频率和场形态最接近上一组的模式。更自动化的方式是“辅助扫描”配合“特征频率搜索方式”设为“靠近上一角度解附近”这样求解器会把上一次的解当作初始猜测不容易跳模。但这招在COMSOL老版本上支持得不太好保险起见还是结合场图人工确认。具体操作时我会把单一角度θ0.5°的模式场图与θ0°的BIC场图放一起比较确认局域场形态一致只是多了一点点向外辐射的小尾巴说明确实是同一个模式。4.4 Q值随破缺参数的变化规律把每个θ对应的Q值提取出来在Origin或Python里画Q-θ曲线。理想情况下双对数坐标里是一条斜率为-2的直线。我举个例子说明趋势具体数值会随结构变化θ度相对Q值量级仅示意0Q趋向无穷数值上受噪声限制1约50002约12504约3008约80如果曲线斜率偏离-2先检查模式是不是追踪错了再检查网格收敛性。如果θ特别小的时候Q值有明显饱和趋势可能已经碰到数值极限或者材料损耗被计入后起了决定性作用。无损仿真里加的PML不理想也会导致Q饱和。另外θ0°时完全BIC的Q值理论上无穷大但数值上因为有限精度、网格离散化、PML不完美虚部不会精确为零。这个“数值Q”可能高达10^7到10^9但别拿它当成真实物理量来报它只说明这里确实接近BIC条件。5. 频域透射谱用可观测的方式验证BIC5.1 模型改造Port、空气层与PML本征模式阶段用的模型在算频域透射谱时要做一些改造。首先为了保证入射波干净地进入结构并透过顶部和底部要足够厚的空气层。我的习惯是在硅柱上下各留至少一个波长的空气也就是300 μm以上其次空气层最外侧不再用简单的散射边界条件而是设置Port边界。COMSOL中端口要设置成周期型端口。具体操作是添加“端口”特征类型选“周期性端口”分别给顶部和底部两个面。顶部Port 1作为入射端口模式类型选平面波指定电场极化方向沿x或y轴底部Port 2是出射端口对应透射波的接收。PML层在这个模型里还保留吗如果在Port外面再套PML概念上会冲突因为Port本身就处理了向外辐射的通道。我的做法是做透射谱时不再额外加PML顶部/底部直接落在Port边界上做本征模式时则不加Port、加PML两者各司其职。5.2 扫频设置从粗扫到精扫高频Q共振的第一个陷阱就是扫频步长。假设你算出的准BIC在0.85 THzQ5000那么共振线宽大约是 (f_0/Q 850,\text{GHz}/5000 \approx 0.17,\text{GHz})。如果你整个扫描从0.4 THz扫到1.2 THz步长取0.005 THz5 GHz那可能一个共振峰只落一两个频点曲线看起来就是平的完全看不到共振。正确的做法是两段式扫描先粗扫一次范围可以宽一点比如0.6到1.0 THz步长5 GHz目的是看透射谱大概在哪个位置有凹陷或尖峰然后锁定峰值区域做精细扫描。精细扫描范围我通常取共振中心±10线宽比如0.849到0.851 THz步长设为线宽的十分之一也就是0.017 GHz左右这样一条光滑的Fano曲线就能出来了。注意频率单位很绕写参数时一定要确认单位。COMSOL里频域扫描的“频率”范围可以用 [Hz] 写也可以用 8.5e11[Hz] 这种形式。建议直接用 GHz 写比如扫频范围写“849[GHz] 到 851[GHz]”COMSOL支持这种单位表达式直观不易错。5.3 Fano线型与洛伦兹拟合透射谱上的BIC共振通常呈Fano线型不是简单对称的峰或谷。透射率表达式可以用标准的Fano公式拟合[ T(\omega) T_0 A,\frac{\left(q\frac{\Gamma}{2}\omega-\omega_0\right)^2}{(\omega-\omega_0)^2\left(\frac{\Gamma}{2}\right)^2} ]这里的 (\omega_0) 是共振角频率(\Gamma) 是线宽能量维度q是Fano参数控制线型不对称程度。拟合得到 (\Gamma) 后Q值用 (Q\omega_0/\Gamma) 计算。如果你不想做Fano拟合还有一个粗略做法直接读取透射谱峰或谷的半高全宽FWHM用[ Q \approx \frac{f_0}{FWHM} ]算出一个数量级估计。但Fano线型本身不对称直接量FWHM误差不小只能用来快速核对正式对比文献时还是建议拟合。COMSOL里导出数据的方法很简单在“一维绘图组”里用“全局”绘制 S21 参数然后右键导出“数据”格式选CSV。拿到数据后用Python的scipy.optimize.curve_fit或者Origin的自定义函数拟合都行。拟合参数初值用手动估计(\omega_0) 设为透射谱上最陡峭位置附近的频率(\Gamma) 先给 0.1 GHz多试几次就能收敛。5.4 与文献结果对比时怎么调整复现文献时最常遇到的就是共振频率对不上。碰到这种情况我一般按照嫌疑顺序排查第一材料折射率。文献如果用的是3.41而不是3.45共振频率可能偏移1%到2%。第二结构尺寸。椭圆的半径、高度、周期差几微米对THz波段影响可能非常显著一定要逐项核对论文的方法段落和数据补充说明。第三衬底影响。我前面建议一开始不做衬底但很多文献的样品是放在衬底上的衬底会显著改变共振位置和辐射特性。如果需要带衬底复现几何里必须加上介质层而且物理机制会更复杂。Q值对不上的时候优先检查材料损耗。THz高阻硅损耗很小但文献中如果特别强调实测Q值远低于仿真Q值那很可能就是材料吸收和工艺粗糙度贡献的仿真无损模型本来就该高于实测。6. 实操中遇到的坑与排查方法6.1 透射谱上找不到共振最常见的原因是什么这个我遇到过太多次了。透射谱扫完一片平滑没有峰也没有谷第一反应不要怀疑物理接口先看对称性破缺参数是不是0。如果θ0°这正好处在BIC点共振模式不辐射透射谱上根本不会有任何特征。这时候把椭圆旋转角度设成2°或者更大重新算透射谱尖峰立刻跑出来。第二个常见原因是入射偏振方向跟模式不耦合。BIC或者准BIC模式有明确的对称性和偏振选择性如果入射电场沿x方向而你要观测的模式沿y方向偏振那就完全激发不出来。把Port里的极化方向改一下看另一条偏振通道的结果。第三才是扫频步长问题。如果破缺参数和偏振都调了还看不到峰就用本征模式算出来的频率做中心只做窄带精细扫描别省那几步。6.2 特征频率搜不到想要的模式该怎么调整特征频率搜索最怕的就是“搜索基准频率”给得不准导致所需模式数够不到目标模式。我一般是先扩大模式数比如一次算40个扫描区间从0.3 THz到1.5 THz看看这些模式里虚部最小的几个分别落在哪个频段。这样虽然算得慢但能对结构的所有共振有个全貌。搜到一批模式后按虚部绝对值排序虚部最小的那些就是高Q候选。如果发现所有模式的虚部都大得离谱先从网格和PML找原因。网格太粗会抬高所有辐射模式的虚部PML厚度不够也会导致泄漏场反射回来破坏模式纯度。有时候还会遇到完全没有虚部很小的模式那说明当前几何参数下就不存在近BIC。可以扫描一下椭圆长短轴的比值或者微调柱高让结构进入能支撑高Q的相区。不要在一个参数配置上死磕。6.3 模式跳变与Q值虚高虚低参数扫描时θ从0°变到8°模式追踪最大的坑就是“跳模”。你可能以为自己在追踪同一个模式实际上某个θ点之后求解器跳到了另一个模式上Q值曲线突然不连续或者频率突变。排查方式还是看场图把每个θ下我关心的模式场图截出来确认局部场分布一致。如果某一步场形态变了那就是跳模了需要把那个点的搜索基准频率手动改成上一个θ计算出的频率重新求解。Q值虚高或虚低的问题则更隐蔽。虚高一般出现在θ极小的情况下数值噪声和PML不完全吸收会让本应无限的Q变成一个样本很大、但不稳定的大数。虚低往往来自网格不收敛或材料损耗。所以报告Q值时我建议至少说明网格收敛情况和是否包含材料损耗这样数字才是公平的。6.4 内存和求时长不够用时我怎么做THz超构表面模型并不算太大但加上PML和精细网格后特征频率求解也可能吃掉几十GB内存。如果你机器不够富裕试试下面几个办法第一个办法是分两步走。先用较粗网格和较少模式数定位目标模式的频率和基本性质再用细网格单独算目标模式。网格加密后求解器只针对一个频段搜索内存和时间都省。第二个办法是减少PML层的厚度和网格密度PML域里的网格不需要跟结构域一样细用映射网格或稍粗的自由网格就行。第三个办法是优先使用直接求解器比如MUMPS或PARDISO高Q模式对迭代求解器的数值误差更敏感迭代解法器容易假收敛。如果你做频域精细扫描时觉得每个频点都重新求解太慢可以试试把扫频范围压得再窄一点。我见过有人从0.1到1.5 THz全线扫步长还取100 MHz机器跑了两天没跑完。正确流程是先本征定位、再窄带精扫把资源集中在该用的地方。最后说一点我自己的习惯这套仿真做完之后我最大的体会是复现BIC的工作量其实不大难的是每一步都清楚自己算出来的是什么。比如看到虚部很小的特征频率要知道那是BIC看到透射谱上的Fano峰要能把它跟某个高Q本征模式对应起来看到Q值随θ的变化曲线要能判断它是不是符合平方反比。把这些环节串起来COMSOL在你手里就不再是“照着文献点按钮”的工具而是可以自由改参数、设计新结构的实验台。还有一个小技巧也顺便分享每次跑完参数扫描我都会把特征频率实部、虚部和Q值一起导出到表格连同对应θ存成一个CSV文件。后面写论文或者做汇报时这些中间数据比最终那张漂亮的透射谱图有用得多。你复现时也可以这么干后面省很多事。