☰
COMSOL双层扭转光子晶体准BIC仿真:偏振调控与远场提取全流程
2026/10/5 4:17:12 网站建设 项目流程

前阵子我一直在做双层扭转光子晶体里的准BIC仿真,想着把"扭转角"变成一个可调的偏振旋钮。翻了一圈社区,发现单层光子晶体里算BIC的教程不少,但一加上"双层相对旋转"这个自由度,很多讨论就直接卡在几何搭建上:两层之间怎么相对转?小角度带来的莫尔周期怎么进COMSOL?BIC的远场偏振又该怎么从本征模里抠出来?这些问题不解决,后面谈远场调控都是空中楼阁。

这篇文章就把我自己跑通这条线的完整过程整理出来,覆盖物理图像、COMSOL建模、特征频率求解、远场偏振提取和性能分析。适合正在做连续体束缚态(BIC)、莫尔光子晶体、偏振调控仿真,或者想用COMSOL复现光学类顶刊工作的读者。先说明一下,我用的软件版本是 COMSOL 6.4,借助波动光学模块的三维电磁波频域接口完成全部计算,后面提到的所有操作流程在6.3上也基本适用,只是求解器细节略有差异。

1. 扭转双层板里的BIC是怎么"漏"出来的:先建立正确的物理图像

1.1 连续谱束缚态:表面上是"不辐射",本质是"对称性禁戒"

很多刚接触BIC的同学会把它理解成"一个完美不向外辐射的光学模式",这个说法方向对,但容易误导。严格意义上的BIC确实是一个本征频率落在辐射连续谱之内、但能量依然局域在结构内部、不向远方泄漏的模式。问题是:连续谱那么宽,为什么光会愿意待在一个有限大小的结构里?原因通常不是"漏不出去",而是"对称性不允许它漏"。

以光子晶体板为例,一个模式要向外辐射,它携带的动量必须能匹配自由空间的传播波,同时它的场分布还要和外部平面波有非零耦合系数。如果模式在面内具有某种对称性,比如对中心点的180度旋转反对称,而外部平面波在这个方向上是对称的,那两者投影为零,辐射通道就被关闭了。这个关闭过程不需要任何材料吸收,纯粹是几何对称性带来的干涉相消。只要加工精度无限好、结构无限理想,这个模式的Q因子就是无穷大。

但真实世界没有无限理想的结构。当我们引入微扰,比如改变孔的尺寸、移动孔的位置,或者把上下两层板相对转一个角度,原本被对称性禁止的辐射通道就会被打开一条缝。这时候原来的BIC变成了准BIC(quasi-BIC),Q因子从无穷大掉到一个有限但可能仍然很大的值,同时向外辐射出携带特定振幅和相位的电磁波。这个"泄露出什么偏振的光"就取决于微扰是怎么破坏对称性的。看到这你应该明白了:微扰的结构形式,决定了远场偏振;微扰的强度,决定了Q因子和辐射强度。这为偏振调控提供了清晰的物理抓手。

1.2 扭转角把对称性破坏变成连续可调过程

如果只是要破坏对称性产生准BIC,其实有很多办法,比如把圆孔改成椭圆孔、把某个孔往旁边移动几十纳米。那为什么偏偏要用"扭转双层光子晶体"这个方案?因为扭转角是一个非常干净、连续可调的几何自由度,而且它带来的物理效果比"单纯打破C2对称"更丰富。

想象两片相同的方形晶格光子晶体板,沿z方向叠放,中间隔一层空气。当旋转角θ=0时,双层结构保持很好的面内对称性,Γ点附近可以存在对称性保护的BIC。当θ从零逐渐增大,两层板之间的相对位置在莫尔超晶格的意义下不再是简单的周期平移,而是产生了一个长周期的面内调制。这个调制等效于给原始模式加上了一个可调的动量耦合,把原本闭合的等频面扭曲,也让模式在倒空间里的辐射通道出现非对称的泄漏。

