企业官网建设流程全解析

把文献里那套THz超构表面的BIC结构，用COMSOL完整复现一遍，是件看着容易、做起来容易卡壳的事。我第一版照着论文尺寸画完模型、扫完透射谱，曲线平得像条直线，BIC连影子都没有；后来改成先做本征模式分析、再回头找透射特征，几分钟就定位了问题。这篇文章就把整套思路、参数配置和踩坑过程写清楚，给已经会COMSOL基本建模、但想复现光学/太赫兹超构表面文献里“高Q共振、准BIC”现象的同学做参考。

1. 动手之前，先把BIC超构表面这个模型拆明白

1.1 先搞清楚BIC现象在仿真里长什么样

连续谱束缚态（Bound States in the Continuum，BIC）这个名字听起来很绕，但物理图像其实不难。常见的“困住光”的办法是法布里-珀罗腔，靠两端反射镜把光来回反射，Q值由镜子反射率决定。BIC是另一条完全不同的路：结构是开放周期阵列，光却能一直待在结构里不向外辐射，不是因为有什么墙挡住它，而是这个模式的对称性跟自由空间辐射通道“正交”，电磁波在远场上相消干涉，向外辐射的通道等于被关闭了，理论上Q趋于无穷大。

放在COMSOL里，这个特性能直接体现在本征频率上：我们算出来的特征频率通常写成 ( f = \text{Re}(f) + i,\text{Im}(f) ) 的形式，虚部代表模式损耗。普通泄漏模式虚部很大，BIC模式的虚部则趋近于零。所以复现BIC，最直接的目标就是找到那个虚部极小的本征模式，再通过对称性破缺把它变成能观测到的准BIC（quasi-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就变成准BIC，Q值随θ增大而降低。

一些参考尺寸如下，这些是依据常见文献量级取值的，你也可以按自己手上的论文调整：

• 周期 (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之间，取决于具体电阻率和工艺。很多文献直接取 (n=3.45)，也就是相对介电常数 (\varepsilon_r \approx 11.9)。仿真初期我建议先按无损材料算，也就是设实部为11.9，虚部为零。这样得到的Q值纯粹是辐射Q，方便跟理论趋势对比。

这里有个单位上的坑：COMSOL的电磁波频域接口默认频率单位是Hz，你写特征频率搜索值时，如果想用THz，直接写“0.8[THz]”是可以的，COMSOL认识这个单位；如果不放心，就统一写 8e11 [Hz]，永远不出错。更关键的是折射率与介电常数的换算，千万不能用 (n=3.45) 就直接填 (\varepsilon_r=3.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 μm，PML厚度给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 Hz，Q就是 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的直线。我举个例子说明趋势（具体数值会随结构变化）：

如果曲线斜率偏离-2，先检查模式是不是追踪错了，再检查网格收敛性。如果θ特别小的时候Q值有明显饱和趋势，可能已经碰到数值极限，或者材料损耗被计入后起了决定性作用。无损仿真里加的PML不理想也会导致Q饱和。

另外，θ=0°时完全BIC的Q值理论上无穷大，但数值上因为有限精度、网格离散化、PML不完美，虚部不会精确为零。这个“数值Q”可能高达10^7到10^9，但别拿它当成真实物理量来报，它只说明这里确实接近BIC条件。

5. 频域透射谱：用可观测的方式验证BIC

5.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 THz，Q=5000，那么共振线宽大约是 (f_0/Q = 850,\text{GHz}/5000 \approx 0.17,\text{GHz})。如果你整个扫描从0.4 THz扫到1.2 THz，步长取0.005 THz（5 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文件。后面写论文或者做汇报时，这些中间数据比最终那张漂亮的透射谱图有用得多。你复现时也可以这么干，后面省很多事。