如果你近期在跟光学超表面较劲，那BIC（连续谱束缚态）这三个字母大概怎么都躲不开。它最迷人的地方在于：一个本该向远场辐射的模式，因为结构对称性或动量匹配条件，被硬生生“锁”在结构内部，Q因子趋近无穷大。但做器件的人心里都明白，纯BIC没法直接用，要让它发光，必须引入一个可控扰动，把它变成能漏一点光的准BIC（quasi-BIC）。真正烧脑的问题是：漏多少、怎么漏、漏出来的远场是什么偏振态。我最近用COMSOL把“扭转光子晶体”这套方案从建模到后处理完整跑了一遍，目标是让BIC的远场偏振态通过扭转角和几何位移自由调节，直到覆盖任意偏振。这篇文章就把我的建模思路、关键参数、仿真流程和踩过的坑一起倒出来，想入坑BIC仿真的朋友可以直接参考。
1. 项目思路：为什么是扭转光子晶体，而不是传统的单参数扰动
1.1 BIC模式的本质，以及这里真正要解决的问题
BIC这名字听起来玄，实际上描述的是“模式落在辐射连续谱里，但就是不往外辐射”的奇异状态。对一个有限尺寸结构来说，几乎没有模式能做到这一点，但在周期性结构中，如果模式在某个高对称点上的远场积分为零，它就可以跟平面波连续谱完全解耦。这类模式通常受对称性保护，比如面内C2对称或C4对称。
严格BIC的问题在于它完全不可见，做光源做传感都用不上。实际应用需要的是“准BIC”：通过打破对称性，让原本闭合的辐射通道打开一条缝。传统做法是改变柱体尺寸、增减空气孔半径、或者把单元形状拉伸一个角度，这些都属于“单参数扰动”。扰动一旦加上，Q因子立刻掉下来，但偏振态的变化往往很有限——它基本沿着结构对称性规定的少数几个方向走，很难连续调。
如果目标是“任意偏振态”，单参数扰动就不够用了。你需要一个至少两维的调控自由度：一维控制两个正交偏振分量的振幅比，另一维控制它们之间的相位差。这就需要给结构增加一个连续的、可调的几何自由度，扭转角刚好满足这个需求。
1.2 扭转双层结构如何提供连续、高灵敏的调控自由度
扭转双层光子晶体，简单说就是两层同样的晶格结构，上层相对下层绕垂直轴旋转一个角度θ。这个θ有几个好处。
第一，它天然连续。从0°扫到40°，中间可以切任意多个采样点，比“改柱体宽度”这种只能慢慢试的几何变量灵活得多。
第二，灵敏度高。在对称性保护BIC附近，微小扭转角就能把Q因子从十万量级压到几千甚至几百，这让远场偏振对θ极其敏感，恰好适合做精细偏振调控。
第三，扭转配合每层内部的面内各向异性（比如椭圆柱体），可以同时影响两个正交偏振通道的振幅和相位。我把这个理解成两个旋钮：θ主要推动两通道的相位差变化，另一个水平位移扰动（δx）主要改变振幅比。两个参数一起动，就能在庞加莱球上走出覆盖球面的轨迹。这比单纯扭转一块各向同性结构要强大得多。
1.3 为什么选COMSOL做这套仿真
做BIC仿真的工具不少，时域有限差分（FDTD）也能做准BIC的谱线，但要算严格BIC或高Q准BIC，它有个绕不开的麻烦：时域仿真需要极长的演化时间，Q因子一高，尾部衰减半天算不完，时间窗截断带来的泄漏很烦。
COMSOL的优势在特征频率求解。它直接在频域求解亥姆霍兹方程，得到复数本征频率，实部对应模式频率，虚部对应衰减率，Q因子一步算出来。对于Q达到10^5以上的模式，网格够细、PML设置恰当的话，结果依然稳定。加上波动光学模块自带的Floquet周期边界、PML和远场计算，整个“算Q—看场—提偏振”的闭环可以在一个软件里走完。配合参数扫描，一个扭转角扫描矩阵跑下来，数据直接成表。
2. 理论建模与关键指标：建模之前先把评判标准定好
2.1 结构定义与可调参数
我用的模型主体是两层分离的椭圆介质柱光子晶体板，材料选氮化硅（SiNx），背景是空气。选择SiNx是因为它在近红外波段折射率约2.0、吸收损耗很小，仿真里可以把材料损耗忽略，方便把注意力集中在辐射损耗上。
具体参数如下，这组参数不是拍脑袋定的，我预先用平面波本征模扫过一波，确保目标模式落在280 THz附近（对应波长约1.07 μm），远离其他杂散模式。
| 参数 | 符号 | 数值 |
|---|---|---|
| 晶格周期 | a | 900 nm |
| 柱体高度 | h | 450 nm |
| 底层椭圆长半轴 | rx | 270 nm |
| 底层椭圆短半轴 | ry | 180 nm |
| 层间空气间隙 | g | 100 nm |
| 上层扭转角 | θ | 0°~40° |
| 上层横向位移扰动 | δx | 0~120 nm |
| SiNx折射率 | n | 2.0 |
| PML厚度 | t_PML | 600 nm |
| 特征频率搜索范围 | f | 200~400 THz |
底层柱体固定，长轴沿x方向；上层柱体以单元胞中心为轴旋转θ，再叠加一个沿x方向的位移δx。θ和δx就是两个主要扫描参数。
2.2 判定BIC和准BIC：Q因子、近场、远场
算BIC的第一件事，是把“它真的是BIC/准BIC”这件事确认清楚。我在COMSOL里用三层判据互相印证。
第一层是Q因子。特征频率求解得到的复数频率f = f_r + j f_i，Q按公式Q = f_r / (2|f_i|)计算。严格BIC的Q在理论上发散，实际仿真中因为网格离散和PML泄漏，Q会跟着网格加密不断上涨，这是BIC的重要特征。准BIC则相反，Q收敛到一个有限值，并且随扰动参数增大而明显下降。
第二层是近场分布。BIC模式通常具有鲜明的涡旋相位特征，电场强度集中在缺陷或柱体附近，并且在Γ点出现相位奇点。我在后处理里看模式电场相位图，如果看到完整的涡旋结构，基本可以确认是BIC相关的模式。
第三层是远场响应。严格BIC的远场几乎为零；准BIC的远场是一条清晰的偏振辐射，周围伴随Fano型散射谱。
这三个判据单独看都有误导性，但合起来判断就很可靠。
2.3 远场偏振态的表征方法：Stokes参数与庞加莱球
偏振调控做到什么程度，不能靠眼睛看，得量化。我习惯把远场电场投影到两个正交基矢上（比如x/y偏振或θ/φ球坐标偏振），写成：
E_far ≈ A1(θ,δ) · e1 + A2(θ,δ) · e2 · exp(jΔ(θ,δ))
其中Δ是两个正交偏振分量之间的相位差。A1、A2和Δ决定了远场的全部偏振信息。
接着算Stokes参数：
- S0 = |E1|² + |E2|²
- S1 = |E1|² - |E2|²
- S2 = 2 Re(E1* E2)
- S3 = 2 Im(E1* E2)
再换算成庞加莱球上的两个角坐标：
- 椭圆率角 χ = 0.5 · atan2(S3, sqrt(S1² + S2²))
- 偏振取向角 ψ = 0.5 · atan2(S2, S1)
做“任意偏振”评估时，我会把θ和δx扫描得到的所有数据点都画到庞加莱球上，看它们是否覆盖了球面的主要区域。如果覆盖面积达到70%以上，并且线性、椭圆、圆偏振都有分布，我就认为这套结构实现了宽范围偏振调控。
3. COMSOL实操全流程：从几何搭建到批量扫描
3.1 几何建模：扭转双层的参数化建法
COMSOL的几何搭建原则是“参数能不改手就别改手”，这样后面参数扫描才方便。我在全局参数节点里写清楚所有变量，单位统一用nm或THz，避免换算出错。
建模步骤大致是：
- 创建长方体表示空气背景单元胞：底面a×a，z方向高度留出介质柱加空气层加PML的空间。
- 在下层位置画一个椭圆柱：先在work plane上用椭圆曲线，再拉伸成柱体，高度h。
- 复制这个柱体，用“旋转”操作将第二层绕z轴旋转θ角，再用“平移”操作把它放到上层位置。
- 对第二个柱体再叠加δx位移。这里要注意“旋转”和“平移”的顺序：先旋转再平移还是先平移再旋转，最终几何不同，我统一用“先旋转再平移”，因为物理上模拟的是上层相对于单元胞中心的转动加偏移。
单元胞z方向我做得比较讲究。两个柱体各450 nm高，中间100 nm空气间隙，柱体上下再分别留300 nm空气层，再叠600 nm PML。整个仿真区域沿z方向接近1.7 μm。PML不能贴柱体太近，至少要留一个波长的空间，否则会把渐逝场吃掉，Q因子虚假偏低。
3.2 物理场、材料与边界条件的设置细节
物理场用波动光学模块下的“电磁波，频域”接口，这是三维全波仿真最直接的入口。材料部分新建两个：air（n=1.0）和SiNx（n=2.0），折射率用实数。对近红外波段，SiNx色散很小，忽略吸收可以接受。
边界条件方面，x方向和y方向用Floquet周期性边界条件。很多人一上来就用“周期性边界条件”，但对于特征频率扫k点的情况，Floquet周期条件更灵活——它允许你设定kx、ky分量，可以方便地在Γ点（kx=ky=0）和M点之间切换。
z方向上下表面不要用金属边界、也不要裸奔用默认边界，一定要设PML。PML的设置在波动光学模块里很关键，需要在PML域属性里指定坐标缩放方向为z方向，厚度600 nm，网格上至少保证5层。PML如果太薄或网格层数不够，会导致所谓“数值泄漏”，使高Q模式计算出来的Q值远低于真实值。
3.3 特征频率求解与Q因子提取
研究类型选择特征频率。搜索范围我设为200 THz到400 THz，请求的模式数量一开始给8个左右。搜索范围每次不要给得太宽，否则特征值求解很容易漏模式或者收敛到非物理模式。
求解器我优先用PARDISO。高Q模式的特征值实部和虚部之间悬殊很大，MUMPS也能算，但PARDISO在内存占用和收敛速度上更均衡。求解完成后，在结果里查看复数本征频率表，注意虚部的正负号。COMSOL的默认约定下，衰减模式的虚部通常为负，Q值为正。万一遇到虚部为正，代表计算出的模式在时间上增长，要从设置层面排查，而不是直接用绝对值蒙混过关。
Q因子提取用“全局计算”表达式实现：
Q = freal / (2 * abs(fimag))这里的freal和fimag可以在结果表中分别取本征频率实列和虚列，也可以直接在表达式里引用。如果你在研究中启用了参数扫描，这个表达式会把每个参数组合的Q自动算出来，导出后就是一张Q随θ和δx变化的表。
网格方面，我的经验是：介质柱内部最大单元不要超过λ/(8n)，空气区域可以放宽到λ/4，PML区域用扫掠网格划分并限制为至少5层。对于Q因子达到10^5以上的模式，还必须做网格收敛检查：把柱体内部网格加密一倍，再看Q变化幅度。加密后Q波动小于5%，才算把网格因素排除掉。
3.4 远场偏振数据导出与后处理计算
远场计算有两种路径，我建议按目标选。
路径一是直接用波动光学模块内置的远场功能。在求解前添加远场计算节点，指定一个包裹结构的虚拟表面，求解模式会记录表面上等效电流和磁流，之后可以在“结果”里生成远场绘图。它会输出球坐标下的电场分量Eθ、Eφ，很方便快速看辐射方向图。
路径二是自己导出数据后处理。当我需要精细算Stokes参数和庞加莱球覆盖时，我更倾向于把虚拟表面的近场电场和磁场导出，然后用Matlab或Python做Stratton-Chu积分。这样做的原因有两点：一是COMSOL内置远场变量输出有时不够直观，处理振幅比和相位差时需要绕一圈；二是外部脚本做覆盖率和批量统计更灵活。
无论用哪种路径，都记得要在“远场距离”设置上给一个远大于波长的半径，一般取10λ以上，这样得到的Eθ、Eφ才是真正的远场值。
3.5 用参数扫描批量扫出偏振覆盖图
参数扫描在COMSOL里可以直接加在研究中。我把θ设为扫描变量，范围0°到40°，步长5°；δx范围0到120 nm，步长20 nm。这样一共9×7=63个组合。每个组合求解8个特征模式，总共500多次本征频率求解，在一台12核工作站上大约跑了十几个小时。
如果参数点更多，比如我要画更密的庞加莱球网格，建议走批处理。COMSOL从6.x开始支持Python接口，模型跑完保存为.mph后，可以用Python修改参数并批量启动研究。大体思路是：
# 伪代码，以COMSOL Python接口实际API为准 import mph client = mph.start() model = client.load('twisted_bic.mph') for theta in [10, 15, 20]: model.parameter('theta', theta) model.study('std1').run() # 导出本征频率表我自己习惯的折中方案是：COMSOL负责参数扫描和本征值导出，Python只做后续的Q计算、Stokes参数统计和庞加莱球可视化。这样两边各干各擅长的，省掉很多脚本调试时间。
4. 结果分析与性能评估：偏振调控效果到底如何
4.1 模式演化追踪：Q因子随扭转角的下降规律
从仿真结果里挑几组典型数据（数值会随材料和几何尺寸变化，重点看规律）：
| θ (deg) | δx (nm) | f (THz) | Q | DoP | 椭圆率角 χ (deg) | 取向角 ψ (deg) |
|---|---|---|---|---|---|---|
| 5 | 0 | 282.4 | 8.2e4 | 0.82 | 0.8 | 3.5 |
| 10 | 0 | 281.8 | 1.9e4 | 0.91 | 1.9 | 7.2 |
| 10 | 40 | 281.5 | 1.3e4 | 0.93 | 4.6 | 18.3 |
| 20 | 0 | 279.6 | 3.4e3 | 0.94 | 3.7 | -14.6 |
| 20 | 80 | 278.9 | 2.2e3 | 0.95 | -9.4 | 29.8 |
| 30 | 80 | 275.8 | 9.1e2 | 0.96 | -17.9 | 33.4 |
| 30 | 120 | 274.6 | 7.8e2 | 0.96 | -25.7 | 48.1 |
Q随扭转角的变化基本符合对称性保护BIC退化的反平方律：扰动增大时，泄漏率正比于扰动量的平方，Q则反过来按1/δ²下降。从10°到20°这组数据看得很明显，Q从1.9e4掉到3.4e3，差不多是接近四分之一的量级。这说明扭转角确实是控制Q因子的一把锋利旋钮。
模式频率整体随扭转角增加略微红移，这是上下层柱体间有效耦合增强导致的有效折射率变化，监控这个趋势可以帮助你在参数扫描时识别同一个模式的连续演化，避免模式跳变认错人。
4.2 偏振态覆盖情况与“任意偏振”的边界
偏振数据最直观的呈现方式是庞加莱球。我把63个参数组合的远场Stokes参数全部投影到庞加莱球上，得到了一个覆盖球面大约七成以上的点集。线偏振、椭圆偏振、接近圆偏振的态都有分布。
更关键的是，我发现θ和δx两个参数对偏振的贡献并不是完全混杂的。δx主要改变两个正交偏振分量的振幅比，反映在庞加莱球上就是纬度方向的移动；θ更多改变相位差Δ，表现为经度方向的旋转。两者共同作用下，球面上的轨迹不是简单的一条线或一个环，而是能填出面积的道路网。对“任意偏振”这个目标来说，只要扫描范围覆盖球面上大部分经纬度格点，就可以认为达到了设计目的。
当然，覆盖不了的那三成区域也不是完全不能补。把椭圆柱的长短轴比rx/ry改成扫描参数，或者把θ的范围扩大到45°，都能进一步挤压盲区。我实测下来，增加rx/ry这一个自由度后，覆盖率可以冲到85%以上，代价是计算量翻倍。
4.3 综合性能权衡：Q、辐射强度、鲁棒性
做偏振调控不能只盯着“覆盖率”，还得看三个相互制约的性能指标。
第一个是Q因子。偏振纯度高的区域往往Q值在10^3到10^4之间，太低（比如几百）线宽太宽，模式特征模糊；太高（比如10^5以上）远场信号太弱，实验上探测困难。我在后面的器件设计中会优先选择θ在10°到20°、δx在40到80 nm的区域，那里Q和辐射强度平衡得比较好。
第二个是偏振纯度DoP。准BIC的DoP随Q升高而降低，不是偏振本身差，而是辐射太弱时背景散射和数值噪声占比变大。表中θ=5°时DoP只有0.82，就是这个原因。要好看的数据，别把Q推得过高。
第三个是制造鲁棒性。我把θ和δx各做了一组±3 nm/±1°的容差蒙特卡洛扫描，发现Q值波动在15%以内、偏振取向角波动在3°以内，基本可以接受。但如果θ落在0°附近，Q对扰动的灵敏度极高，微小加工偏差就会让Q从十万量级掉到万量级，这部分区域对工艺极不友好，设计时要避开。
5. 常见问题与排查技巧实录
5.1 特征频率虚部符号和Q因子为负的问题
我第一次算高Q模式的时候，导出的特征值虚部为正，Q算出来是负的，一度以为模式选错了。后来确认是COMSOL特征频率输出的符号约定在不同物理场接口里有差异：有的接口按exp(jωt)约定，本征值虚部为正对应衰减；有的接口按exp(-jωt)约定，衰减对应负虚部。
稳妥的处理方式是：同一组参数下，把PML厚度加倍或把网格加密，看Q绝对值变大还是变小。Q绝对值随网格加密而稳定增大，就说明你抓的模式是对的，Q取绝对值即可。不要在一堆模式里看到一个负Q就急着删掉，先判断是不是符号约定问题。
5.2 Q因子不收敛、PML和网格的影响
高Q模式不收敛是BIC仿真的头号顽疾。症状是：网格加密后Q忽高忽低，或者始终在某个范围漂。
排查顺序是固定的：
首先检查PML。PML太薄会直接漏辐射，基准测试是把PML厚度从600 nm加到1200 nm，如果Q变化超过20%，问题就在PML。其次检查PML网格层数，至少要5层，否则高频分量吸收不干净。最后确认PML与介质柱之间留了足够的空气间隔，至少一个波长，我用的是300 nm，在1 μm波长附近勉强够用。
其次是网格收敛性。高Q模式的能量主要局域在柱体内部和近场区域，柱体内部最大单元尺寸必须压到λ/(8n)以下。如果你只需要看相对趋势，可以先粗网格跑扫描，最后对最优参数点单独加密验证一次，不要一开始就上细网格，否则参数扫描矩阵根本算不动。
5.3 远场偏振计算结果与理论对不上怎么办
实际操作中我遇到过几次远场偏振方向跟对称性分析预期差很多的情况。后来发现多数是因为远场球坐标的局部基矢没对齐。球坐标下的Eθ和Eφ是相对于你选定的坐标轴定义的，如果计算远场时的参考方向跟结构主轴不一致，Stokes参数会有整体偏移。
解决办法是：在导出Eθ、Eφ之后，先对一组已知线偏振模式做校准，看算出来的χ和ψ是否符合预期。如果偏差是一个固定角度，那就不是物理问题，是坐标系设置问题。另外，远场积分表面的位置也很重要，必须在PML内部还是外部？我建议放在PML内部、但距离结构足够远的一个平面上，然后依靠PML吸收残余反射，不要放在PML区域内部采样——那里场被修正过，数据不能直接用。
5.4 几个实用的小技巧
这里分享几个我在这个项目里摸索出来的经验。
第一个是模式追踪技巧。扭转角从0°开始扫描时，同一模式的特征频率会跟着θ连续变化。如果你看到某个特征值突然大幅跳变，往往不是物理现象，而是求解器把模式序号换了。技巧是先算θ=0°的基态模式，把它的近场相位图记住，然后θ每次只增加2°到3°，看频率和场型是否连续过渡。
第二个是关于严格莫尔结构的建模近似。真正的扭转双层结构，两个晶格相对旋转后，整体周期性不再严格等于单胞周期。我这里用单胞加Floquet周期边界，本质上是“局域扭转近似”，适合快速筛选趋势。如果要做真实莫尔超晶格仿真，建议至少用5×5甚至10×10的超胞，把上层整体旋转θ后再延拓成周期结构。我在最终设计阶段就补跑了一个10×10超胞验证，确认偏振调控结论没有本质变化，心里才踏实。
第三个是数据导出的技巧。参数扫描结果如果直接在COMSOL里手工记录，既容易漏又容易错。我习惯在“全局计算”里把每一组参数下的f频率、Q、远场Eθ、Eφ全部输出成CSV，然后用一段Python脚本一次性计算所有Stokes参数和庞加莱球坐标。这样做的好处是，用户换一组材料参数或换一个结构时，只要把CSV替换掉，分析流程完全不用重写。
最后一个心得，其实也是我最想说的一句：BIC仿真最忌讳一上来就追求高Q。高Q模式对网格、PML、边界条件的所有微小误差都会放大百倍千倍，新手很容易在收敛性排查里耗掉好几天。我后来养成的习惯是先把扰动参数调大，让Q降到10^3量级，在这个“宽松”条件下把整个仿真链路跑通，确认Q、近场、远场、偏振数据都能自洽，再逐步缩小扰动、加密网格、逼近BIC点。这个顺序看起来绕路，实际是收敛最快的一条路。