更关键的是,扭转角不是一个"有/无"的二值参数,而是一个可以连续变化的角度。这意味着我们可以连续地调制辐射通道的开口大小和相位关系,而不是只能选"对称结构"或"某种固定非对称结构"。在实际仿真里,这意味着只要建立好一套参数化几何模型,扫描θ从0度到十几度,就能看到Q因子连续下降、远场偏振连续演变的过程。这种连续性对实验和器件设计都很有价值,因为扭转角在加工上是一个比较好控制的量。

当然,扭转角也有让仿真头疼的一面:它破坏了单晶胞的严格平移对称性。严格说,只有双层晶格矢量的某个公倍数周期才是结构的真实周期,这个周期往往大得离谱。后面我会单独讲怎么在COMSOL里处理这个"周期失配"问题,这也是很多模型卡壳的根源。

2. COMSOL建模第一步:晶胞参数、材料参数和扭转角的周期匹配

2.1 基础参数怎么定:从周期800nm出发的一套常用组合

做光学BIC仿真,第一步是确定一个可以出结果的基础构型。我参考了近红外波段硅基光子晶体里常见的设计,使用如下参数作为起点:

参数符号数值说明
晶格周期a800 nm正方形晶格,格点间距
单层板厚t220 nm硅膜厚度,保证平板导模在目标波段附近
层间空气间隙h_gap120 nm双层板之间的垂直间距
孔半径r140 nm方形晶格上刻蚀圆柱空气孔
硅折射率n_si3.48无损耗介质近似,忽略吸收
背景折射率n_air1.0空气
扫描频率范围f180–260 THz对应约1150–1660nm波段

这个参数组合的好处是:膜厚和周期决定的导模频率落在近红外通信波段附近,便于后续和实验对接;空气孔半径适中,网格不会因为孔壁太薄而爆炸;双层之间120nm的间隙既让两层保持光学耦合,又不会强到把整个模式压缩成纯腔模。

材料方面,我在初算阶段把硅设成纯实数折射率,不引入材料吸收。这样做是有意的:BIC准BIC的Q因子非常高,在无吸收模型里可以单独评估辐射泄漏对Q的贡献。如果一上来就加入材料损耗,虚部会被材料吸收主导,反而看不清几何调控的效果。等到研究工艺容差时再按需打开吸收项。

2.2 扭转角不是直接在边界条件里转一个角度

这是很多初学者最容易掉进去的坑。有人觉得,反正只是两层板相对旋转,那我就在COMSOL里把第二层几何直接旋转θ,然后晶胞四个侧面依然用周期性条件不就行了?如果θ等于零,这么干没问题;但一旦θ不为零,第二层晶格矢量和第一层晶格矢量不再对齐,同一套周期性边界条件不可能同时满足两个晶格的平移对称性。强行算出来的"BIC",本质上是把非周期结构强加了一个虚构周期,得到的是伪模式。

正确的处理思路,是把结构看成两层晶格形成的莫尔超晶格,计算域取超晶格的最小重复单元。理想情况下,小角度扭转的莫尔周期差不多是 a / (2 sin(θ/2)),θ=1°时周期能达到45微米量级,直接三维建模根本不可行。实际工程里要避免这种全莫尔建模。我采用的办法是选一组"有理近似转角",让两层晶格在一个可接受大小的超胞里重新对齐。

举个例子:如果选择第一层晶格矢量为 (a,0) 和 (0,a),第二层相对旋转 θ,那么只要 θ 满足旋转矩阵四个角度的三角函数都是整数比,就可以在一个有限超胞内闭合周期。常见的选择是一层旋转一种类似 2×2 超胞的转角,比如 θ=90°退化为正交错位,但这个角度太极端,物理上不是准BIC的理想区间。另一个更实用的路线是做"近似扭转":取一个超胞,第一层和第二层的孔洞阵列分别按各自的晶格画好,四周用周期性条件,但在几何里保留微小角度失配,只要超胞尺寸适当,误差可以控制在可接受范围。

我在实际模型里更多是绕开这个矛盾,只计算有限的几个代表性转角,比如 θ=2°、5°、10°。做法是:先在第一层几何里画好孔阵列并旋转整体部件,再把第二层以相同原点旋转同样角度但不改变周期方向,而是让两层板都保留自己的严格周期,再把整个超胞视为一个新的周期单元。这种方法对真实莫尔结构是一种近似,但用来研究"扭转角如何调制Q因子和偏振"是够用的。重点在于:在论文中描述模型时必须写明你的扭转周期假设,否则审稿人一定会问超胞选择依据。

