光子晶体板BIC的动量空间偏振拓扑：从建模到拓扑电荷提取

这段时间一直在琢磨光子晶体板里的连续谱束缚态BIC恰好读到PhysRevLett.124系列里关于动量空间偏振拓扑的那几篇工作自己也动手算了远场偏振态。这东西有个很有意思的点BIC本身是“看不见”的它不向外辐射Q值趋近无穷但只要你绕开它的动量位置画一个闭合回路远场偏振方向就会像绕涡旋一样整体转一圈拓扑电荷妥妥地写在那。这篇文章不打算复述论文而是把我从建模、仿真到后处理这条线完整捋一遍重点放在远场偏振怎么算、拓扑电荷怎么抽、以及计算过程中哪些地方容易翻车。适合正在复现BIC拓扑特性的研究生也适合刚接触动量空间偏振分析、想搞清楚“绕数”到底怎么从数据里提取出来的朋友。1. 在动笔算远场之前先用对称性摸清BIC的家底很多人一上来就铺网格跑大仿真结果算出来的“BIC”其实是数值假象。我自己的经验是无论你是想复现文献里的涡旋图案还是想在新结构里找BIC第一步都应该是对称性分析而不是数值计算。1.1 BIC在动量空间里为什么“隐形”先理清机制。光子晶体板的辐射通道对应倒格矢和衍射级次布洛赫模式只要满足相位匹配条件原则上都能耦合到外部平面波产生辐射损耗。BIC特殊就特殊在在某个动量点 ( \mathbf{k}_0 ) 上虽然模式频率处在辐射连续谱里但辐射通道的耦合系数恰好归零于是模式变得完全局域化。这个“归零”不是碰运气通常背后有对称性保护。常见的情况是模式在某个镜面对称操作下具有奇宇称而辐射平面波在该操作下是偶宇称两个场的空间重叠积分为零辐射就被禁戒了。这也是为什么绝大多数对称保护BIC会出现在高对称点比如 ( \Gamma ) 点或者 ( M ) 点。动量空间里BIC是一个孤立奇点你在能带图上看到的是一个嵌在连续谱中的离散态频率上很难和普通模式区分但Q值峰会非常尖锐。1.2 偏振态为什么要远场看近场里有大量倏逝波分量和局域场的复杂相位信息直接分析近场电场矢量很容易被表面等离激元、材料吸收或者网格数值噪声干扰。远场不一样远场只留下可以传播的辐射分量每个衍射级次对应一束方向明确的平面波电矢量可以分解成s偏振和p偏振两个正交分量。这样每一束辐射波都对应一个明确定义的偏振椭圆即琼斯矢量的实部取向、椭圆率或者斯托克斯参数。动量空间里每个 ( \mathbf{k} ) 点对应不同的辐射方向也就对应一组远场偏振态。把这些偏振态画在动量平面( k_x, k_y )上就会形成一幅“偏振地图”。BIC的位置恰好就是这幅地图里偏振方向发生涡旋缠绕的中心点。1.3 PhysRevLett.124系列工作给出的计算坐标系那几篇PRL文章包括研究动量空间偏振奇点和拓扑电荷与能带拓扑对应关系的工作给我的启发不只是在物理图像上关键是它们把计算坐标约定讲得很清楚。比如远场偏振方位角 ( \chi(\mathbf{k}) \frac{1}{2}\arg(S_1 i S_2) ) 的定义、s/p分量的基矢选取、以及拓扑电荷 ( q \frac{1}{2\pi} \oint d\mathbf{k} \cdot \nabla_\mathbf{k} \chi(\mathbf{k}) ) 的正负号约定。不同文献对坐标系的定义可能相反这直接导致拓扑电荷符号的差异。你如果不先把坐标约定定死后面算出来的电荷是 1 还是 -1 完全没法跟文献比对。2. 结构建模与仿真参数哪些量直接决定偏振拓扑图对称性分析只能告诉你哪里可能有BIC真要算远场偏振还是得落到具体结构和数值工具上。我以一块均匀介质光子晶体板为例把参数选择和背后的理由一并写清楚。2.1 晶格、孔形、材料厚度怎么选模型主体是一块悬浮在空气中的介质平板周期性地打孔形成二维光子晶体。这是研究BIC最经典的平台文献里大量使用硅、氮化硅、砷化镓等材料。晶格类型三角晶格或四方晶格都可以。三角晶格在 ( \Gamma ) 点通常有六重对称性的模式四方晶格则是四重对称。不同的转动对称性决定了高对称点上允许存在的模式分化和不可约表示。孔形圆孔最简单但圆孔在四方晶格下可能出现某些方向上的模式简并。如果想打开简并、让BIC的涡旋结构更“干净”可以考虑椭圆孔或矩形孔通过打破某个镜面对称来调控辐射通道。平板厚度厚度决定了波导模式的截止频率和色散。通常在 ( \lambda \sim 900 ) nm 附近工作的氮化硅板厚度取 300~400 nm 比较合适。厚度太薄模式被推到光锥外太厚高阶模出现给远场偏振分析带来额外分支。占空比孔半径与晶格常数之比 ( r/a ) 在 0.25~0.35 之间。占空比太小能带比较平BIC附近模式密度低占空比太大孔壁太薄制造和数值建模都容易出问题。作为参考我常用的一组参数是晶格常数 ( a 700 ) nm氮化硅折射率 2.02孔半径 ( r 0.3a )板厚 ( t 350 ) nm上下包层均为空气。这套参数下目标频段内通常能同时看到 ( \Gamma ) 点的对称保护BIC和一个可调位置的偶然BIC。2.2 动量空间网格密度别让涡旋从网格缝里溜走远场偏振计算最核心的数值输入是动量空间中的密集采样。BIC附近的偏振方向变化非常剧烈离奇点越近偏振方位角在单位弧长上的变化越大。如果动量网格太粗糙绕一圈的采样点不够展开后的相位可能直接跳过 ±π 分支拓扑电荷对不上。我的做法是分两步走粗扫找BIC位置在全布里渊区用均匀网格步长 ( \Delta k 0.01\times(2\pi/a) )计算Q值热图锁定高Q峰。细扫围绕BIC以BIC位置为圆心半径取 ( 0.02 \sim 0.1 \times(2\pi/a) )沿角度方向均匀取 128~256 个点径向取 20~30 层。这样单次细扫就有几千个 ( \mathbf{k} ) 点的仿真任务量。好在每个点是独立的本征问题并行化非常方便。2.3 警惕“伪BIC”如何用对称性检验计算结果数值仿真可能出现两类“伪BIC”一类是有限尺寸效应引起的光子晶体板横向尺寸不够大时模式泄漏到边界Q值被截断看起来像很窄的共振但不是真正零辐射另一类是数值损耗带来的仿真软件引入的人工吸收导致所有模式都有有限Q某个模式刚好Q值异常高容易被误判。判别方法很简单检查模式场的对称性。真正由对称性保护的BIC其电场或磁场分量应当具有严格的奇偶性。你可以在模式场输出后手动计算它在镜面操作下的变换看投影到辐射平面波基矢上的积分是否为零。这个检查在建模阶段做一次能省掉后面大量的无效计算。3. 远场偏振的四步计算流程这节是全文的核心操作部分。我按照自己实际跑通的流程按步骤写下来每一步的目的、输入、输出都会说明。3.1 第一步能带求解与模式辨识计算每一个动量点 ( \mathbf{k} ) 上的本征频率和本征场分布。我习惯用基于平面波展开的严格耦合波分析RCWA因为它在处理周期性结构时的数值收敛性有保证而且可以直接给出衍射级次的振幅这是后面远场投影的基础。求解时注意开放系统是有辐射损耗的本征频率是复数虚部对应辐射衰减率。为了准确得到Q值和辐射通道信息必须求解复数本征值问题而不是把周期边界条件当成普通光子晶体能带那样只算实数频率。很多初学者用FDTD跑单个光源激发再通过场衰减拟合出Q值这种方式算BIC可以但拿不到各个衍射级次的复振幅做不了远场偏振分析。想算偏振拓扑尽量用频域方法RCWA或者有限元本征模都可以。3.2 第二步辐射场投影得到本征场后下一步是把模式场分解到可传播的衍射级次上。假设模式场的横向分量是 ( \mathbf{E}{\parallel}(\mathbf{r}) )按照布洛赫定理拆成各阶倒格矢分量。对每个倒格矢 ( \mathbf{G} )辐射波矢 ( \mathbf{k}{\parallel} \mathbf{k} \mathbf{G} )。若 ( |\mathbf{k}_{\parallel}| k_0 )真空波数这个分量落在光锥内可以辐射否则是倏逝波不贡献远场。辐射场的电场矢量可以分解为[ \mathbf{E}{\text{rad}} \propto \sum{\mathbf{G}} c_{\mathbf{G}} \hat{\mathbf{e}}_{\mathbf{G}} ]其中 ( c_{\mathbf{G}} ) 是从本征场投影而来的复系数( \hat{\mathbf{e}}{\mathbf{G}} ) 是相应衍射级次的偏振基矢。具体分解时把基矢取为垂直于 ( \mathbf{k}{\parallel} ) 的两正交方向一个定义为s偏振垂直于入射面另一个是p偏振平行于入射面。3.3 第三步构造动量空间的琼斯矢量与偏振椭圆参数对每个 ( \mathbf{k} ) 点把两个正交偏振分量的复振幅组合成琼斯矢量[ \mathbf{J}(\mathbf{k}) \begin{pmatrix} E_s(\mathbf{k}) \ E_p(\mathbf{k}) \end{pmatrix} ]接下来可以根据琼斯矢量计算斯托克斯参数。我自己最常用的三个量是( S_0 |E_s|^2 |E_p|^2 )总强度( S_1 |E_s|^2 - |E_p|^2 )( S_2 2\operatorname{Re}(E_s E_p^*) )注意这里不需要 ( S_3 )。因为从二维动量平面辐射到远场的电场矢量本质上是一个实矢量在垂直于传播方向的平面内的振动椭圆率信息固然可以从 ( S_3 ) 得到但拓扑电荷关心的是偏振主轴在动量空间的旋转也就是方位角 ( \chi )它由 ( S_1, S_2 ) 决定[ \tan 2\chi \frac{S_2}{S_1} ]之所以用这个式子而不是直接算 ( \arctan(E_p/E_s) )是因为 ( 2\chi ) 的周期是 ( \pi )不存在 ( \pi/2 ) 角度突变相位展开更稳定。3.4 第四步绘制绕BIC闭合回路的偏振演化拓扑电荷本质上是闭合积分的结果所以最后一步是选取一条围绕BIC的闭合路径通常取以BIC动量坐标为中心的圆半径在细扫范围内。沿圆周每个采样点读出 ( \chi(\theta) )这里 ( \theta ) 是圆周上的极角。正常情况下( \chi ) 会随着 ( \theta ) 单调变化绕一整圈变化量为 ( 2\pi q )。把 ( \chi(\theta) ) 画出来你会看到一条斜率为 ( q ) 的直线如果 ( \chi ) 展开成连续函数或者一条来回折叠的锯齿线如果做了模 ( \pi ) 操作。我建议把原始数据和展开后的数据分开存原始数据用于检查偏振分布是否存在断点展开后的数据用于计算斜率。直径方向上出现斜率突变、锯齿跳变不规则的地方往往是数值采样不足或者模式交叉区域的信号。4. 如何从远场数据里稳定抽取拓扑电荷拓扑电荷的计算说穿了是绕数计算公式不复杂但数值实现里到处是坑。4.1 幅角积分与相位展开的工程细节严格做法是[ q \frac{1}{2\pi} \oint_C d\mathbf{k} \cdot \nabla_{\mathbf{k}} \chi(\mathbf{k}) ]圆路径上离散化为[ q \approx \frac{1}{2\pi} \sum_{i1}^{N} [\chi(\theta_{i1}) - \chi(\theta_i)] ]其中每一步的角度差需要做 unwrap即限制在 ( (-\pi/2, \pi/2] ) 区间内。因为偏振方位角是模 ( \pi ) 的物理量不是模 ( 2\pi ) 的所以计算差值时看的是 ( \pi ) 域上的跳变。具体实现上我推荐先用 ( \tan 2\chi S_2/S_1 ) 得到每个点的 ( \chi )然后对 ( 2\chi ) 做标准的 ( 2\pi ) 展开最后再除以2。这样能避免 ( \chi ) 本身跨 ( \pm\pi/2 ) 边界的跳变问题。4.2 符号约定为什么你算出的电荷和文献相反我在复现过程中遇到一个非常典型的诡异结果围绕同一个BIC用论文里的公式算出来拓扑电荷是 1用我自己的脚本算却是 -1而且两边看起来都“自洽”。根源在偏振基矢的定义顺序。有人把s偏振放在基矢第一分量有人放p偏振有人定义方位角是相对于 ( k_x ) 轴的夹角有人相对于 ( k_y ) 轴。你交换 ( s/p ) 两个分量的顺序或者把x/y轴互换一次( \chi ) 会变成 ( -\chi ) 或 ( \pi/2 - \chi )绕数符号就可能反转。我的建议是在你自己的代码里固定一套约定并写清楚。比如基矢顺序固定为 ( (E_s, E_p) )方位角 ( \chi ) 定义为 ( E_s/E_p ) 复比值相位的一半动量坐标系与晶格 ( a_1, a_2 ) 方向重合观察方向取从元器件上方看向动量平面的方向。然后拿一个已知拓扑电荷的文献案例比如 ( \Gamma ) 点的对称保护BIC通常报道为 ( q 1 )做校准。校准通过后再算你的新结构。4.3 多模式叠加区的拓扑含义在动量空间的高对称线上两个模式的偏振分量可能会交叉导致某个 ( \mathbf{k} ) 点上 ( S_0 0 )也就是完全消偏振。从拓扑角度讲这些点是偏振矢量场上额外的奇点可能是涡旋源也可能是鞍点。在抽取BIC的拓扑电荷时闭合路径不能穿过这些点否则幅角跳跃无法处理。实际操作中我会在细扫阶段顺便输出各模式的重叠积分。如果闭合路径附近有其他模式与目标BIC模式的频率差小于线宽就缩小回路半径或换一个对称方向避开模式交叉区域。5. 我踩过的数值坑远场相位在奇点附近是怎么崩掉的这节写几个真实遇到并最终解决的问题比前面的理论更有直接参考价值。5.1 高Q本征模式的收敛阈值BIC附近模式Q值极高辐射损耗极小这意味着本征频率虚部非常小。RCWA或者有限元的迭代求解器在这种参数下很容易收敛到一个“伪模式”上——频率实部对虚部完全错误Q值被低估几个数量级。我踩过最严重的一次三层网格加密后Q值从 ( 10^5 ) 跳到 ( 10^7 )再加密又跳回 ( 10^4 )完全不稳定。最后查下来是倒格矢截断阶数不够高角向谐波分量的倏逝场没被收敛。对于氮化硅这种折射率约2的材料倒格矢截断至少取 ( N_G \geq 21 )即每个方向 ±10 阶才比较稳折射率再高的材料要更多阶。5.2 极化率与数值噪声的伪装电荷另一个大坑是在远离BIC的动量区域远场强度很弱这时候数值噪声会影响偏振方向的判断。你会看到( \chi ) 场里出现很多成对的小型涡旋电荷正负抵消看起来像是物理上有趣的结构。判断方法是检查远场强度 ( S_0 )。真实奇点处 ( S_0 0 )数值噪声引起的赝涡旋通常出现在 ( S_0 ) 很小但不完全为零的位置。把 ( S_0 ) 当作权重来画 ( \chi )噪声涡旋会自动淡出物理涡旋则保持清晰。5.3 中心BIC点本身怎么处理( \Gamma ) 点的对称保护BIC对应 ( \mathbf{k} 0 ) 附近 ( E_s, E_p ) 同时趋于零这个点本身不该参与差值计算。我会在闭合路径上排除距离奇点过近的采样点通常要求采样半径大于三个网格间距防止数值奇异值污染积分。此外如果结构不是完美镜面对称例如基板或衬底引入微小折射率梯度BIC会从 ( \Gamma ) 点偏移到附近某个 ( \mathbf{k}_0 )闭合回路要相应地以 ( \mathbf{k}_0 ) 为圆心。不要拿 ( \Gamma ) 点当BIC位置否则附近偏振图景是缺了一块的绕数必然算错。6. 结果的可信度检验与后续扩展算完远场偏振拓扑电荷不能直接写进报告里。我总结了三个自检维度方便你判断结果是物理可靠还是数值幻觉。6.1 三个自检判据第一拓扑守恒。对理想无损耗平板整个布里渊区内所有涡旋的拓扑电荷之和应为零净拓扑荷守恒。我算过三角晶格光子晶体板在 ( \Gamma ) 点有一个 ( q1 ) 的对称保护BIC必然在布里渊区边界或者某个高对称线上对应有一个 ( q-1 ) 的伴随奇异点。只看到单个电荷而不找它的配对意味着你采样的动量范围窄了。第二Q值幂率。把BIC附近沿某个方向偏移 ( \delta k )Q值应当按 ( Q \propto (\delta k)^{-2} ) 发散。收敛的数值结果应该能拟合出指数 -2 附近的值。如果拟合指数偏离很多说明模式不是真正的BIC而是被数值损耗或者结构不对称性破坏掉的准BIC拓扑电荷的物理意义会变质。第三晶格对称性一致性。四方晶格下( \Gamma ) 点的BIC周围偏振场应该表现出与晶格群一致的四重旋转对称性。计算结果如果出现明显不对称第一怀疑网格没有沿对称方向对齐第二怀疑倒格矢截断破坏了模式简并方向。6.2 从偏振地图还能读出什么Q-BIC与偏振分束做这一步不只是为了验证BIC存在。动量空间的偏振拓扑结构可以直接用于设计偏振分束器件某个 ( \mathbf{k} ) 方向对应纯s偏振另一个方向对应纯p偏振入射光按动量分开并被不同衍射级次携带就是一个不依赖传播距离的空间偏振分光方案。另外通过微调结构参数让BIC从完全偶合转向准BIC即人为引入一定辐射损耗可以在保持拓扑结构的同时获得有限Q值这类Q-BIC共振可以用来增强光学非线性、激光或者传感灵敏度而它的偏振响应继承了母体BIC的拓扑性质天然具备对偏振和入射角的敏感选择性。6.3 我的一点个人建议最近几次做这类计算我的体会有三一是把对称性分析当成本地固件看待别迷信仿真软件一键结果二是相位展开的代码一定要单独写好并测试它比本征求解器更容易出隐蔽错误三是养成把原始远场数据、展开后数据、最终拓扑电荷结果分别存档的习惯审稿或者组内讨论时拿得出手也方便不同约定下交叉验证。动量空间里的偏振拓扑就像一副地图BIC是地图上标着台风眼的风暴中心远场偏振方向是绕着中心旋转的流线。只要把动量网格铺得足够密、把相位展开做得足够稳、把坐标约定钉得足够牢那一圈圈涡旋图案就一定会老老实实地在你面前展开。