做BICs和光子晶体超表面计算有一段时间了最让我头疼也最上头的就是复现文献里那种“围绕BIC动量点的远场偏振矢量涡旋图”。一句话说清楚的话拓扑BICs连续谱束缚态在动量空间是一个偏振奇点远场偏振矢量会围绕它旋转一圈转多少圈就是拓扑荷的数值。可文献通常只给结论和漂亮的示意图不会告诉你COMSOL里到底怎么建模型、怎么扫动量空间、怎么把一个本征模的复电场变成一张能积出拓扑荷的偏振图。这篇我把整个流程从物理概念到COMSOL操作、再到后处理算拓扑荷的完整链条拉通包括我在这个过程中踩过的坑。适合正在用COMSOL做超表面、光子晶体能带与BIC研究的研究生和工程师也适合想复现文献里偏振矢量图、确认自己算的是不是真正BIC的人。1. 拓扑BICs的物理图景为什么远场偏振矢量图能定义“电荷”1.1 连续谱束缚态不是“没有损耗”很多人第一次接触BICs会把它误解成“品质因子无穷大的共振模”。这说法对了一半但容易让人忽略核心物理也容易在COMSOL里犯错误。我建议把BIC理解为一个本可以在连续辐射通道中漏泄出去的模式却因为某种破坏性干涉、对称性或拓扑约束向外辐射的幅度严格为零。也就是说这个模和自由空间里传播波之间没有耦合通道因此在无限大结构里它不衰减Q因子严格发散。在光子晶体超表面里面BIC通常落在某个衍射级次恰好可以开启/关闭的临界点附近或者落在周期性结构中某个布洛赫波矢下辐射本征矢出现零投影的位置。常见的有两类对称保护BIC出现在高对称点比如正方晶格Γ点靠晶格群的不可约表示和远场平面波的对称性失配来关断辐射还有偶然BIC也叫Friedrich–Wintgen型BIC出现在非高对称动量点是两个或多个模式的远场辐射路径之间发生相消干涉导致的。无论哪种你去看它的远场偏振分布时都会发现围绕BIC的那个动量点偏振方向旋转了一圈旋转圈数就是拓扑荷。这一点是BIC和普通高Q峰的本质区别——普通高Q峰即使Q很高周围偏振场也没有涡旋结构。1.2 偏振涡旋与拓扑荷的物理含义把超表面每个出射方向的透射或辐射场记为Jones矢量也就是复的横向电场分量Ex, Ey。每给一个动量点k(kx, ky)就对应一个Jones矢量。把这个矢量在k空间里画成箭头围绕BIC所在点走一圈闭合路径箭头方向会以一个固定规律旋转。如果绕一圈后偏振椭圆的长轴方向转过π拓扑荷就是±1如果转过mπ拓扑荷就是±m。这个“电荷”不是电磁荷它衡量的是偏振矢量场在动量空间里的绕数。它的价值在于鲁棒性——拓扑荷是个整数不能连续变化所以只要BIC周围的远场偏振矢量场连续变化围绕BIC闭合回路的拓扑荷就不会因为结构参数微扰而改变。这就是为什么很多BIC能够在结构对称性破缺后依然存活只是从Γ点挪到了非零的动量位置。真正改变拓扑荷需要让一对正负拓扑荷的BIC在动量空间相遇、相互湮灭这也是近年来拓扑BIC研究里最有趣的现象之一。1.3 你需要从仿真中提取什么要在COMSOL里把这件事落地核心需求可以拆成三块。第一找BIC的位置。通常做法是在动量空间扫描计算某个带的Q因子Q因子发散的那个点就是BIC候选点。第二围绕BIC取一个闭合k空间路径比如以BIC为圆心的圆在这个路径上等间距取点。第三对每个k点得到复数Jones矢量也就是横向电场分量的幅度和相位然后画偏振矢量图算拓扑荷。整个过程在COMSOL里的困难点不是单个步骤有多难而是每一步都有容易做得“差不多”但又不对的地方。特别是Q因子发散点周围偏振变化通常非常剧烈网格、频率扫描间隔、相位展开稍微处理不好算出来的拓扑荷就不是整数这时候第一反应不该是怀疑物理模型而是去排查数据处理链路。2. 拓扑荷的数学定义与远场偏振图的构成逻辑2.1 Jones矢量与偏振椭圆取向角一个沿z方向传播的平面波其电场横向分量写成复数形式Ex |Ex| exp(iφx)Ey |Ey| exp(iφy)这个Jones矢量Ex, Ey就携带了偏振的全部信息。实际画偏振矢量图时我们通常画偏振椭圆长轴方向因为这个方向在近线性偏振区域里最直观、也最容易连续追踪。偏振椭圆长轴与x轴的夹角ψ不是简单的arg(Ey/Ex)因为椭圆偏光的长轴方向还需要考虑两个分量的相位差和幅度比。用到斯托克斯参数会更方便S1 |Ex|² − |Ey|² S2 2 Re(Ex* · Ey)这里的Ex*是Ex的共轭。长轴取向角由下式唯一确定ψ 0.5 * atan2(S2, S1)为什么要用atan2而不是atan因为S1和S2同时为负时ψ可能在第二象限普通atan只给四分之一象限的值会直接把偏振矢量图搅乱。atan2是唯一没有这个歧义的写法。2.2 环绕积分与离散采样计算定义好ψ(k)之后拓扑荷就是一个关于动量空间的闭合积分q (1 / π) ∮_C ∇_k ψ(k) · dk等价地在离散采样的情况下沿逆时针方向对一圈N个采样点把相邻点的ψ差加起来注意相位展开要按2ψ做最后再除以2πq (1 / 2π) ∮_C d[arctan2(S2, S1)]这里为什么是除以2π而不是π因为ψ本身已经是从atan2(…)里取出来的半个角度arctan2(S2,S1)的相角变化一整圈是2π而ψ对应这个相角的一半所以拓扑荷等于这个相角变化除以2π。换成ψ的话就是说绕一圈ψ累计变化πq。真正写代码时我习惯先构造一个复数Z S1 i·S2 (|Ex|² − |Ey|²) i·2·Re(Ex* Ey)然后用np.unwrap(np.angle(Z))做相位展开最后取总变化除以2π。这样比直接展开ψ少了边界判断也更稳。2.3 符号约定与边界情况符号问题是新手最容易翻车的地方也是我不会回避的部分。同样的物理系统不同的文献可能给出q1或q−1的BIC原因往往只是路径方向约定不同沿着动量空间逆时针环绕为正算出来是什么就是什么但有的作者用的是顺时针为正有的作者还额外对y轴做了镜像翻转。我的建议是先固定自己的约定写进代码注释里然后拿一个已知BIC当样例验证符号。比如规范文献中某对称保护BIC的拓扑荷是1你的流程复现出来如果是−1就把路径方向或旋转矩阵翻转一下然后整个项目保持统一。这不算作弊因为拓扑荷的符号本来就依赖坐标方向约定重要的是同一种约定下符号一致、可复现。还有一个值得注意的边界情况如果闭合路径上恰好采样到了某种圆偏振点C点这时候偏振椭圆退化成圆长轴方向没有定义ψ会出现跳变。BIC本身通常不是C点但它的附近可能有C点出现采样时如果某个k点非常靠近C点就会导致拓扑荷计算出现伪跳变。遇到这种情况要么增加采样密度要么让闭合路径偏离C点位置。3. COMSOL超表面模型搭建几何、材料与边界条件的取舍3.1 几何与参数设计示例参数在COMSOL里建光子晶体超表面模型我习惯先用简单几何把整个流程跑通再逐步增加复杂度。以硅纳米柱正方晶格为例一组可以直接起步的参数是晶格周期a 800 nm圆柱高度h 400 nm圆柱半径r 150 nm硅的折射率n_si 3.45环境是空气n_air 1。这个结构在近红外波段有一个Γ点附近的中频带适合做对称保护BIC演示。如果你要研究偶然BIC可以后面把圆柱改成椭圆柱短半轴沿x或y方向通过改变椭圆度让BIC离开Γ点。注意一个原则先算无衬底的板在整个流程验证无误后再加二氧化硅衬底。加了衬底后上下不再对称很多原本在对称结构里清晰的模式会发生劈裂也会引入新的衍射级次新手很难判断Q因子发散是真实BIC还是模式劈裂造成的假象。3.2 物理场接口与周期边界设置模型导航里选三维物理场用“波动光学 电磁波、频域”也就是ewfd接口。研究类型有两种选择定位BIC用“特征频率”研究提取远场偏振矢量图用“频域”研究加周期端口。很多人在同一个模型里试图一次搞定所有事情我的经验是分两个研究步骤更清爽先用特征频率扫动量空间定位BIC锁定位置后再单独跑频域端口模型提取偏振数据。边界条件方面x和y方向用周期性条件类型选“Floquet周期性”布洛赫波矢分量kx和ky定义为全局参数方便扫描。在COMSOL的周期性条件里两个方向上的波矢可以有不同符号习惯我个人在设置时会把x方向写成exp(-i·kx·x)对应的约定这样和大多数光子晶体文献的布洛赫定理形式一致。z方向用散射边界条件或完美匹配层PML。具体操作上我建议把波矢分量用归一化形式声明KX 0归一化布洛赫波矢范围通常取[−0.5, 0.5]KY 0在周期性条件的波矢设置里把x方向布洛赫波矢填成KX * pi / a注意COMSOL里波矢单位可能需要对应弧度/米y方向同理。这样扫描时只需改KX和KY数值上更直观。3.3 网格与求解器的关键配置网格是COMSOL模拟里最影响BIC结果的部分。BIC的核心特征远场投影为零数值上如果网格不够精细这个零会被打散出现伪泄漏。我的配置习惯是圆柱内部用最大单元尺寸λ_eff/8左右空气区域最大单元λ0/10靠近圆柱和边界处的边界层用三层以上。虽然计算量大一些但特征频率法扫几十个k点是可以接受的。关键点在于周期性边界的网格必须严格匹配。也就是x0面和xa面上的网格分布要一致y方向同理。网格不匹配会在本应周期性连续的位置引入非物理反射导致Q因子被严重低估。求解器方面特征频率研究要指定待找模式个数一般填“附近模式数”8-12个并锁定搜索频率范围比如从f_center/2到2×f_center。我通常会先在Γ点附近找到一个结构场集中在柱体内部的模式确认它的电场图与理论的TM-like或TE-like模式对应然后把这个模式作为后续扫描的追踪对象。COMSOL的参数化扫描支持使用上一步解做初始值这个功能在扫动量空间时一定要打开否则每个k点重新迭代容易跳到别的带上。4. 锁定BIC的动量扫描策略从Q因子发散点到采样路径4.1 用特征频率扫描找Q因子尖峰BIC位置不一定在可猜测的高对称点尤其研究偶然BIC时定位它是整个流程中最耗时也最关键的一步。我的做法是先粗后细的两级扫描。先在动量空间做一个较低分辨率的网格例如kx∈ −0.3, 0.3 、ky∈[−0.3, 0.3]步长0.05。对每个k点用特征频率求解找到目标模式计算两个量实部频率 fr Re(f)辐射阻尼比 γi |Im(f)| / Re(f)然后画一张伪色彩图以kx为横轴、ky为纵轴、颜色的对数为Q值或1/γiBIC位置会出现一个极其尖锐的亮区。如果扫到的Q值超过10^6就说明附近极可能有BIC此时把动量范围缩小到这一步找到的极值点周围步长改到0.001量级继续细化。为什么这个方案有效因为BIC附近Q因子随动量距离呈倒二次方发散Q ~ 1/Δk²所以即使在粗网格上只要网格点距离BIC足够近Q值就会出现可辨识的异常峰。如果整个网格扫完了都没看到任何一个点的Q比周围高一个量级基本可以排除这个波段存在BIC需要调整参数重新来找。4.2 如何区分真实BIC与PML伪模态特征频率求解带PML的结构PML里也会出现一堆非物理模式这些模式的复频率分布很有迷惑性。新手照着“Q值最高”去选很容易选到一个完全局域在PML里的伪模。我识别真实BIC的经验有两条结合起来很好用。第一条是看电场分布真实的结构BIC电场能量几乎全部集中在超表面柱体附近电磁场在PML内部的能量占比应该接近于零。在COMSOL里可以做一个全局计算统计PML区域和结构区域的电场能量密度积分。如果PML区域的能量占比超过10%基本不用继续分析了。第二条是检验模式对网格细化的收敛性。真实BIC的频率和Q值会随着网格加密单调趋近某个极限值而PML伪模会乱跑。在最终定位BIC前把最佳网格粗细细化和细化两个级别各算一次如果BIC点的Q值变化小于一个数量级说明网格可以接受如果Q值还在显著上升说明之前的网格不够继续加密。4.3 扫描窗口与采样点数的实际经验找到BIC动量坐标后还要设计计算拓扑荷用的闭合路径。路径半径的选择是门学问。半径太大路径会穿过带边缘或碰到其他模式偏振矢量图变得不连续半径太小对网格和数值精度的要求陡然上升Q值巨大用频域法提取偏振时共振峰窄到难以精确采样。我常用的规则是路径包围范围内不包含其他BIC或C点半径取BIC到最近动量奇点距离的1/5到1/3。路径上采样点数的选择32点是起点推荐64点。为什么不能太少因为围绕BIC一圈的偏振矢量旋转在距离奇点很近的地方变化率最高若采样点之间跳过了一个超过90°的取向角变化后续相位展开就会出现偏差。若后续发现ψ沿路径出现不自然的大跳变先增加采样点重算不要急着改网格。为了一次扫描得到对称的闭合路径我习惯用参数化控制。定义扫描参数t从0到2π圆上的kx(t) kx0 R·cos(t)ky(t) ky0 R·sin(t)在全局参数里用表达式声明扫描t的取值。这样得到的路径天然首尾相连方便拓扑荷积分。5. 从COMSOL结果到拓扑荷偏振矢量提取与数值积分全流程5.1 提取远场复电场分量在这一步我推荐的方案是用频域端口分析而不是直接读本征模式的远场节点。原因很现实周期结构超表面的本征模式自带辐射通道但COMSOL不同版本里远场节点的定义差异较大尤其是在周期性结构里“远场”这个词容易和PML外推混淆很难保证你拿到的偏振方向就是物理上朝某一衍射级次辐射的平面波偏振。而周期端口通过S参数给出的正交偏振分量定义清晰稳定可复现性最好。模型调整为z方向上面和下面各加一个周期端口方向相反。端口模式一般有两个正交偏振模式分别对应x偏振和y偏振具体要看端口模式的极化方向。从上方端口注入一个x偏振平面波下方端口输出的两个正交偏振透射系数的复数S参数就分别是Jones矢量的Ex和Ey分量Ex S21_xxx偏振入射x偏振透射Ey S21_yxx偏振入射y偏振透射这里的关键是共振频率处的取值。对闭合路径上的每个k点先由特征频率法得到该k点模式频率fr然后在频域研究里以fr为中心扫描一个窄带范围一般设为fr ± 10×Δf其中Δf是共振线宽可用Q估算Δf fr/Q。扫描频率点数建议100-200个高Q情况下适当增加。取|S21_xx|²|S21_yx|²最大也就是透射峰处的复S参数作为该k点的Jones矢量。所有k点都采完之后把幅度相位数据导出成CSV或文本文件交给Python处理。5.2 用Python计算偏振椭圆取向角拿到每个k点的复电场(Ex, Ey)后计算拓扑荷的代码逻辑很简洁。下面这段是我一直在用的核心流程import numpy as np # coords: 采样点角度theta按逆时针排列0 - 2*pi # Ex, Ey: 复数Jones矢量分量跟theta一一对应 def topological_charge(theta, Ex, Ey, clockwiseFalse): # 斯托克斯参数 S1 np.abs(Ex)**2 - np.abs(Ey)**2 S2 2.0 * np.real(Ex * np.conj(Ey)) # 构造用于相位积分的复数 Z S1 1j * S2 phase np.angle(Z) phase_unwrapped np.unwrap(phase) # 总相位变化 dphi phase_unwrapped[-1] - phase_unwrapped[0] q dphi / (2.0 * np.pi) # 若需要顺时针则取反 return q if not clockwise else -q这段代码里最需要注意的是为什么用角度的相位展开而不是按ψ。因为Z的定义里面已经包含了偏振椭圆的所有信息它的幅角就是斯托克斯矢量在赤道平面上的方位角展开后的总变化直接除以2π就是拓扑荷。用ψ展开还需要人为乘2增加一步就增加一份出错可能。有些数据组如果采样密度太低np.unwrap本身可能也救不回来。此时我会先打印phase的变化轨迹看是否存在超过π的邻近点跳变。如果有先把采样点加密或者检查是不是碰到了C点。5.3 环路积分的相位展开技巧相位展开在实际数据处理中地位很高我可以提供一个提高鲁棒性的技巧在做np.unwrap之前先把相位数组延拓一个周期。也就是把theta从0到2π的数据复制一份接到原来的末尾形成0到4π的长序列在这份加长序列上做unwrap然后取2π到4π段的总变化。这样做的目的是避免展开算法在首尾边界处产生错误判断。而且我总是保证闭合路径首尾两个点实际上是同一个物理点。也就是说采样点数组里第0个和第N个点坐标必须完全一致带扫描时把t0和t2π都放进去这样积分闭合性才有保证。用参数化扫描时很容易做到但手动列k点表格时经常漏了最后一个点或者最后一个点跟初始点差了一个微小偏移最后算出来的q可能会有0.01量级的误差看起来不是整数一时半会儿又查不出原因。6. 我踩过的坑与校验习惯网格、相位缠绕与符号错乱6.1 三个导致拓扑荷非整数的常见原因拓扑荷连续计算出来的结果不是整数几乎可以确定是数据链路某个环节出了问题。根据我的经验最常见的有三个原因。第一个是采样密度不足。BIC附近偏振变化极快如果路径上相邻k点间的偏振方向角改变超过π/2相位展开就会判读错误。这时候把采样点从32加到64或者128问题往往就消失了。第二个是频率脱离共振中心。频域扫描提取Jones矢量时如果取的频率稍微偏离共振中心透射场还混着背景连续谱的贡献得到的Jones矢量不纯粹偏振方向就被污染了。所以每个k点都要精确找到透射峰值的位置不能拍脑袋用BIC中心频率代替所有k点的采样中心。第三个是相位参考面不一致。COMSOL的S参数相位会随端口参考面位置线性漂移如果你在几何建模时端口位置设置得不规则或者参考面之间距离没有统一那么得到的Ex和Ey会同时包含一个随波长变化的公共相位这个公共相位不影响拓扑荷因为S1和S2是二次型公共相位会抵消。但如果Ex和Ey分别用的参考面不一致就会直接影响Re(Ex*·Ey)从而污染S2。这是我第一次算出来q0.3时浪费了两天才找到的罪魁祸首。6.2 校验方法改变半径、改变点数、对比近场投影一套流程能不能信就要看它稳不稳。我常用的三层校验方法如下。第一层换闭合路径半径。围绕同一个BIC用R、2R、4R三个半径分别计算拓扑荷q应该严格相等。如果4R出了变化大概率是大半径路径触碰到了其他奇点或模式如果小半径R变了则是数值精度不够需要加密网格或增加采样点。第二层换闭合路径形状。把圆形路径换成矩形路径只要仍然环绕同一个BIC一圈结果也必须相同。这个校验能抓出因采样点分布不均导致的采样混叠。第三层和理论预期对照。对于对称保护BIC从对称性分析可以确定某个带的BIC的拓扑荷符号对于偶然BIC文献报道过很多典型结构的拓扑荷。先用已知结构标定自己的代码再跑新结构而不是一上来就处理完全未知的模型这能省下大量的调试时间。6.3 适合新手的最小复现路径如果你是第一次做这项仿真计算我强烈建议不要直接在完整的三维超表面模型上硬刚。先做一个更小的复现实验把流程跑通再说。最小复现路径我推荐这样使用第3节给的硅圆柱正方晶格只算Γ点的对称保护BIC。围绕Γ点取半径R0.02(2π/a)的圆采样32个点用频域端口法提取Jones矢量并计算拓扑荷。由于对称保护BIC的理论拓扑荷在大多数约定下是1少数文献给−1取决于坐标约定如果你的流程输出恰好是±1说明整条链路是通的。然后你再换非零动量的偶然BIC或者换成更复杂的结构此时所有调试精力都可以放在物理本身而不是数据处理流程。我自己的调试习惯还有一个先在一条直线路径上测试。固定kyBIC的ky把kx从BIC左侧扫到右侧画出每个点的ψ看它是不是连续变化。直线路径虽然不是闭合回路但能从视觉上快速暴露相位跳变、采样不够、端口参考面问题等早期异常。确认直线路径没问题后再升级到闭合圆周路径这样出问题时能快速定位到底在哪一步。