2.3 计算域配置:PML、周期性条件和网格收敛

几何建好之后,计算域按如下方式布置:z方向从下往上依次是下层空气缓冲层(厚度约一个波长)、下层硅板、空气间隙、上层硅板、上层空气缓冲层,最外侧再放完美匹配层(PML)吸收泄漏辐射。横向四个侧面使用周期性条件,k点固定在Γ点(kx=ky=0)。研究BIC最先看Γ点,因为许多对称性保护的BIC都在这个高对称点附近出现,而且Γ点附近远场提取比较直观。

PML厚度至少要达到工作波长的1到2倍。我试过只加薄薄一层PML,结果计算出来的模式虚部被PML自身反射污染,Q因子上限被压制在10^4量级,怎么调结构都上不去。后来把PML厚度加到2微米以上,并设置成"从物理域向外拉伸"的坐标变换类型,Q因子上限才恢复到10^6以上。

网格方面,硅板及其周围的场是重中之重。我的收敛标准设置为:连续两次加密网格之后,目标模式实部频率变化小于0.02%,虚部变化小于5%。实际操作中,先在硅板厚度方向用5层扫掠网格,平板面内最大单元尺寸取 λ/(6n_si) 左右;空气孔边界加一层边界层网格;空气间隙区域的网格尺寸可以放宽到硅板内部的两倍。这样一套网格在COMSOL 6.4下,单次特征频率求解大概需要几十GB内存,视超胞大小而定。如果资源紧张,可以先算零扭转角双层结构来调节网格,再把网格迁移到扭转模型上。

3. 特征频率扫描与远场偏振提取:如何把本征模变成偏振数据

3.1 找BIC的实用判据:实部频率和虚部Q因子

在COMSOL波动光学模块里,直接计算准BIC特征的入口是"电磁波,频域"接口下的特征频率研究。求解器输出的本征频率是一个复数,实部对应共振频率,虚部对应辐射损耗。模式Q因子的表达式是:

  • Q = Re(f) / (-2 * Im(f))

这么算出来的Q包含了所有辐射和吸收通道。当材料无损耗时,Im(f)越小,Q越高,说明模式的能量泄漏越弱。对于理想BIC,Im(f)在数学上等于零,COMSOL里只会得到一个受数值精度限制的极小虚部;对于准BIC,Im(f)是一个明显的非零值,Q在10^3到10^6之间浮动。

我做参数扫描时的习惯是:先固定所有几何参数,让特征频率研究搜索一个较宽的范围,把Γ点附近的模式全部列出来;然后按照Q因子从高到低排序,挑出虚部最小的几个模式观察场图。如果某个模式的电场能量高度集中在两层硅板内部,且空气域里的向上和向下辐射几乎为零,这就是值得追踪的准BIC候选。要特别注意模式追踪问题:参数扫描过程中模式会交叉,单纯按序号取解会出错。我一般会给求解器一个包含模式实频的初值范围,并配合COMSOL 6.4的"按上一步解作为初值"选项,让同一个物理模式在连续参数下被稳定跟踪。

3.2 从近场到远场:傅里叶分解与COMSOL远场域两种路径

本征模式自带的是近场分布,远场偏振需要额外提取。对周期性平板结构,我倒推荐用傅里叶分解的思路,比直接插远场域更直观。过程如下:在紧贴上层板的上表面提取切向电场分量 Ex(x,y) 和 Ey(x,y),做二维空间傅里叶变换,得到不同倒格矢分量对应的幅值和相位。在Γ点,辐射连续谱里只有零阶衍射通道(即直上直下方向)能传播到远场,所以取这个零阶分量的 Ex 和 Ey,就得到了远场偏振的两个复振幅。

在COMSOL里实现这个操作,可以使用"派生值"功能里的表面积分配合相位因子。更通用的方式是把表面电场数据通过LiveLink for MATLAB或Python脚本导出,在外部做FFT并换算Stokes参数。如果你就是想全部留在COMSOL界面里,那可以选择在空气缓冲层外层放置一个半球形"远场域",软件会按散射场公式把边界电场投影到远场。不过对周期结构本征模,这个选项容易受到边界截断和缓冲层厚度的影响,所以我个人更喜欢"表达式导出+外部FFT"的路径。

