做超表面纳米光子学方向的人大概都绕不开连续谱束缚态BICs这个话题。我第一次在COMSOL里复现文献中的准BIC共振时最头疼的不是模型建不出来而是结构参数明明对齐了算出来的Q值却比论文里低了两个数量级——后来才明白问题出在我对“动量空间”这个自由度的理解不够。这篇内容我不打算按软件手册的顺序讲COMSOL基本操作而是从BICs在光子晶体超表面中的动量空间行为出发把动态调节的完整仿真流程拆开了讲覆盖单位晶胞建模、特征频率模式筛选、Q因子提取、带结构扫描和常见的数值坑。适合刚开始用COMSOL做超表面研究、已经能跑通频域透射谱但想进一步调BIC的研究生和工程师。1. 先搞明白BIC在动量空间里是什么状态动态调节调的又是哪个自由度1.1 对称性保护BIC与偶然BIC的本质区别连续谱束缚态这个名字听起来玄但它要表达的东西很直接一个模式的频率落在了连续辐射谱的范围里却不向外辐射能量。放在光子晶体超表面的语境里就是结构支持某个本征模式它的频率和周围空间能传播的平面波频率重叠但因为对称性或干涉相消的原因模式和所有辐射通道的耦合都归零能量就被“锁”在结构内不往外跑。文献里最常见的分类是两类。第一类叫对称性保护BIC通常出现在布里渊区高对称点比如Γ点。它的存在不依赖精细调节只要结构保持某种面内镜像对称、旋转对称模式远场的极性和可辐射平面波不匹配就能形成BIC。类比一下就像一个人站在旋转门正中间四个门扇对称分布他往哪个方向推门都推不开——对称性让所有通道都失效了。第二类叫偶然BIC也叫Friedrich-Wintgen型BIC它依赖两个模式之间的干涉相消需要在某个特定波矢、特定结构参数下才会出现调节精度要求更高。这两类BIC在动量空间里的表现很不一样。对称性保护BIC的位置基本钉死在Γ点参数变了它也不会挪窝偶然BIC则可能出现在Γ点以外结构参数一变它在动量空间里的位置还会漂移。做动态调节之前必须先把目标BIC归好类因为后面选参数、扫波矢、解释Q值变化趋势都要靠这张“身份卡”。1.2 动态调节的三种物理路径理解了BIC的本质再看“动态调节”就清晰了——BIC一旦受到扰动有辐射通道被打开就会变成准BICquasi-BICQ因子从一个理论上无穷大的数值跌到有限值。调节的手段大致有三种对应到COMSOL里就是三套不同的参数扫描流程。第一种是结构几何扰动最常用。比如把纳米圆盘改成椭圆或者在圆盘旁边加一个小孔、小柱用一个无量纲参数δ描述破缺程度。对称性保护BIC对几何扰动最敏感δ从0变到0.1的时候Q值可能从百万级降到百级跨了好几个数量级。第二种是入射角扫描也就是在动量空间里换个方向激发。超表面的能带结构里BIC存在于特定的平行波矢(kx,ky)处实验中通过改变斜入射角度可以让激发条件刚好落在BIC附近谱线从Fano线型逐渐收窄Q值在特定角度达到峰值。第三种是材料动态调谐比如改变周围介质的折射率、给结构加载电压调控载流子浓度能带整体移动BIC在动量空间的位置也会跟着漂。从仿真角度看第一种用参数化扫描结构尺寸第二种用Floquet周期条件里的波矢扫描第三种用材料折射率参数扫描。三者经常组合使用但核心逻辑都是一样的先找一个完好结构确认BIC存在再引入可控扰动观察Q因子和共振波长的变化规律。2. COMSOL建模从单位晶胞到可参数化光子晶体超表面2.1 几何、材料与周期的参数化设置COMSOL里做周期结构超表面标准路线是RF模块下的电磁波物理场研究类型选“特征频率”或“频域”。模型只要一个单位晶胞就够了水平方向用周期边界垂直方向加空气间隔和吸收边界这和实验里整个阵列对应起来理解晶胞无限重复等效于一个无限大的周期平面。拿一个最经典的例子做底子正方形晶格周期a 700 nm硅纳米盘半径r 180 nm高度h 200 nm背景是空气。硅的折射率取3.48在近红外波段色散不大可以先按常数算。COMSOL里参数表的设置我习惯这样组织参数名值说明a700[nm]晶格周期r180[nm]圆盘半径h200[nm]圆盘高度delta0几何扰动幅度扫描用kx0平行波矢x分量扫描用ky0平行波矢y分量扫描用n_Si3.48硅的折射率n_bg1.0背景折射率f0c/(0.7[um]*2.4)特征频率搜索基准值最后一项f0是我自己的习惯先估算一个大概的共振频率特征频率研究里把搜索范围放在它附近能省掉大量无关模式。估算思路很简单模式在硅和空气中混合分布等效折射率大概在2.2到2.6之间硅盘的有效波长大约是700nm除以2.4共振频率就在c/(0.7μm×2.4)附近算出来约等于1.78E14 Hz。实际跑的时候可以放宽一点左右各留20%的余量。把结构尺寸、周期、扰动参数全部建成全局参数后面做参数化扫描会非常省事不用每次改几何。2.2 Floquet周期条件和PML的放置逻辑模型里最关键的两个边界设置一个是水平方向的周期条件一个是垂直方向的开放边界。水平方向把单位晶胞的x和y方向的相对面分别设成周期条件Periodic Condition类型选Floquet周期性。COMSOL里会要求填k矢量也就是平行方向的Bloch波矢。这里要注意周期条件只能设一对面对x方向一对、y方向一对一共两对每个周期条件里都要输入kx和ky。当kxky0时对应正入射也就是Γ点想算带结构或者斜入射就得把这两个参数变成扫描变量。垂直方向也就是晶胞的上方和下方要考虑能量可以往外辐射。最稳妥的做法是加一层空气间隔再接完美匹配层PML。我给一个具体尺寸建议硅盘上下各加0.5λ的空气间隔PML厚度做成1λλ取模型中心频率对应的自由空间波长。PML本身在特征频率计算里会引入一批非物理模式能量集中在PML区域里后面选模式的时候要明确把它们剔掉。判别方法也简单画电场分布图看能量是不是一团一团地“糊”在PML层里——只要是这种模式一律不选。有些做法是用散射边界条件代替PML在纯BIC状态下问题不大因为模式本身不辐射但算准BIC时Q值往往在1E3到1E5之间散射边界对弱辐射的微小反射也会污染结果我实际对比下来还是PML更稳定。2.3 网格策略为什么BIC模拟对网格一致性要求极高网格可能是BIC仿真里最容易被低估的一环。BIC的Q值本质上是数值里虚部的倒数Q越高意味着模式往外漏的能量越小——小到和网格离散误差同一个量级。如果网格太粗数值损耗会“淹没”真实的辐射损耗算出来的Q光子晶体的Q值远低于理论值。反过来网格波动也会制造虚假的泄漏通道Q表现得比实际低。我的经验是硅盘内部网格最大尺寸取有效波长的1/8到1/10空气区域最多取λ/8。以700nm周期、硅折射率3.48为例硅内部光的有效波长约是201nm那硅盘内网格最大尺寸要控制在25nm左右空气区域控制在90nm左右。这个量级下的3D模型自由度大概在几十万到上百万一台16GB内存的工作站还能跑得动。但网格策略真正的核心是“一致”。做参数扫描的时候如果扰动参数一改网格也跟着自动加密Q值的变化里就混进了网格带来的变化干扰很大。我的做法是把几何扰动设计成不改变网格拓扑的形式比如圆盘变椭圆只靠一个参数控制半径网格重画但密度分布基本一致。如果结构里加小孔这种导致局部网格剧变的操作要专门对孔周边加固定尺寸约束保证孔在不同尺寸下网格密度近似。收敛性检查也别省把全局网格尺寸整体缩小一半选同一个模式看Q值变化偏差超过5%就说明网格还没收敛加密之后重算。3. 用特征频率研究锁定BIC模式复数特征值中读出Q因子的方法3.1 特征频率研究的设置与模式筛选模型建好后研究类型选“特征频率”。这一步的目的是直接求解模式的本征频率复数特征值里实部是共振频率虚部是辐射损耗的阻尼Q因子等于两者比值。求解器设置里有几个关键选项。特征频率搜索范围我给一个宽区间通常在估算的f0前后扫5E13 Hz每次求解10到20个模式。这个数量确实会带来大量无关模式但BIC必须找了才知道在哪里宁可多算不可漏算。物理场设置里电磁波特征频率研究默认解的是角频率单位是rad/s这一点对后面算Q尤其重要先确认设置里特征频率单位是角频率还是频率别算完发现差了6个数量级。算完之后后处理是对着一组特征频率列表逐个看电场分布。BIC模式的典型电场图是能量极强地局域在硅盘内部或者表面向外场衰减极快PML里几乎没有响应。对比之下普通辐射模式在PML区域还会有明显的振荡图样。我在最开始跑的时候犯过糊涂看到某个模式的电场分布很规整一度以为它是BIC后来一查Q值只有200再对比频域透射谱发现根本对应不上——直觉不能替代数值判据。3.2 Q因子计算公式与数值工具验证Q因子的计算式子很直接Q ω_real / (2|ω_imag|)如果COMSOL输出的是角频率直接代入实部和虚部单位会自动约掉如果是频率同样可以用f_real/(2|f_imag|)两种方式结果一致。虚部的符号一般是负的表示能量损耗取绝对值即可。举个例子假设某一个准BIC模式特征频率实部是1.414E15 rad/s虚部是7.07E9 rad/s那Q就是1.414E15除以1.414E10算出来约等于10万。这是一个可参考的量级完全对称结构下的BICQ值往往能到1E7甚至更高数字上的限制主要来自数值精度而加了几个百分比的扰动之后Q快速回落到1E3到1E5这是频域实验里能测得比较准的范围。特征频率法算的Q值不是一个结果就结束了我通常会在同一次参数扫描里配合频域研究做交叉验证。频域扫描里画透射谱谱线在共振峰附近呈现Fano线型共振峰的宽度Δλ对应Q λ0/Δλ。把两种方法算出来的Q放在一张双对数图里对比吻合的区间才能放心写进论文里。频域谱的好处是它直接对应实验可测的量特征频率法快但抽象两条路互相印证能避免特征频率模式下数值噪声造成的假结果。3.3 判断“这个模式是不是BIC”的三个判据初学者问得最多的问题就是跑完特征频率之后怎么从一长串模式里认出哪个是BIC。我给三个判据按重要性排序。第一Q值数量级。BIC对应的是模式虚部接近零Q在数值噪声允许的范围内趋于无穷。对称结构下Q值至少比周围普通辐射模式高两个数量级以上如果所有模式Q值都停在几百那大概率还没找到BIC或者网格不够细。第二动量空间的发散性。把kx设成0时高Q的那个模式往kx的正方向扫几个小步长观察Q的变化。BIC如果是对称性保护的kx一旦离开Γ点Q会立刻掉下来如果BIC在某个非零kx处那么Q会随着kx接近那个点持续上升。动量空间里Q必须有一个“发散的尖刺”这是BIC最硬的特征。第三近场与远场的能量分配。画电场模的切面图BIC的模式能量局域在结构内周围和远场几乎看不到与结构内同量级的场强然后用全局积分算Poynting矢量的净辐射功率BIC模式的净辐射功率应该趋近于数值零。这三条都满足才能有把握地说这是BIC而不是普通高Q谐振。4. 动态调节实操从BIC到quasi-BIC的Q调控以及动量空间色散的绘制4.1 打破对称性的扰动方式与Q∝δ^-2规律的复现对称性保护BIC的好处是它不依赖精细调节只要对称性在Γ点就有一个完美的BIC。那么动态调节的第一种玩法就出来了引入一个很小的几何扰动δ破坏对称性原来的BIC变成准BICQ值随δ增大而下降。具体操作上我建议用一个最简单的扰动模型——把圆盘改成椭圆。定义两个半径R_x和R_yR_x r(1δ)R_y r(1-δ)参数δ就是椭圆度。当δ0时是完美圆盘BIC在Γ点δ从0.001开始增大对称性逐级破坏。此时C2旋转对称性仍然存在但面内的镜像对称性被破坏远场辐射通道被打开。扫描方式上COMSOL里用参数化扫描直接把δ设为扫描参数从0.001到0.05步长可以取对数分布比如0.001、0.002、0.005、0.01、0.02、0.05。每个δ值都算一次特征频率记录目标模式的实部和虚部。计算完之后把Q值对δ画在双对数坐标里理论上会得到一条斜率为-2的直线这就是经典的Q∝δ⁻²标度律。为什么是-2可以这样理解δ破坏对称性之后模式向连续谱泄漏的通道打开泄漏强度正比于“扰动带来的远场偶极矩”的平方偶极矩本身正比于δ所以辐射率正比于δ²Q因子反比于δ²。如果你的扰动方式选择的是打破C2对称性比如在圆盘中心加一个偏心小孔标度律可能变成δ⁻⁴这也是一种值得复现的现象但和δ⁻²的物理来源不同。这一步算完你手上就有了第一张关键图Q随δ变化的关系图。这张图能直接说明动态调节的范围——想得到Q1E4的准BICδ取多少才合适。4.2 扫描Bloch波矢绘制带结构/动量空间色散图第二种动态调节的思路是在动量空间中移动激发点这要求把BIC的完整色散关系算出来。做法是固定几何参数δ在Floquet周期条件的k矢量里填入扫描变量kx和ky然后沿布里渊区特定路径扫描。以正方格子为例标准路径是Γ(0,0)到X(π/a,0)再到M(π/a,π/a)。每个k点都做一次特征频率求解记录目标模式的频率实部和Q。在COMSOL里这一步对应的还是参数化扫描只不过这次扫描的是kx、ky。间距可以先用粗扫描比如沿路径取20个点找Q的尖峰位置确定BIC的大致坐标后再在那个点附近加密用小步长把Q曲线的峰值形状分辨出来。把所有k点对应的特征频率画在一起横轴取k在路径上的位置纵轴取频率得到的就是带结构。BIC的位置在图上表现为一条带的下边缘或者带内交叉处出现一个点但光看带结构和普通带边区分度不高——关键是把Q值作为第三维度或者颜色映射标进去。我习惯把log10(Q)值用颜色映射在色散曲线上BIC点附近颜色会亮得发白因为Q高到几百万普通模式色散段颜色灰暗一眼就能认出BIC在哪里。这张图就是“动量空间中的BIC位置”的可视化成果。实际操作里要注意特征频率求解在每个k点都要跑一遍完整的3D特征值问题20个点就是20次求解单次几分钟的话总计小半天。为了效率第一阶段用较粗糙网格粗扫定位之后再加密网格细扫BIC附近这样可以省不少时间。如果你有两个模式在带宽内靠得很近还要小心模式交叉处特征频率的实部可能发生“换支”后处理时要按Q值连续性而不是频率连续性来追踪同一个模式。4.3 用频域研究补充远场响应Fano线型与透射谱特征频率研究算的是模式的参数但实验上看到的通常是共振谱。为了把仿真推向实验可测需要换到频域研究计算单位晶胞的透射和反射。这一步的模型在几何上不变边界设置要调整一下水平方向仍然用Floquet周期条件垂直方向把PML撤掉换成两个端口Port一个作为入射端口一个作为透射端口。端口设置里要指定激发的是哪个衍射级次因为当晶格周期较大、入射角度斜的时候可能出现高阶衍射级。一个实用技巧先按最低阶平面波入射0th order跑一版看透射谱里目标共振峰附近有没有额外的“台阶”如果有多半是高阶衍射在作祟那就需要补设高阶端口模式否则透射谱不干净。频域扫描得到透射谱后准BIC的表现是一条不对称的Fano线型共振峰旁边的背景透过率不为零峰谷两侧明显不对称。改变δ参数再扫一组透射谱会看到明显的趋势——δ越大Fano线越宽峰谷不对称性越显著δ趋于0线宽越来越窄最后变成一条几乎隐形的窄线。把谱线半峰全宽用Q f0/Δf换算一下再和特征频率法里提取的Q放在一起对确认两条路线结果一致之后整条仿真链路就算闭环了。5. 实际仿真中的高频问题与后处理经验5.1 数值假模与网格收敛性的甄别跑BIC仿真最难过的事就是算出来的模式看起来像BIC细看却不是。第一个高频问题来自PML。特征频率研究加PML之后PML层里会生成一批本征模式它们的特征频率和结构模式混在一起电场分布又不像结构模式那样有清晰的局域。解决办法是批量绘制电场模分布逐个筛选。跑一次求解十几个模式其实真正有物理意义的只有三四个剩下全在PML或者边缘区域。我自己的习惯是先在模型树上加一个“电场模”的切面图把解出的模式按频率升序排列每个模式截一张图存下来然后批量翻看。这个过程看着机械但能有效避免最后整理数据时拿到一整套不能用的模式列表。第二个高频问题来自网格收敛。Q值如果对网格尺寸特别敏感说明网格还没有收敛到数值真解。我每选一个网格密度就跑一次目标模式的Q连续两次加密后Q值变化小于5%才继续参数扫描。特别提醒一句不同δ值下达到收敛所需的网格密度可能不相同最好在最细网格下抽两三个δ点做校验。第三个问题是模式追踪断裂。参数扫描里随着δ或者波矢变化模式的实部频率可能在相邻扫描点之间跳变导致你在连续追踪Q曲线的过程中突然出现一个空洞或者跳档。解决办法在后处理阶段就预留好特征频率扫描时让COMSOL输出模式序号、实部、虚部三个量后处理软件里按实部排序、按虚部连续性检查必要时手动画“模式指纹”辅助追踪。5.2 数据导出、色散图与论文级后处理把COMSOL算完的结果整理成论文级别的图工作量不比建模小但方法对了是可以流程化的。参数化扫描跑完之后利用COMSOL的“全局计算”节点把特征频率的实部、虚部、Q值一次性导出成表格同时附上扫描参数像δ和kx。导出之后我习惯用Python的matplotlib整理数据简单脚本几行就能完成色散图和Q曲线的绘制。给你一个我常用的后处理骨架import pandas as pd import matplotlib.pyplot as plt import numpy as np df pd.read_csv(eigen_data.csv) df[Q] df[freq_real] / (2 * abs(df[freq_imag])) df[logQ] np.log10(df[Q]) plt.figure(figsize(7,5)) plt.scatter(df[kx_path], df[freq_real], cdf[logQ], cmapviridis, s30) plt.ylabel(Frequency (Hz)) plt.xlabel(k along Γ-X-M path) plt.colorbar(labellog10(Q)) plt.savefig(bandstructure_Q_colored.png, dpi300)这类图的优势一眼可见横轴是布里渊区路径纵轴是频率颜色表达Q的强弱BIC位置的高Q尖峰清清楚楚。如果再补充一组电场分布的截面图就构成了“带结构高Q模式实空间分布”的完整证据链。还有一个小习惯做参数化扫描前先在文件名里把扫描条件写清楚比如ellipse_delta0.01_kx0.0后期拼图或者复盘的时候少走很多弯路。对了频域透射谱的后处理和特征频率数据要打通。透射谱里Fano线型拟合是另一个繁琐环节好在Python的scipy里curve_fit能直接做不对称线型拟合拟合出线宽再换成Q和特征频率法互相验证。这一步能挡住相当一部分因为模式选错导致的伪Q值。最后说一个我自己的心得每次算完一个参数点我会把特征频率实部、虚部连同δ、kx一起存成一行文本攒成CSV之后画Q曲线、色散图都特别省事。COMSOL的导出功能用熟之后这套流程基本可以在两三天内跑完一轮完整的参数扫描。如果复现的时候遇到Q值死活提不上去先沉住气检查网格再看模式有没有选对最后才怀疑模型本身——排查顺序反了往往绕一大圈回到原点。 
长河君 · 2026/10/9 9:41:06 · 阅读 1.2k COMSOL超表面仿真：连续谱束缚态动量空间分析与Q因子动态调节
做超表面纳米光子学方向的人大概都绕不开连续谱束缚态BICs这个话题。我第一次在COMSOL里复现文献中的准BIC共振时最头疼的不是模型建不出来而是结构参数明明对齐了算出来的Q值却比论文里低了两个数量级——后来才明白问题出在我对“动量空间”这个自由度的理解不够。这篇内容我不打算按软件手册的顺序讲COMSOL基本操作而是从BICs在光子晶体超表面中的动量空间行为出发把动态调节的完整仿真流程拆开了讲覆盖单位晶胞建模、特征频率模式筛选、Q因子提取、带结构扫描和常见的数值坑。适合刚开始用COMSOL做超表面研究、已经能跑通频域透射谱但想进一步调BIC的研究生和工程师。1. 先搞明白BIC在动量空间里是什么状态动态调节调的又是哪个自由度1.1 对称性保护BIC与偶然BIC的本质区别连续谱束缚态这个名字听起来玄但它要表达的东西很直接一个模式的频率落在了连续辐射谱的范围里却不向外辐射能量。放在光子晶体超表面的语境里就是结构支持某个本征模式它的频率和周围空间能传播的平面波频率重叠但因为对称性或干涉相消的原因模式和所有辐射通道的耦合都归零能量就被“锁”在结构内不往外跑。文献里最常见的分类是两类。第一类叫对称性保护BIC通常出现在布里渊区高对称点比如Γ点。它的存在不依赖精细调节只要结构保持某种面内镜像对称、旋转对称模式远场的极性和可辐射平面波不匹配就能形成BIC。类比一下就像一个人站在旋转门正中间四个门扇对称分布他往哪个方向推门都推不开——对称性让所有通道都失效了。第二类叫偶然BIC也叫Friedrich-Wintgen型BIC它依赖两个模式之间的干涉相消需要在某个特定波矢、特定结构参数下才会出现调节精度要求更高。这两类BIC在动量空间里的表现很不一样。对称性保护BIC的位置基本钉死在Γ点参数变了它也不会挪窝偶然BIC则可能出现在Γ点以外结构参数一变它在动量空间里的位置还会漂移。做动态调节之前必须先把目标BIC归好类因为后面选参数、扫波矢、解释Q值变化趋势都要靠这张“身份卡”。1.2 动态调节的三种物理路径理解了BIC的本质再看“动态调节”就清晰了——BIC一旦受到扰动有辐射通道被打开就会变成准BICquasi-BICQ因子从一个理论上无穷大的数值跌到有限值。调节的手段大致有三种对应到COMSOL里就是三套不同的参数扫描流程。第一种是结构几何扰动最常用。比如把纳米圆盘改成椭圆或者在圆盘旁边加一个小孔、小柱用一个无量纲参数δ描述破缺程度。对称性保护BIC对几何扰动最敏感δ从0变到0.1的时候Q值可能从百万级降到百级跨了好几个数量级。第二种是入射角扫描也就是在动量空间里换个方向激发。超表面的能带结构里BIC存在于特定的平行波矢(kx,ky)处实验中通过改变斜入射角度可以让激发条件刚好落在BIC附近谱线从Fano线型逐渐收窄Q值在特定角度达到峰值。第三种是材料动态调谐比如改变周围介质的折射率、给结构加载电压调控载流子浓度能带整体移动BIC在动量空间的位置也会跟着漂。从仿真角度看第一种用参数化扫描结构尺寸第二种用Floquet周期条件里的波矢扫描第三种用材料折射率参数扫描。三者经常组合使用但核心逻辑都是一样的先找一个完好结构确认BIC存在再引入可控扰动观察Q因子和共振波长的变化规律。2. COMSOL建模从单位晶胞到可参数化光子晶体超表面2.1 几何、材料与周期的参数化设置COMSOL里做周期结构超表面标准路线是RF模块下的电磁波物理场研究类型选“特征频率”或“频域”。模型只要一个单位晶胞就够了水平方向用周期边界垂直方向加空气间隔和吸收边界这和实验里整个阵列对应起来理解晶胞无限重复等效于一个无限大的周期平面。拿一个最经典的例子做底子正方形晶格周期a 700 nm硅纳米盘半径r 180 nm高度h 200 nm背景是空气。硅的折射率取3.48在近红外波段色散不大可以先按常数算。COMSOL里参数表的设置我习惯这样组织参数名值说明a700[nm]晶格周期r180[nm]圆盘半径h200[nm]圆盘高度delta0几何扰动幅度扫描用kx0平行波矢x分量扫描用ky0平行波矢y分量扫描用n_Si3.48硅的折射率n_bg1.0背景折射率f0c/(0.7[um]*2.4)特征频率搜索基准值最后一项f0是我自己的习惯先估算一个大概的共振频率特征频率研究里把搜索范围放在它附近能省掉大量无关模式。估算思路很简单模式在硅和空气中混合分布等效折射率大概在2.2到2.6之间硅盘的有效波长大约是700nm除以2.4共振频率就在c/(0.7μm×2.4)附近算出来约等于1.78E14 Hz。实际跑的时候可以放宽一点左右各留20%的余量。把结构尺寸、周期、扰动参数全部建成全局参数后面做参数化扫描会非常省事不用每次改几何。2.2 Floquet周期条件和PML的放置逻辑模型里最关键的两个边界设置一个是水平方向的周期条件一个是垂直方向的开放边界。水平方向把单位晶胞的x和y方向的相对面分别设成周期条件Periodic Condition类型选Floquet周期性。COMSOL里会要求填k矢量也就是平行方向的Bloch波矢。这里要注意周期条件只能设一对面对x方向一对、y方向一对一共两对每个周期条件里都要输入kx和ky。当kxky0时对应正入射也就是Γ点想算带结构或者斜入射就得把这两个参数变成扫描变量。垂直方向也就是晶胞的上方和下方要考虑能量可以往外辐射。最稳妥的做法是加一层空气间隔再接完美匹配层PML。我给一个具体尺寸建议硅盘上下各加0.5λ的空气间隔PML厚度做成1λλ取模型中心频率对应的自由空间波长。PML本身在特征频率计算里会引入一批非物理模式能量集中在PML区域里后面选模式的时候要明确把它们剔掉。判别方法也简单画电场分布图看能量是不是一团一团地“糊”在PML层里——只要是这种模式一律不选。有些做法是用散射边界条件代替PML在纯BIC状态下问题不大因为模式本身不辐射但算准BIC时Q值往往在1E3到1E5之间散射边界对弱辐射的微小反射也会污染结果我实际对比下来还是PML更稳定。2.3 网格策略为什么BIC模拟对网格一致性要求极高网格可能是BIC仿真里最容易被低估的一环。BIC的Q值本质上是数值里虚部的倒数Q越高意味着模式往外漏的能量越小——小到和网格离散误差同一个量级。如果网格太粗数值损耗会“淹没”真实的辐射损耗算出来的Q光子晶体的Q值远低于理论值。反过来网格波动也会制造虚假的泄漏通道Q表现得比实际低。我的经验是硅盘内部网格最大尺寸取有效波长的1/8到1/10空气区域最多取λ/8。以700nm周期、硅折射率3.48为例硅内部光的有效波长约是201nm那硅盘内网格最大尺寸要控制在25nm左右空气区域控制在90nm左右。这个量级下的3D模型自由度大概在几十万到上百万一台16GB内存的工作站还能跑得动。但网格策略真正的核心是“一致”。做参数扫描的时候如果扰动参数一改网格也跟着自动加密Q值的变化里就混进了网格带来的变化干扰很大。我的做法是把几何扰动设计成不改变网格拓扑的形式比如圆盘变椭圆只靠一个参数控制半径网格重画但密度分布基本一致。如果结构里加小孔这种导致局部网格剧变的操作要专门对孔周边加固定尺寸约束保证孔在不同尺寸下网格密度近似。收敛性检查也别省把全局网格尺寸整体缩小一半选同一个模式看Q值变化偏差超过5%就说明网格还没收敛加密之后重算。3. 用特征频率研究锁定BIC模式复数特征值中读出Q因子的方法3.1 特征频率研究的设置与模式筛选模型建好后研究类型选“特征频率”。这一步的目的是直接求解模式的本征频率复数特征值里实部是共振频率虚部是辐射损耗的阻尼Q因子等于两者比值。求解器设置里有几个关键选项。特征频率搜索范围我给一个宽区间通常在估算的f0前后扫5E13 Hz每次求解10到20个模式。这个数量确实会带来大量无关模式但BIC必须找了才知道在哪里宁可多算不可漏算。物理场设置里电磁波特征频率研究默认解的是角频率单位是rad/s这一点对后面算Q尤其重要先确认设置里特征频率单位是角频率还是频率别算完发现差了6个数量级。算完之后后处理是对着一组特征频率列表逐个看电场分布。BIC模式的典型电场图是能量极强地局域在硅盘内部或者表面向外场衰减极快PML里几乎没有响应。对比之下普通辐射模式在PML区域还会有明显的振荡图样。我在最开始跑的时候犯过糊涂看到某个模式的电场分布很规整一度以为它是BIC后来一查Q值只有200再对比频域透射谱发现根本对应不上——直觉不能替代数值判据。3.2 Q因子计算公式与数值工具验证Q因子的计算式子很直接Q ω_real / (2|ω_imag|)如果COMSOL输出的是角频率直接代入实部和虚部单位会自动约掉如果是频率同样可以用f_real/(2|f_imag|)两种方式结果一致。虚部的符号一般是负的表示能量损耗取绝对值即可。举个例子假设某一个准BIC模式特征频率实部是1.414E15 rad/s虚部是7.07E9 rad/s那Q就是1.414E15除以1.414E10算出来约等于10万。这是一个可参考的量级完全对称结构下的BICQ值往往能到1E7甚至更高数字上的限制主要来自数值精度而加了几个百分比的扰动之后Q快速回落到1E3到1E5这是频域实验里能测得比较准的范围。特征频率法算的Q值不是一个结果就结束了我通常会在同一次参数扫描里配合频域研究做交叉验证。频域扫描里画透射谱谱线在共振峰附近呈现Fano线型共振峰的宽度Δλ对应Q λ0/Δλ。把两种方法算出来的Q放在一张双对数图里对比吻合的区间才能放心写进论文里。频域谱的好处是它直接对应实验可测的量特征频率法快但抽象两条路互相印证能避免特征频率模式下数值噪声造成的假结果。3.3 判断“这个模式是不是BIC”的三个判据初学者问得最多的问题就是跑完特征频率之后怎么从一长串模式里认出哪个是BIC。我给三个判据按重要性排序。第一Q值数量级。BIC对应的是模式虚部接近零Q在数值噪声允许的范围内趋于无穷。对称结构下Q值至少比周围普通辐射模式高两个数量级以上如果所有模式Q值都停在几百那大概率还没找到BIC或者网格不够细。第二动量空间的发散性。把kx设成0时高Q的那个模式往kx的正方向扫几个小步长观察Q的变化。BIC如果是对称性保护的kx一旦离开Γ点Q会立刻掉下来如果BIC在某个非零kx处那么Q会随着kx接近那个点持续上升。动量空间里Q必须有一个“发散的尖刺”这是BIC最硬的特征。第三近场与远场的能量分配。画电场模的切面图BIC的模式能量局域在结构内周围和远场几乎看不到与结构内同量级的场强然后用全局积分算Poynting矢量的净辐射功率BIC模式的净辐射功率应该趋近于数值零。这三条都满足才能有把握地说这是BIC而不是普通高Q谐振。4. 动态调节实操从BIC到quasi-BIC的Q调控以及动量空间色散的绘制4.1 打破对称性的扰动方式与Q∝δ^-2规律的复现对称性保护BIC的好处是它不依赖精细调节只要对称性在Γ点就有一个完美的BIC。那么动态调节的第一种玩法就出来了引入一个很小的几何扰动δ破坏对称性原来的BIC变成准BICQ值随δ增大而下降。具体操作上我建议用一个最简单的扰动模型——把圆盘改成椭圆。定义两个半径R_x和R_yR_x r(1δ)R_y r(1-δ)参数δ就是椭圆度。当δ0时是完美圆盘BIC在Γ点δ从0.001开始增大对称性逐级破坏。此时C2旋转对称性仍然存在但面内的镜像对称性被破坏远场辐射通道被打开。扫描方式上COMSOL里用参数化扫描直接把δ设为扫描参数从0.001到0.05步长可以取对数分布比如0.001、0.002、0.005、0.01、0.02、0.05。每个δ值都算一次特征频率记录目标模式的实部和虚部。计算完之后把Q值对δ画在双对数坐标里理论上会得到一条斜率为-2的直线这就是经典的Q∝δ⁻²标度律。为什么是-2可以这样理解δ破坏对称性之后模式向连续谱泄漏的通道打开泄漏强度正比于“扰动带来的远场偶极矩”的平方偶极矩本身正比于δ所以辐射率正比于δ²Q因子反比于δ²。如果你的扰动方式选择的是打破C2对称性比如在圆盘中心加一个偏心小孔标度律可能变成δ⁻⁴这也是一种值得复现的现象但和δ⁻²的物理来源不同。这一步算完你手上就有了第一张关键图Q随δ变化的关系图。这张图能直接说明动态调节的范围——想得到Q1E4的准BICδ取多少才合适。4.2 扫描Bloch波矢绘制带结构/动量空间色散图第二种动态调节的思路是在动量空间中移动激发点这要求把BIC的完整色散关系算出来。做法是固定几何参数δ在Floquet周期条件的k矢量里填入扫描变量kx和ky然后沿布里渊区特定路径扫描。以正方格子为例标准路径是Γ(0,0)到X(π/a,0)再到M(π/a,π/a)。每个k点都做一次特征频率求解记录目标模式的频率实部和Q。在COMSOL里这一步对应的还是参数化扫描只不过这次扫描的是kx、ky。间距可以先用粗扫描比如沿路径取20个点找Q的尖峰位置确定BIC的大致坐标后再在那个点附近加密用小步长把Q曲线的峰值形状分辨出来。把所有k点对应的特征频率画在一起横轴取k在路径上的位置纵轴取频率得到的就是带结构。BIC的位置在图上表现为一条带的下边缘或者带内交叉处出现一个点但光看带结构和普通带边区分度不高——关键是把Q值作为第三维度或者颜色映射标进去。我习惯把log10(Q)值用颜色映射在色散曲线上BIC点附近颜色会亮得发白因为Q高到几百万普通模式色散段颜色灰暗一眼就能认出BIC在哪里。这张图就是“动量空间中的BIC位置”的可视化成果。实际操作里要注意特征频率求解在每个k点都要跑一遍完整的3D特征值问题20个点就是20次求解单次几分钟的话总计小半天。为了效率第一阶段用较粗糙网格粗扫定位之后再加密网格细扫BIC附近这样可以省不少时间。如果你有两个模式在带宽内靠得很近还要小心模式交叉处特征频率的实部可能发生“换支”后处理时要按Q值连续性而不是频率连续性来追踪同一个模式。4.3 用频域研究补充远场响应Fano线型与透射谱特征频率研究算的是模式的参数但实验上看到的通常是共振谱。为了把仿真推向实验可测需要换到频域研究计算单位晶胞的透射和反射。这一步的模型在几何上不变边界设置要调整一下水平方向仍然用Floquet周期条件垂直方向把PML撤掉换成两个端口Port一个作为入射端口一个作为透射端口。端口设置里要指定激发的是哪个衍射级次因为当晶格周期较大、入射角度斜的时候可能出现高阶衍射级。一个实用技巧先按最低阶平面波入射0th order跑一版看透射谱里目标共振峰附近有没有额外的“台阶”如果有多半是高阶衍射在作祟那就需要补设高阶端口模式否则透射谱不干净。频域扫描得到透射谱后准BIC的表现是一条不对称的Fano线型共振峰旁边的背景透过率不为零峰谷两侧明显不对称。改变δ参数再扫一组透射谱会看到明显的趋势——δ越大Fano线越宽峰谷不对称性越显著δ趋于0线宽越来越窄最后变成一条几乎隐形的窄线。把谱线半峰全宽用Q f0/Δf换算一下再和特征频率法里提取的Q放在一起对确认两条路线结果一致之后整条仿真链路就算闭环了。5. 实际仿真中的高频问题与后处理经验5.1 数值假模与网格收敛性的甄别跑BIC仿真最难过的事就是算出来的模式看起来像BIC细看却不是。第一个高频问题来自PML。特征频率研究加PML之后PML层里会生成一批本征模式它们的特征频率和结构模式混在一起电场分布又不像结构模式那样有清晰的局域。解决办法是批量绘制电场模分布逐个筛选。跑一次求解十几个模式其实真正有物理意义的只有三四个剩下全在PML或者边缘区域。我自己的习惯是先在模型树上加一个“电场模”的切面图把解出的模式按频率升序排列每个模式截一张图存下来然后批量翻看。这个过程看着机械但能有效避免最后整理数据时拿到一整套不能用的模式列表。第二个高频问题来自网格收敛。Q值如果对网格尺寸特别敏感说明网格还没有收敛到数值真解。我每选一个网格密度就跑一次目标模式的Q连续两次加密后Q值变化小于5%才继续参数扫描。特别提醒一句不同δ值下达到收敛所需的网格密度可能不相同最好在最细网格下抽两三个δ点做校验。第三个问题是模式追踪断裂。参数扫描里随着δ或者波矢变化模式的实部频率可能在相邻扫描点之间跳变导致你在连续追踪Q曲线的过程中突然出现一个空洞或者跳档。解决办法在后处理阶段就预留好特征频率扫描时让COMSOL输出模式序号、实部、虚部三个量后处理软件里按实部排序、按虚部连续性检查必要时手动画“模式指纹”辅助追踪。5.2 数据导出、色散图与论文级后处理把COMSOL算完的结果整理成论文级别的图工作量不比建模小但方法对了是可以流程化的。参数化扫描跑完之后利用COMSOL的“全局计算”节点把特征频率的实部、虚部、Q值一次性导出成表格同时附上扫描参数像δ和kx。导出之后我习惯用Python的matplotlib整理数据简单脚本几行就能完成色散图和Q曲线的绘制。给你一个我常用的后处理骨架import pandas as pd import matplotlib.pyplot as plt import numpy as np df pd.read_csv(eigen_data.csv) df[Q] df[freq_real] / (2 * abs(df[freq_imag])) df[logQ] np.log10(df[Q]) plt.figure(figsize(7,5)) plt.scatter(df[kx_path], df[freq_real], cdf[logQ], cmapviridis, s30) plt.ylabel(Frequency (Hz)) plt.xlabel(k along Γ-X-M path) plt.colorbar(labellog10(Q)) plt.savefig(bandstructure_Q_colored.png, dpi300)这类图的优势一眼可见横轴是布里渊区路径纵轴是频率颜色表达Q的强弱BIC位置的高Q尖峰清清楚楚。如果再补充一组电场分布的截面图就构成了“带结构高Q模式实空间分布”的完整证据链。还有一个小习惯做参数化扫描前先在文件名里把扫描条件写清楚比如ellipse_delta0.01_kx0.0后期拼图或者复盘的时候少走很多弯路。对了频域透射谱的后处理和特征频率数据要打通。透射谱里Fano线型拟合是另一个繁琐环节好在Python的scipy里curve_fit能直接做不对称线型拟合拟合出线宽再换成Q和特征频率法互相验证。这一步能挡住相当一部分因为模式选错导致的伪Q值。最后说一个我自己的心得每次算完一个参数点我会把特征频率实部、虚部连同δ、kx一起存成一行文本攒成CSV之后画Q曲线、色散图都特别省事。COMSOL的导出功能用熟之后这套流程基本可以在两三天内跑完一轮完整的参数扫描。如果复现的时候遇到Q值死活提不上去先沉住气检查网格再看模式有没有选对最后才怀疑模型本身——排查顺序反了往往绕一大圈回到原点。