1. 仿真目标与整体思路一套模型串起线性手性和非线性增强1.1 这套仿真到底在算什么先给还没完全进入状态的同学交个底标题里这些东西其实是在讨论同一个物理问题——一个介质超表面在什么几何条件下能同时具备两个性质一个是线性光学层面的“本征手性BIC”另一个是非线性光学层面的“三次谐波增强”。如果你最近在调研手性超表面的文献大概率会看到“本征手性连续谱束缚态”这个说法。这里的“本征”两个字是指结构本身在二维平面上就是手性的不需要斜入射不需要外加磁场更不需要把结构做成三维螺旋只要圆偏振光垂直打上去透射谱里就会自然出现左旋圆偏光和右旋圆偏光的差异也就是CD信号。而BIC呢通俗讲是一种被“锁”在结构里的模式理论上它和外部辐射通道完全解耦所以没有辐射损耗Q因子可以趋近无穷但由于它不往外漏光你直接测远场是看不见它的。这时候“手性”的作用就体现出来了——它通过破坏结构的面内镜像对称性把部分BIC从“完全暗态”变成“准BIC”让能量既能被高度局域增强又能在远场被激发和读取。这一个逻辑链条直接决定了整套仿真的思路。我见过不少刚入门的同学一上来就把全部目标堆到一起想一次性得到所有图结果参数乱成一锅粥。我做这套模型的经验是先把它拆成线性本征问题、线性透射问题、非线性谐波问题三个层次每个层次对应特定的计算方法和后处理手段。本文后面所有章节基本就是沿着这条主线展开的。1.2 为什么把三次谐波和手性BIC放在同一个体系里三次谐波THG是非线性光学里最经典的“探针”效应之一。它的物理来源是三阶非线性极化辐射强度理论上正比于局域基频电场强度的六次方——注意是六次方。这意味着如果你能把结构里的局域场增强10倍THG信号理论上能增强10^6倍量级这是一个非常恐怖的非线性增强杠杆。普通超表面的局域场增强靠的是等离激元或者米氏共振增强倍数往往只有几倍到几十倍。但BIC不一样特别是近BIC模式下理论Q因子可达10^5甚至更高对应的局域场增强可以达到两个数量级以上。所以一个自然的想法是把“高Q的准BIC模式”当作非线性增强的平台而“本征手性”则给这个平台增加了偏振选择性——也就是说这种增强只对某一旋向的圆偏光有效。从应用角度看这个组合能解决一个实际问题手性非线性器件通常要求结构同时具备“强手性响应”和“强非线性转换效率”。单独做手性超表面共振峰处的CD信号能做到很高但普通共振Q因子低非线性效率上不去单独做高Q非线性超表面结构往往是对称的对手性没有响应。本征手性BIC走的是另一条路——在手性单元结构里构造BIC等于是在同一套几何里同时调用对称性破坏和高Q束缚态两套物理机制。对仿真来说这意味着你必须同时获得五类数据二维能带图用来确认BIC的存在与位置Q因子图用来量化模式束缚强度手性透射曲线用来验证线性的偏振选择性电场的近场分布图用来定位增强热点和分析模式对称性远场偏振图则用来确认辐射光的偏振状态并指导后续实验测量。1.3 一个容易理解的结构设计出发点既有“手性旋钮”又有“泄漏阀门”在我参考过的多篇介质超表面文献里构造本征手性BIC最常用的手段是取一个原本具有镜像对称性的高Q超表面单元然后在几何上引入一个小扰动比如把方柱的角削掉一块、在圆盘上挖一个斜向椭圆孔、把矩形柱改为Z形或S形等。这个“削角”或者“开孔”的动作干了两件事它破坏了面内镜像对称使模式对圆偏振产生手性选择性同时让原本和外部连续谱完全解耦的BIC产生了一个微弱的辐射通道把Q因子从无穷拉回到有限但依然很高的值。这个逻辑可以用一个很直白的类比来说明把BIC想象成一栋只有一扇门的密室BIC状态就是门反锁、里面的人死活出不来。你如果只想看到这个人永远看不到因为他不透气。现在你在墙上开了一条很小的缝这门缝就是“手性扰动”——人终于能往外走一点Q因子仍然很高因为他大部分时间还是被关着的但你已经有办法感知到他的存在和活动状态了。门缝开得越大信号越强但Q下降得也越快共振增强打折扣。所以几乎所有相关研究都会做一个参数扫描比如从0度到30度去扫这个结构几何参数观察Q因子、CD信号和THG增强之间的权衡。因此仿真时我给读者的第一个建议就是不要急着复现某个具体论文的全部结果先把你的基础结构建好确保能在一或两个参数变化下系统性地从BIC过渡到准BIC再做后续优化。2. 几何建模、材料选型与COMSOL参数化设置2.1 晶胞参数和材料参数怎么起步如果你手里没有现成的器件尺寸从一个相对通用的近红外介质超表面配置起步是比较稳妥的。我用过比较多的是非晶硅或者晶体硅纳米柱波长范围放在1300到1700纳米。这个区间最方便的地方在于硅的折射率约3.4到3.5损耗比较低单元的尺寸大概在亚波长量级网格量和计算时间都比较可控。举个例子假设周期P900纳米柱高H600纳米基础的方形柱或者矩形柱边长约450纳米当你需要它进入本征手性区时最直接的做法是把方形柱的一组对角切掉或者在一个圆盘上引入一个“S形”的镂空。这些几何操作都可以在COMSOL的几何序列里用布尔运算一次性完成。材料设置上记得硅的折射率要按波长变化来设定。如果直接填固定折射率在宽谱透射扫描中会丢掉很多色散细节尤其是BIC共振本身带有很强的模式色散折射率不对会导致共振位置和Q因子都出现偏移。折射率数据可以从常见的光学常数库中提取再在COMSOL材料节点里用插值函数引入。2.2 在COMSOL里建立手性单元并定义扫参进入COMSOL波动光学模块之前先想清楚你要用的是二维还是三维模型。说句实话二维能带图、手性透射曲线和Q因子图很多情况下可以在合理的二维化模型里快速出结果但如果你的结构是介质柱阵列电场在高度方向上有显著分布那我建议直接上三维因为二维模型会丢掉很多垂直方向的物理细节。我实测下来三维模型配合周期性边界条件在网格数量控制得当的情况下计算负担并不像想象中那么可怕。三维建模的关键步骤是几何操作和参数化建立一个长方体单元长宽均为P高度为H。在中心添加一个需要做成手性形状的介质柱体。最简单的做法是建立一个长方体然后在两个对角处添加两个小圆柱或三角块用差集把它们“切掉”。把用于切角的块的尺寸或者旋转角设为全局参数比如d_corner。这个参数最好独立于柱子的整体尺寸方便后面扫参时只改变“不对称度”。在单元上下面添加空气层厚度至少为工作波长的1.5倍以上给倏逝场留有足够衰减空间同时设置周期性端口。模型的四周Z方向用周期性边界条件顶面和底面用端口边界条件。第一次跑的时候可以把周期P和柱高H设为扫描参数但我经验是优先扫“不对称度”这一个参数。因为它才是控制扰动强度、决定Q因子和手性大小的核心旋钮。等你能观察到随扰动参数增大、Q因子下降且CD增强的趋势后再去优化周期和柱高放大效果效率会高很多。2.3 网格和求解器设置上的几个关键点BIC模式之所以难算是因为它的辐射损耗非常小透射谱上的共振峰极窄需要很高的频率分辨率才能解析。COMSOL的频域扫描如果只是均匀步长扫描经常会漏掉共振点更常见的问题是网格太粗时有限元数值泄露会人为增加辐射损耗让BIC模式根本看不出高Q特征。我的建议是分两次求解。第一次用稍微均匀的网格和中等频率步长比如2纳米间隔扫出透射谱上的大概共振位置确定共振位置后再在共振波长附近做局部细化扫描步长缩到0.1到0.01纳米量级同时把网格在介质柱内部和上下表面的边界层加密。注意在BIC模式里电场能量主要集中在介质内部和表面附近空气区网格密度过高意义不大反而会拖慢速度。这里有一个特别容易被忽略的点——边界条件的取向。手性结构对左右旋圆偏振光的响应不同因此你的端口中必须能独立区分LCP和RCP。在COMSOL中使用周期性端口时可以设置两个正交线偏振的本征模式然后通过设置相位关系来合成圆偏振入射。如果你不希望手工操作端口模式也可以直接用背景场替代——在散射场公式下把入射平面波设为Ex分量相位0度、Ey分量相位正负90度就能得到两种不同旋向。第一次算的时候最好把两种入射的几何都跑一遍确保CD信号的符号和手性结构的方向自洽。3. 二维能带、Q因子提取与BIC模式判定3.1 能带图不要只扫Γ点很多人在验证BIC的时候习惯性地只在布里渊区中心点即Γ点附近做频率扫描看到某个模式Q因子很高就下结论说这是BIC。这其实不够严谨尤其当结构引入手性扰动之后原本在Γ点的对称保护BIC可能会移动到布里渊区内部成为驻波状态而不是仍然严格在Γ点。更稳妥的做法是做一个“二维能带图”也就是沿布里渊区边界走一圈路径通常对正方晶格取Γ-X-M-Γ甚至完整遍历第一布里渊区内的波矢。在这个路径上求解特征频率就能看出模式是如何随波矢变化的。一个真正的BIC在能带图上的表现是在特定的波矢点附近出现一条和其他能带没有交叉的平坦带或者高Q模式带且该点在频率上不与任何辐射连续谱重合。在COMSOL里做能带图方法是周期性边界条件里定义波矢分量kx和ky用参数化扫描改变它们对每个波矢做一次特征频率求解然后把得到的本征频率和波矢位置画出来。生成能带图的脚本思路其实很简单难的是确保每条带上的模式都是物理真实的模式。我的排查方式是在每个波矢点保存一个电场分布快照检查模式是否为周期性边界下的布洛赫模式而非局部边界伪影。3.2 Q因子怎么从仿真里提取BIC本身Q是无穷大的但你实际能计算到的只是准BIC。Q因子的提取有两种常用途径。第一种透射谱拟合。在频域扫描得到透射曲线后把共振峰用Fano线型或者Lorentz线型做拟合从线宽估Q λ₀ / Δλ。这个方法直观但对网格要求高窄峰下如果频率扫描步长不够小很容易低估Q。第二种本征值法也是我更推荐的做法。在COMSOL的特征频率求解器中设置包含复数本征值也就是考虑阻尼情况每得到一个模式的复数频率ω ω_r iω_i就可以得到Q ω_r / (2|ω_i|)。这个数值在有限的网格精度下会有一个人为上限一般是10^5到10^7量级。如果你把扰动参数调小Q计算结果不再上升甚至出现振荡那就说明网格已经到极限了需要加密或者换用更高阶的形函数。我一开始做这类仿真时习惯直接看透射谱的峰宽后来写论文才发现审稿人更认本征值提取的Q曲线。两种方法同步使用很有必要因为它们能从不同角度校验模式是否真实。3.3 怎么判断你找的模式真的是手性BIC判断方法可以从两个维度去做。第一个维度是场分布真正的BIC模式通常具有独特的近场对称性其远场辐射分量为零你在近场电场图中会看到一个高局域的涡旋状或节点状结构。当加入手性扰动后准BIC的电场图中会出现不对称化趋势某些热点向切角方向偏移。第二个维度是Q因子随扰动参数的幂律变化理论上一维扰动导致的准BIC其Q和扰动参数的平方成反比即Q ∝ 1/δ²。你可以把不同切角尺寸对应的Q画在双对数坐标里看斜率是否有趋近2的线性区间。有了这两个判断你才能比较有把握地说“我看到的这个高Q模式确实来源于BIC泄漏后的准BIC而不是普通的米氏共振。”4. 手性透射曲线与远场偏振图4.1 圆偏振入射设置的前后一致性问题手性透射曲线本身不难算但初次接触的人最大的困惑是圆偏振的定义。我在仿真代码里习惯把迎着光传播方向看的电场旋转方向作为判据。假如光沿Z方向传播LCP的定义是电场的x分量超前y分量90度数学上写作E E0(x̂ iŷ)。RCP则对应E E0(x̂ - iŷ)。但是要特别提到的是不同文献的坐标定义可能反过来因此你仿真得到的CD正负号并不是物理的第一要素关键是你在透射曲线、远场偏振和后续THG曲线里要保持同一套定义。否则你可能会看到一套矛盾的结果——比如透射谱显示右手圆偏透过率高而远场偏振图却显示共振辐射是左旋圆偏这种冲突九成不是物理问题而是你的坐标定义在某一环节翻转了。实际计算时一次性在模型中加两种旋向的入射背景场可能会占用更多内存。我通常分两次扫描一次LCP、一次RCP用脚本把S参数导出后在绘图阶段再合并成手性透射图。4.2 透射T_LCP和T_RCP能读出什么在手性透射曲线图里完整的物理信息不仅包括同偏振透射T_LL和T_RR还包括交叉偏振透射T_LR和T_RL。对纯本征手性超表面且不考虑非对称反射的情况下交叉偏振项和同偏振项的差异可以用于计算椭偏参数和圆二色性。仿真中如果只设置了单一端口模式一些软件会自动忽略交叉偏振分量这是不够的。我都是用两个正交模式和它们的相位叠加来定义端口矢量这样S参数矩阵可以包含全部四个透射分量。一个典型的本征手性BIC透射谱长什么样在共振波长附近你会看到某一旋向的透射谱中出现一个非常窄的类Fano凹陷或者尖峰而另一旋向的透射谱几乎不变。CD曲线在共振处出现一个符号明确、峰值显著的信号。这里有个容易被误导的地方很多新手会把CD峰位直接等同于BIC共振峰位其实由于Fano干涉背景的存在CD峰的极值位置通常会和透射凹陷位置有微小偏移。不要慌这是离散态和连续背景干涉的正常结果不是计算错误。4.3 远场偏振图到底怎么画从仿真的角度看远场偏振图通常有两种画法。一种是直接在COMSOL的远场节点里计算某个频率下超表面整个周期阵列在远场球面上的电场角分布。另一种是手动推导取近场平面的数据做平面波展开然后对每个衍射级次或者每个k方向的辐射分量提取偏振态。如果你使用COMSOL远场节点最需要注意的问题是远场计算需要一个封闭的积分表面而且电流源J和磁流源M的提取需要保证在结构外部没有其他源。对周期结构来说我习惯在一个完整周期单元周围建立一个立方体积分表面把实际的切向场都提取出来再在该表面上计算远场贡献。提取偏振图的具体操作可以把远场球面上的电场写成球坐标分量E_theta和E_phi然后投影到传播方向的垂直平面分解成左旋和右旋分量的振幅。可以用归一化的Stokes参数S1/S0和S2/S0来定量描述辐射场的偏振态。在准BIC共振频率处你会观察到辐射场以某一旋向为主偏振椭圆度接近1而偏离共振时背景辐射主要是线偏振成分椭圆度降下来。这个现象非常直观也算是本征手性BIC在远场留下的“指纹”。5. 三次谐波增强的原理与建模仿真流程5.1 电场图为什么是三步走里的第一步计算THG之前我强烈建议你先把基频共振处的电场图认认真真看几遍。THG信号正比于局域电场的六次方这个关系决定了并不是结构所有区域都对远场THG有贡献只有那些场强极高的“热点”才是真正贡献非线性信号的区域。画电场图时有几个实用建议对三维模型至少画两个截面一个是z H/2处的水平截面观察面内模式分布和手性扰动导致的电场不对称另一个是沿着x方向的垂直截面观察模式在垂直方向上的局域化程度。颜色标尺用对数坐标因为BIC模式下的场增强可能达到几十到几百倍线性坐标下低场区域会全部变成蓝色看不出热点位置。BIC的电场图有明显特点模式能量高度局域电场强度在柱体内部远大于外部空气区。当你加入手性扰动后场分布中会出现明显的旋转特征。我常见的现象是在某些相角下电场的旋向会切角附近形成“涡旋结构”。这些涡旋往往就是非线性增强的源点。如果这种区域恰好和材料的高非线性系数区域重叠THG效率就会非常好。5.2 三次谐波仿真的“两趟计算”严格全波模拟THG可以理解为一次三阶非线性极化过程基频电场E(ω)在介质中产生非线性极化P_NL(3ω) ε₀χ(3)E(ω)³这个极化在三倍频处充当一个等效电流源向外辐射三倍频光。工程上我推荐把计算分成两步走。第一步是纯线性的基频场求解在COMSOL的频域求解器中将一个周期性超表面结构放在端口激励下频点在基频波长获得整个结构的电场分布Eω(r)。第二步是把非线性极化当作一个已知的“等效源”在三倍频频率处再做一次线性电磁场求解——此时介质的折射率要用三倍频波长处的实际折射率并且把源项设置为从第一步电场推导出的P_NL。直接在COMSOL中一次性同时求解基频和三倍频是一套相当复杂的多物理场耦合多数时候收敛困难和内存开销都会让你怀疑人生。分两步做每一步物理目标单一排查问题也容易。如果你只是想要一个相对THG效率估算而且是非磁性的非手性材料体系也可以走捷径把每个像素的局部增强因子|E_ω(r)/E₀|取六次方再在体积内积分用这个量来表征THG的远场变化趋势。很多文献里的THG增强倍数图其实都源于这个简化模型在缺少高精度实验对比时完全够用。5.3 手性THG曲线怎么从两次仿真中提取当你做完了基频的线性仿真后会很自然地拿到两个量LCP入射下的电场分布E_L(r)和RCP入射下的电场分布E_R(r)。用这两组电场去计算THG源的强度对比可以得到一个估算手性非线性响应强度差。如果需要得到“手性THG转换效率曲线”这种更完整的结果就需要在上述两步法的基础上做一个波长扫描在基频波长附近选一组合适的波长点比如从共振中心向两边各取3到5个点每个波长都做一次线性求解和一次三倍频源辐射求解最后把三倍频频率处的透射功率画成一条曲线。这组曲线做两个旋向各跑一遍就能拿到漂亮的LCP-THG和RCP-THG对比图。这里提醒一个取点问题因为Q因子很高时共振极窄你选的波长扫描点在远离共振时THG效率会快速低到数值噪声水平这会给结果带来很大的误差。此时需要注意增加频率点密度或者在计算THG源时使用相对场增强而不是直接用绝对功率值。另外一个坑是三次谐波频率处材料的折射率。很多材料在近红外基频处损耗很低但在三倍频处可能落入吸收带此时必须使用带虚部的复折射率否则你会得到偏高的THG功率估算。实际操作中我还见过有人直接在基频求解时不加三阶项然后在后处理里用“E^3”的公式去画非线性源分布这本身没问题但如果你的结构里电场增强幅度太大局部功率密度可能已经超过了非微扰阈值那理论预估本身就值得怀疑了。6. 常见问题与排查技巧6.1 这些问题我基本都踩过下面我挑几个反复出现的坑位总结一下按自检优先级整理成表方便你遇到问题时按图索骥。现象主要原因排查/解决方案透射谱共振特别宽Q因子达不到文献量级网格不够细端口边界距离结构太近加密柱体内部网格增加上下空气层厚度改用本征值法提取QLCP和RCP透射曲线完全重合几何结构没有破缺镜像对称或圆偏振定义配置错误检查手性扰动参数是否生效核对端口相位关系本征值求解死活找不到BIC模式扫频范围太窄或模式Q太高迭代不收敛扩大搜索频率范围把扰动参数调到较大值先找到趋势能带图某些带上出现莫名“断带”波矢扫描路径定义有误或者模式追踪算法跳带按特定顺序扫描布里渊区路径借助电场分布辅助判断同一模式的延续CD峰位置和透射凹陷位置不一致Fano线型干涉导致的正常偏移不需要修正报道时说明CD峰位和共振中心差别的来源即可三次谐波功率在远离共振时出现负值或噪声THG信号低于数值噪声本底改用相对增强倍率表示或增加网格加密和迭代精度6.2 三个提升仿真效率的小技巧第一手性BIC结构本身并不需要非常复杂的网格。共振模式高度局域在介质内部这时你要做的是在柱体内部使用比空气层密得多的“局部细化”。可以创建一个只覆盖介质柱的域然后在网格序列中单独给它设置最大单元尺寸把它设为波长的十分之一以内。空气区的网格则可以设为波长的四分之一甚至更粗。这样在不牺牲精度的前提下计算时间能压缩到原来的三分之一以下。第二二维能带图和Q因子扫描的计算量通常不大但每个起点都手动设置扫描会很痛苦。写一个简单的循环脚本把k路径、扰动参数和特征频率搜索范围都参数化。脚本输出后你就可以一次跑出所有需要的数据睡一觉起来直接收结果。第三如果电场的后处理结果异常、热点位置在每次网格调整后变化明显这往往意味着你的结构处于“临界状态”不是BIC发生了偏移就是某个参数刚好落在模式简并点上。这时不要单纯追求加密网格而是去扫描一个控制参数比如增大一点切角尺寸看模式是否能稳定。模式不稳定通常意味着你需要调整模型本身而不是网格。6.3 怎样保证结果的物理可解释性仿真做得再漂亮最后也需要能和实验或者文献对比。在这类超表面计算里实际器件总有制造误差、衬底折射率的影响、入射光的发散角效应这些都会让测量到的Q因子远低于理论值。所以有条件的话在仿真后处理阶段加一个参数化分析稍微改变一下周期大小或者柱高尺寸记录Q因子的变化范围这样能帮助理解实验上的容差范围。如果你准备把这套方法复用在其他波段或者其他材料体系比如把近红外硅换成可见光区的氮化硅或者中红外的砷化镓结构尺寸需要按比例缩放。核心逻辑是保持单元尺寸和波长比基本不变同时注意新材料的折射率实部和虚部即损耗特性是否满足高Q条件。材料折射率虚部只要达到10的负三次方量级就已经可能把Q因子压在几千以内这会直接改变你的“非线性增强杠杆”效果。 COMSOL仿真手性BIC超表面：从线性模式到三次谐波增强
1. 仿真目标与整体思路一套模型串起线性手性和非线性增强1.1 这套仿真到底在算什么先给还没完全进入状态的同学交个底标题里这些东西其实是在讨论同一个物理问题——一个介质超表面在什么几何条件下能同时具备两个性质一个是线性光学层面的“本征手性BIC”另一个是非线性光学层面的“三次谐波增强”。如果你最近在调研手性超表面的文献大概率会看到“本征手性连续谱束缚态”这个说法。这里的“本征”两个字是指结构本身在二维平面上就是手性的不需要斜入射不需要外加磁场更不需要把结构做成三维螺旋只要圆偏振光垂直打上去透射谱里就会自然出现左旋圆偏光和右旋圆偏光的差异也就是CD信号。而BIC呢通俗讲是一种被“锁”在结构里的模式理论上它和外部辐射通道完全解耦所以没有辐射损耗Q因子可以趋近无穷但由于它不往外漏光你直接测远场是看不见它的。这时候“手性”的作用就体现出来了——它通过破坏结构的面内镜像对称性把部分BIC从“完全暗态”变成“准BIC”让能量既能被高度局域增强又能在远场被激发和读取。这一个逻辑链条直接决定了整套仿真的思路。我见过不少刚入门的同学一上来就把全部目标堆到一起想一次性得到所有图结果参数乱成一锅粥。我做这套模型的经验是先把它拆成线性本征问题、线性透射问题、非线性谐波问题三个层次每个层次对应特定的计算方法和后处理手段。本文后面所有章节基本就是沿着这条主线展开的。1.2 为什么把三次谐波和手性BIC放在同一个体系里三次谐波THG是非线性光学里最经典的“探针”效应之一。它的物理来源是三阶非线性极化辐射强度理论上正比于局域基频电场强度的六次方——注意是六次方。这意味着如果你能把结构里的局域场增强10倍THG信号理论上能增强10^6倍量级这是一个非常恐怖的非线性增强杠杆。普通超表面的局域场增强靠的是等离激元或者米氏共振增强倍数往往只有几倍到几十倍。但BIC不一样特别是近BIC模式下理论Q因子可达10^5甚至更高对应的局域场增强可以达到两个数量级以上。所以一个自然的想法是把“高Q的准BIC模式”当作非线性增强的平台而“本征手性”则给这个平台增加了偏振选择性——也就是说这种增强只对某一旋向的圆偏光有效。从应用角度看这个组合能解决一个实际问题手性非线性器件通常要求结构同时具备“强手性响应”和“强非线性转换效率”。单独做手性超表面共振峰处的CD信号能做到很高但普通共振Q因子低非线性效率上不去单独做高Q非线性超表面结构往往是对称的对手性没有响应。本征手性BIC走的是另一条路——在手性单元结构里构造BIC等于是在同一套几何里同时调用对称性破坏和高Q束缚态两套物理机制。对仿真来说这意味着你必须同时获得五类数据二维能带图用来确认BIC的存在与位置Q因子图用来量化模式束缚强度手性透射曲线用来验证线性的偏振选择性电场的近场分布图用来定位增强热点和分析模式对称性远场偏振图则用来确认辐射光的偏振状态并指导后续实验测量。1.3 一个容易理解的结构设计出发点既有“手性旋钮”又有“泄漏阀门”在我参考过的多篇介质超表面文献里构造本征手性BIC最常用的手段是取一个原本具有镜像对称性的高Q超表面单元然后在几何上引入一个小扰动比如把方柱的角削掉一块、在圆盘上挖一个斜向椭圆孔、把矩形柱改为Z形或S形等。这个“削角”或者“开孔”的动作干了两件事它破坏了面内镜像对称使模式对圆偏振产生手性选择性同时让原本和外部连续谱完全解耦的BIC产生了一个微弱的辐射通道把Q因子从无穷拉回到有限但依然很高的值。这个逻辑可以用一个很直白的类比来说明把BIC想象成一栋只有一扇门的密室BIC状态就是门反锁、里面的人死活出不来。你如果只想看到这个人永远看不到因为他不透气。现在你在墙上开了一条很小的缝这门缝就是“手性扰动”——人终于能往外走一点Q因子仍然很高因为他大部分时间还是被关着的但你已经有办法感知到他的存在和活动状态了。门缝开得越大信号越强但Q下降得也越快共振增强打折扣。所以几乎所有相关研究都会做一个参数扫描比如从0度到30度去扫这个结构几何参数观察Q因子、CD信号和THG增强之间的权衡。因此仿真时我给读者的第一个建议就是不要急着复现某个具体论文的全部结果先把你的基础结构建好确保能在一或两个参数变化下系统性地从BIC过渡到准BIC再做后续优化。2. 几何建模、材料选型与COMSOL参数化设置2.1 晶胞参数和材料参数怎么起步如果你手里没有现成的器件尺寸从一个相对通用的近红外介质超表面配置起步是比较稳妥的。我用过比较多的是非晶硅或者晶体硅纳米柱波长范围放在1300到1700纳米。这个区间最方便的地方在于硅的折射率约3.4到3.5损耗比较低单元的尺寸大概在亚波长量级网格量和计算时间都比较可控。举个例子假设周期P900纳米柱高H600纳米基础的方形柱或者矩形柱边长约450纳米当你需要它进入本征手性区时最直接的做法是把方形柱的一组对角切掉或者在一个圆盘上引入一个“S形”的镂空。这些几何操作都可以在COMSOL的几何序列里用布尔运算一次性完成。材料设置上记得硅的折射率要按波长变化来设定。如果直接填固定折射率在宽谱透射扫描中会丢掉很多色散细节尤其是BIC共振本身带有很强的模式色散折射率不对会导致共振位置和Q因子都出现偏移。折射率数据可以从常见的光学常数库中提取再在COMSOL材料节点里用插值函数引入。2.2 在COMSOL里建立手性单元并定义扫参进入COMSOL波动光学模块之前先想清楚你要用的是二维还是三维模型。说句实话二维能带图、手性透射曲线和Q因子图很多情况下可以在合理的二维化模型里快速出结果但如果你的结构是介质柱阵列电场在高度方向上有显著分布那我建议直接上三维因为二维模型会丢掉很多垂直方向的物理细节。我实测下来三维模型配合周期性边界条件在网格数量控制得当的情况下计算负担并不像想象中那么可怕。三维建模的关键步骤是几何操作和参数化建立一个长方体单元长宽均为P高度为H。在中心添加一个需要做成手性形状的介质柱体。最简单的做法是建立一个长方体然后在两个对角处添加两个小圆柱或三角块用差集把它们“切掉”。把用于切角的块的尺寸或者旋转角设为全局参数比如d_corner。这个参数最好独立于柱子的整体尺寸方便后面扫参时只改变“不对称度”。在单元上下面添加空气层厚度至少为工作波长的1.5倍以上给倏逝场留有足够衰减空间同时设置周期性端口。模型的四周Z方向用周期性边界条件顶面和底面用端口边界条件。第一次跑的时候可以把周期P和柱高H设为扫描参数但我经验是优先扫“不对称度”这一个参数。因为它才是控制扰动强度、决定Q因子和手性大小的核心旋钮。等你能观察到随扰动参数增大、Q因子下降且CD增强的趋势后再去优化周期和柱高放大效果效率会高很多。2.3 网格和求解器设置上的几个关键点BIC模式之所以难算是因为它的辐射损耗非常小透射谱上的共振峰极窄需要很高的频率分辨率才能解析。COMSOL的频域扫描如果只是均匀步长扫描经常会漏掉共振点更常见的问题是网格太粗时有限元数值泄露会人为增加辐射损耗让BIC模式根本看不出高Q特征。我的建议是分两次求解。第一次用稍微均匀的网格和中等频率步长比如2纳米间隔扫出透射谱上的大概共振位置确定共振位置后再在共振波长附近做局部细化扫描步长缩到0.1到0.01纳米量级同时把网格在介质柱内部和上下表面的边界层加密。注意在BIC模式里电场能量主要集中在介质内部和表面附近空气区网格密度过高意义不大反而会拖慢速度。这里有一个特别容易被忽略的点——边界条件的取向。手性结构对左右旋圆偏振光的响应不同因此你的端口中必须能独立区分LCP和RCP。在COMSOL中使用周期性端口时可以设置两个正交线偏振的本征模式然后通过设置相位关系来合成圆偏振入射。如果你不希望手工操作端口模式也可以直接用背景场替代——在散射场公式下把入射平面波设为Ex分量相位0度、Ey分量相位正负90度就能得到两种不同旋向。第一次算的时候最好把两种入射的几何都跑一遍确保CD信号的符号和手性结构的方向自洽。3. 二维能带、Q因子提取与BIC模式判定3.1 能带图不要只扫Γ点很多人在验证BIC的时候习惯性地只在布里渊区中心点即Γ点附近做频率扫描看到某个模式Q因子很高就下结论说这是BIC。这其实不够严谨尤其当结构引入手性扰动之后原本在Γ点的对称保护BIC可能会移动到布里渊区内部成为驻波状态而不是仍然严格在Γ点。更稳妥的做法是做一个“二维能带图”也就是沿布里渊区边界走一圈路径通常对正方晶格取Γ-X-M-Γ甚至完整遍历第一布里渊区内的波矢。在这个路径上求解特征频率就能看出模式是如何随波矢变化的。一个真正的BIC在能带图上的表现是在特定的波矢点附近出现一条和其他能带没有交叉的平坦带或者高Q模式带且该点在频率上不与任何辐射连续谱重合。在COMSOL里做能带图方法是周期性边界条件里定义波矢分量kx和ky用参数化扫描改变它们对每个波矢做一次特征频率求解然后把得到的本征频率和波矢位置画出来。生成能带图的脚本思路其实很简单难的是确保每条带上的模式都是物理真实的模式。我的排查方式是在每个波矢点保存一个电场分布快照检查模式是否为周期性边界下的布洛赫模式而非局部边界伪影。3.2 Q因子怎么从仿真里提取BIC本身Q是无穷大的但你实际能计算到的只是准BIC。Q因子的提取有两种常用途径。第一种透射谱拟合。在频域扫描得到透射曲线后把共振峰用Fano线型或者Lorentz线型做拟合从线宽估Q λ₀ / Δλ。这个方法直观但对网格要求高窄峰下如果频率扫描步长不够小很容易低估Q。第二种本征值法也是我更推荐的做法。在COMSOL的特征频率求解器中设置包含复数本征值也就是考虑阻尼情况每得到一个模式的复数频率ω ω_r iω_i就可以得到Q ω_r / (2|ω_i|)。这个数值在有限的网格精度下会有一个人为上限一般是10^5到10^7量级。如果你把扰动参数调小Q计算结果不再上升甚至出现振荡那就说明网格已经到极限了需要加密或者换用更高阶的形函数。我一开始做这类仿真时习惯直接看透射谱的峰宽后来写论文才发现审稿人更认本征值提取的Q曲线。两种方法同步使用很有必要因为它们能从不同角度校验模式是否真实。3.3 怎么判断你找的模式真的是手性BIC判断方法可以从两个维度去做。第一个维度是场分布真正的BIC模式通常具有独特的近场对称性其远场辐射分量为零你在近场电场图中会看到一个高局域的涡旋状或节点状结构。当加入手性扰动后准BIC的电场图中会出现不对称化趋势某些热点向切角方向偏移。第二个维度是Q因子随扰动参数的幂律变化理论上一维扰动导致的准BIC其Q和扰动参数的平方成反比即Q ∝ 1/δ²。你可以把不同切角尺寸对应的Q画在双对数坐标里看斜率是否有趋近2的线性区间。有了这两个判断你才能比较有把握地说“我看到的这个高Q模式确实来源于BIC泄漏后的准BIC而不是普通的米氏共振。”4. 手性透射曲线与远场偏振图4.1 圆偏振入射设置的前后一致性问题手性透射曲线本身不难算但初次接触的人最大的困惑是圆偏振的定义。我在仿真代码里习惯把迎着光传播方向看的电场旋转方向作为判据。假如光沿Z方向传播LCP的定义是电场的x分量超前y分量90度数学上写作E E0(x̂ iŷ)。RCP则对应E E0(x̂ - iŷ)。但是要特别提到的是不同文献的坐标定义可能反过来因此你仿真得到的CD正负号并不是物理的第一要素关键是你在透射曲线、远场偏振和后续THG曲线里要保持同一套定义。否则你可能会看到一套矛盾的结果——比如透射谱显示右手圆偏透过率高而远场偏振图却显示共振辐射是左旋圆偏这种冲突九成不是物理问题而是你的坐标定义在某一环节翻转了。实际计算时一次性在模型中加两种旋向的入射背景场可能会占用更多内存。我通常分两次扫描一次LCP、一次RCP用脚本把S参数导出后在绘图阶段再合并成手性透射图。4.2 透射T_LCP和T_RCP能读出什么在手性透射曲线图里完整的物理信息不仅包括同偏振透射T_LL和T_RR还包括交叉偏振透射T_LR和T_RL。对纯本征手性超表面且不考虑非对称反射的情况下交叉偏振项和同偏振项的差异可以用于计算椭偏参数和圆二色性。仿真中如果只设置了单一端口模式一些软件会自动忽略交叉偏振分量这是不够的。我都是用两个正交模式和它们的相位叠加来定义端口矢量这样S参数矩阵可以包含全部四个透射分量。一个典型的本征手性BIC透射谱长什么样在共振波长附近你会看到某一旋向的透射谱中出现一个非常窄的类Fano凹陷或者尖峰而另一旋向的透射谱几乎不变。CD曲线在共振处出现一个符号明确、峰值显著的信号。这里有个容易被误导的地方很多新手会把CD峰位直接等同于BIC共振峰位其实由于Fano干涉背景的存在CD峰的极值位置通常会和透射凹陷位置有微小偏移。不要慌这是离散态和连续背景干涉的正常结果不是计算错误。4.3 远场偏振图到底怎么画从仿真的角度看远场偏振图通常有两种画法。一种是直接在COMSOL的远场节点里计算某个频率下超表面整个周期阵列在远场球面上的电场角分布。另一种是手动推导取近场平面的数据做平面波展开然后对每个衍射级次或者每个k方向的辐射分量提取偏振态。如果你使用COMSOL远场节点最需要注意的问题是远场计算需要一个封闭的积分表面而且电流源J和磁流源M的提取需要保证在结构外部没有其他源。对周期结构来说我习惯在一个完整周期单元周围建立一个立方体积分表面把实际的切向场都提取出来再在该表面上计算远场贡献。提取偏振图的具体操作可以把远场球面上的电场写成球坐标分量E_theta和E_phi然后投影到传播方向的垂直平面分解成左旋和右旋分量的振幅。可以用归一化的Stokes参数S1/S0和S2/S0来定量描述辐射场的偏振态。在准BIC共振频率处你会观察到辐射场以某一旋向为主偏振椭圆度接近1而偏离共振时背景辐射主要是线偏振成分椭圆度降下来。这个现象非常直观也算是本征手性BIC在远场留下的“指纹”。5. 三次谐波增强的原理与建模仿真流程5.1 电场图为什么是三步走里的第一步计算THG之前我强烈建议你先把基频共振处的电场图认认真真看几遍。THG信号正比于局域电场的六次方这个关系决定了并不是结构所有区域都对远场THG有贡献只有那些场强极高的“热点”才是真正贡献非线性信号的区域。画电场图时有几个实用建议对三维模型至少画两个截面一个是z H/2处的水平截面观察面内模式分布和手性扰动导致的电场不对称另一个是沿着x方向的垂直截面观察模式在垂直方向上的局域化程度。颜色标尺用对数坐标因为BIC模式下的场增强可能达到几十到几百倍线性坐标下低场区域会全部变成蓝色看不出热点位置。BIC的电场图有明显特点模式能量高度局域电场强度在柱体内部远大于外部空气区。当你加入手性扰动后场分布中会出现明显的旋转特征。我常见的现象是在某些相角下电场的旋向会切角附近形成“涡旋结构”。这些涡旋往往就是非线性增强的源点。如果这种区域恰好和材料的高非线性系数区域重叠THG效率就会非常好。5.2 三次谐波仿真的“两趟计算”严格全波模拟THG可以理解为一次三阶非线性极化过程基频电场E(ω)在介质中产生非线性极化P_NL(3ω) ε₀χ(3)E(ω)³这个极化在三倍频处充当一个等效电流源向外辐射三倍频光。工程上我推荐把计算分成两步走。第一步是纯线性的基频场求解在COMSOL的频域求解器中将一个周期性超表面结构放在端口激励下频点在基频波长获得整个结构的电场分布Eω(r)。第二步是把非线性极化当作一个已知的“等效源”在三倍频频率处再做一次线性电磁场求解——此时介质的折射率要用三倍频波长处的实际折射率并且把源项设置为从第一步电场推导出的P_NL。直接在COMSOL中一次性同时求解基频和三倍频是一套相当复杂的多物理场耦合多数时候收敛困难和内存开销都会让你怀疑人生。分两步做每一步物理目标单一排查问题也容易。如果你只是想要一个相对THG效率估算而且是非磁性的非手性材料体系也可以走捷径把每个像素的局部增强因子|E_ω(r)/E₀|取六次方再在体积内积分用这个量来表征THG的远场变化趋势。很多文献里的THG增强倍数图其实都源于这个简化模型在缺少高精度实验对比时完全够用。5.3 手性THG曲线怎么从两次仿真中提取当你做完了基频的线性仿真后会很自然地拿到两个量LCP入射下的电场分布E_L(r)和RCP入射下的电场分布E_R(r)。用这两组电场去计算THG源的强度对比可以得到一个估算手性非线性响应强度差。如果需要得到“手性THG转换效率曲线”这种更完整的结果就需要在上述两步法的基础上做一个波长扫描在基频波长附近选一组合适的波长点比如从共振中心向两边各取3到5个点每个波长都做一次线性求解和一次三倍频源辐射求解最后把三倍频频率处的透射功率画成一条曲线。这组曲线做两个旋向各跑一遍就能拿到漂亮的LCP-THG和RCP-THG对比图。这里提醒一个取点问题因为Q因子很高时共振极窄你选的波长扫描点在远离共振时THG效率会快速低到数值噪声水平这会给结果带来很大的误差。此时需要注意增加频率点密度或者在计算THG源时使用相对场增强而不是直接用绝对功率值。另外一个坑是三次谐波频率处材料的折射率。很多材料在近红外基频处损耗很低但在三倍频处可能落入吸收带此时必须使用带虚部的复折射率否则你会得到偏高的THG功率估算。实际操作中我还见过有人直接在基频求解时不加三阶项然后在后处理里用“E^3”的公式去画非线性源分布这本身没问题但如果你的结构里电场增强幅度太大局部功率密度可能已经超过了非微扰阈值那理论预估本身就值得怀疑了。6. 常见问题与排查技巧6.1 这些问题我基本都踩过下面我挑几个反复出现的坑位总结一下按自检优先级整理成表方便你遇到问题时按图索骥。现象主要原因排查/解决方案透射谱共振特别宽Q因子达不到文献量级网格不够细端口边界距离结构太近加密柱体内部网格增加上下空气层厚度改用本征值法提取QLCP和RCP透射曲线完全重合几何结构没有破缺镜像对称或圆偏振定义配置错误检查手性扰动参数是否生效核对端口相位关系本征值求解死活找不到BIC模式扫频范围太窄或模式Q太高迭代不收敛扩大搜索频率范围把扰动参数调到较大值先找到趋势能带图某些带上出现莫名“断带”波矢扫描路径定义有误或者模式追踪算法跳带按特定顺序扫描布里渊区路径借助电场分布辅助判断同一模式的延续CD峰位置和透射凹陷位置不一致Fano线型干涉导致的正常偏移不需要修正报道时说明CD峰位和共振中心差别的来源即可三次谐波功率在远离共振时出现负值或噪声THG信号低于数值噪声本底改用相对增强倍率表示或增加网格加密和迭代精度6.2 三个提升仿真效率的小技巧第一手性BIC结构本身并不需要非常复杂的网格。共振模式高度局域在介质内部这时你要做的是在柱体内部使用比空气层密得多的“局部细化”。可以创建一个只覆盖介质柱的域然后在网格序列中单独给它设置最大单元尺寸把它设为波长的十分之一以内。空气区的网格则可以设为波长的四分之一甚至更粗。这样在不牺牲精度的前提下计算时间能压缩到原来的三分之一以下。第二二维能带图和Q因子扫描的计算量通常不大但每个起点都手动设置扫描会很痛苦。写一个简单的循环脚本把k路径、扰动参数和特征频率搜索范围都参数化。脚本输出后你就可以一次跑出所有需要的数据睡一觉起来直接收结果。第三如果电场的后处理结果异常、热点位置在每次网格调整后变化明显这往往意味着你的结构处于“临界状态”不是BIC发生了偏移就是某个参数刚好落在模式简并点上。这时不要单纯追求加密网格而是去扫描一个控制参数比如增大一点切角尺寸看模式是否能稳定。模式不稳定通常意味着你需要调整模型本身而不是网格。6.3 怎样保证结果的物理可解释性仿真做得再漂亮最后也需要能和实验或者文献对比。在这类超表面计算里实际器件总有制造误差、衬底折射率的影响、入射光的发散角效应这些都会让测量到的Q因子远低于理论值。所以有条件的话在仿真后处理阶段加一个参数化分析稍微改变一下周期大小或者柱高尺寸记录Q因子的变化范围这样能帮助理解实验上的容差范围。如果你准备把这套方法复用在其他波段或者其他材料体系比如把近红外硅换成可见光区的氮化硅或者中红外的砷化镓结构尺寸需要按比例缩放。核心逻辑是保持单元尺寸和波长比基本不变同时注意新材料的折射率实部和虚部即损耗特性是否满足高Q条件。材料折射率虚部只要达到10的负三次方量级就已经可能把Q因子压在几千以内这会直接改变你的“非线性增强杠杆”效果。