判断远场偏振的核心是复振幅比 Ex/Ey 的幅值比和相位差。举个例子:如果 Ex/Ey 的相位差为0或π,远场是线偏振;如果相位差接近±90度且幅值比接近1,远场接近圆偏振;其他情况是椭圆偏振。把这些比值画到庞加莱球上,就能直观看到参数扫描如何移动模式的偏振态。

3.3 第二个自由度:扭转角之外还需要一个参数才能覆盖庞加莱球

标题里说的"任意偏振态BIC",不能靠扭转角一个参数实现。单个自由度的参数扫描只能让偏振态在庞加莱球上画出一条一维轨迹,想要覆盖球面,至少需要两个相互独立的几何自由度:一个负责改变辐射通道的振幅分配,另一个负责改变两个正交偏振分量的相位差。我在模型里选择的是"扭转角+孔对的横向偏移"组合:扭转角主要调节模式自旋相关的耦合强度,孔对偏移则调节面内对称性的破缺程度。

实际操作时,我会把孔对偏移量 d 也设成全局参数,和扭转角一起做二维参数化扫描。d从0变化到20nm,θ从2°变化到12°。在两组参数交叉扫描下,观察到的远场偏振态能够覆盖庞加莱球的多个象限,而不是局限在一条赤道线上。这个结果说明,通过双参数设计确实可以实现"按需选取偏振态"的效果。对实验来说,扭转角和孔偏移都是可在版图和工艺里实现的几何量,不会增加太多制造难度。

4. 性能分析怎么量化:Q因子、偏振纯度和工艺容差

4.1 Q因子随扭转角的变化规律

在我使用的硅基双层孔洞结构里,Q因子随扭转角的变化呈现一个基本趋势:θ=0附近的严格对称构型对应最高的Q因子,但此时远场辐射接近零,远场偏振无从谈起;θ增大后,辐射通道打开,Q因子下降几个数量级,同时远场信号增强。具体数值依赖超胞近似方式,但趋势是稳定的。

更值得关注的是Q因子下降的速率。小角度时,Q与θ近似满足反平方关系,即θ从2°增加到4°,Q可能下降约4倍。这说明扭转角是非常灵敏的调控手段。如果你想要一个Q=10^5以上的高Q准BIC,角度最好控制在3°以内;如果更看重远场辐射强度和偏振信号的清晰度,角度可以放到5°以上。这个取舍没有标准答案,取决于器件的应用场景。

此外,层间间隙 h_gap 也强烈影响Q因子。间隙增大,两层耦合减弱,BIC可以从"双层共同保护"的模式退化到"单层独立模式",Q因子的变化曲线会出现阶梯状跳变。扫描h_gap时不要只算一个值,建议在100nm到200nm之间多取几个点,看Q曲线是否平滑。如果出现突变,先检查是不是模式追踪跳到了别的模式。

4.2 偏振纯度与远场偏振椭圆

偏振纯度可以定量表示为远场中目标偏振态所占比重。线性偏振纯度常用两正交分量的功率比描述,圆偏振纯度则用圆偏振度描述。我在后处理里定义如下Stokes参数:

  • S0 = |Ex|^2 + |Ey|^2
  • S1 = |Ex|^2 - |Ey|^2
  • S2 = 2 * Re(Ex * conj(Ey))
  • S3 = -2 * Im(Ex * conj(Ey))

归一化到S0后,(S1/S0, S2/S0, S3/S0)就是庞加莱球上的坐标。S3/S0接近+1或-1表示高纯度圆偏振,S1/S0接近±1表示高纯度线偏振。我在双参数扫描中发现,通过调整θ和d,S3可以从-0.85变化到+0.85,已经能覆盖很宽的圆偏振范围;在一些特定的参数点,S3超过0.95,说明远场辐射几乎就是纯圆偏振。

