这段时间一直在琢磨光子晶体板里的连续谱束缚态（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_{i=1}^{N} [\chi(\theta_{i+1}) - \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 ) 点有一个 ( q=+1 ) 的对称保护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是地图上标着台风眼的风暴中心，远场偏振方向是绕着中心旋转的流线。只要把动量网格铺得足够密、把相位展开做得足够稳、把坐标约定钉得足够牢，那一圈圈涡旋图案就一定会老老实实地在你面前展开。