偏振性能还要结合Q因子一起看。一个模式如果Q太低,频谱线宽过大,偏振参数在共振峰附近的频率色散会很剧烈;Q太高则信号弱,探测困难。所以实际器件设计要同时盯住两个指标:目标频率处的偏振纯度和Q因子是否落在合理区间。我个人更倾向于选择一个Q在10^4到10^5之间的工作点,此时远场信号足够强,偏振纯度也可以保持在90%以上。

4.3 鲁棒性分析:几何扰动下的表现

仿真结束前一定要做容差分析,否则实验复现时大概率翻车。我关注三个扰动源:孔半径误差、层间距误差和扭转角加工偏差。对孔半径r做±2%扫描,发现Q因子变化幅度可达20%以上,但远场偏振态的变化相对温和,S3坐标偏移通常在0.1以内。层间距h_gap因为本身是垂直方向参数,对偏振的影响更小,但对Q的影响同样显著。扭转角如果偏离设定值0.5°,偏振态在庞加莱球上的位置会沿原有轨迹移动一小段距离,不会跳跃到无关区域。

这个鲁棒性结果其实提供了非常实用的设计建议:不要把一个器件的工作点定在参数空间的极端尖点处,比如S3恰好等于0.999的位置。那里虽然看起来偏振纯度极高,但微小加工误差会把它拉回0.95甚至更低。选择偏球面上一些平滑区域的参数点,工艺容差会大得多。

5. 我在这个项目里踩过的几个坑,以及对应的调参策略

5.1 周期失配导致的假模式

我第一次尝试双层扭转模型时,为了省事,把两层的周期性条件同时设在同一个正方形晶胞上,第二层只是几何上转了5°,结果算出一个Q高达10^7的"BIC"。当时还挺高兴,后来检查场图才发现,这个模式只在超胞的一小块区域里局域,四周电场不连续,明显是周期边界和几何结构冲突造成的伪模式。

排查方法很简单:把计算得到的电场数据按周期条件平移复制,拼成2×2或3×3的超胞看连续性。如果电场在晶胞边界处出现台阶或跳变,说明你的周期性假设不成立。如果结构内的电场在拼起来后连续流畅,这个模式才可能是真BIC。这个检查花费不多,但能省下后面所有无效的参数扫描时间。

5.2 特征频率虚部不收敛

BIC模式的特点就是虚部极小,这对特征频率求解器的收敛性提出了很高的要求。我遇到的问题是:在粗网格下,虚部表现为一个较小但稳定的值,细化网格后虚部不是单调减小,而是先减小再反弹,导致Q因子始终在10^5附近徘徊,加细网格也突破不了。

后来定位到两个原因。一是PML厚度不足,反射的数值波淹没了真实的辐射损耗;二是网格在两层板之间的空气间隙里太粗,导致上下层之间的电磁耦合计算不准确。解决办法是:先固定PML厚度做网格收敛测试,PML厚度设为2微米以上;空气间隙区域单独设置最大单元尺寸不超过100nm。调整之后,虚部随网格加密呈现正常收敛趋势,Q因子上限才真正释放出来。

5.3 远场偏振相位参考点问题

偏振分析中最隐蔽的坑,是本征模的绝对相位是任意的。同一个模式在两次求解中可能整体相位相差90度,如果不固定参考相位,算出来的S2和S3会来回跳,误以为偏振不稳定。

我的处理方式是不要直接看Ex和Ey的绝对相位,而是关注两个分量之间的相对相位差,以及它们的幅值比。在后处理表达式里写成 phase(Ex)-phase(Ey) 和 abs(Ex)/abs(Ey) 这两个相对量,就不会受整体相位漂移影响。还有一种更稳妥的做法是选定结构内的某个固定点作为参考,输出该点的相位,再把其他所有场量都减去这个参考相位。这样即使不同参数点之间模式相位不连续,偏振判读仍然一致。

最后再分享一个我自己的体会:这类扭转光子晶体的仿真,最忌讳一上来就追求"全参数自动化扫描"。先把θ=0的双层结构调到能复现文献里的BIC,确认网格和边界条件都干净了,再逐步引入扭转角和第二个自由度。每一步都做一次场图和Q值检查,后面出问题时就很容易定位是物理模型的问题,还是纯数值上的问题。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询