做光子晶体板仿真的人十有八九都会在COMSOL里跑能带、算Q值。可有段时间我一直在琢磨另一个东西例外点Exceptional PointEP。这东西在非厄米光学、拓扑光子学和灵敏传感里都是热点但真正到手去仿资料却少得可怜。我前阵子用COMSOL的波动光学模块搭了一块正方晶格空气孔硅板在孔里引入增益介质把系统弄成非厄米然后在参数平面上反复扫波矢和损耗系数最后终于在复频率图上看到了两个本征模式“撞”在一起的那个点。这篇东西就是整个折腾过程的复盘EP是什么、为什么值得做以及在COMSOL里怎么建模型、怎么扫参、怎么验证一次讲透。1. 光子晶体板EP是什么先解决“为什么要折腾”1.1 光子晶体板不只有能带还有复频率光子晶体板说白了就是一块厚度在波长量级、面内周期性排布的介质板。和无限周期的二维光子晶体不同板在纵向z方向靠全内反射把光限制住所以它既能像波导一样引导光又保留着面内的能带结构。可是一旦模式频率落在光锥上方它就会耦合到自由空间辐射出去这时候本征频率就不能只当实数看了——COMSOL里算出来的特征频率是复数实部决定振荡快慢虚部对应模式寿命或者辐射损耗。很多人跑光子晶体板只画能带实部把虚部丢在一边这其实浪费了很大一部分信息。特别是做传感和激光器的时候虚部甚至比实部更关键高Q模式意味着光子能在板里待更久传感器的分辨率就高激光器则希望所有模式里只有一个“赢家”。而EP这个奇点恰恰是模式的实部和虚部同时在复平面上“纠缠”出来的点所以想研究它必须习惯把能带图看成复平面上的曲面而不是一条平面曲线。1.2 EP两个模式“合体”的复平面奇点先看一个最简单但最常用的模型。两个模式构成的子空间可以用一个2×2非厄米矩阵来描述H ω0 I δ σ_z i γ σ_x这里的δ可以理解成两个模式原本的频率失谐γ表示增益/损耗造成的耦合或衰减。它的本征值是ω± ω0 ± sqrt(δ² - γ²)当δ γ时开根号是实数两个频率的实部不同、虚部相同当δ γ时根号变成虚数两个模式的实部相同、虚部开始分裂。真正形成“一分为二”和“合二为一”的位置就发生在δ γ的那一刻。这时两条本征值曲线在复平面上交于一点同时两个本征矢也合并成一个矩阵无法对角化——这个点就是例外点简称EP。如果把常见的简并点Dirac点拿来对比Dirac点是两个本征值相等但本征矢仍然正交独立EP是本征值和本征矢一起合并基本相当于“两个不同的人突然变成了一个人”。在COMSOL里Dirac点随便扫周期性结构就可能碰到但EP往往需要引入损耗/增益、或者调整几何参数让系统偏离厄米难度明显上一个台阶。一个比较直观的类比是两列同频的波一个在损耗区里跑一个在增益区里跑调到一个特定的平衡点后它们的“对称性身份”会互换。这就是EP在直观图像上最迷人的地方。1.3 EP在光子晶体板里有什么用EP之所以值得专门去COMSOL里找是因为它带来的物理效应非常锋利。最典型的是高灵敏度传感在普通简并点附近微扰导致的频率分裂正比于扰动强度ε但在EP附近频率分裂正比于ε的平方根。ε很小的时候平方根远远大于线性项所以EP传感器理论上能把灵敏度提升一个量级。光子晶体板本身就是非常好的传感平台把EP嵌进去就等于在一个已经很灵敏的平台上再装一个“放大器”。此外EP在激光器里也有用。光子晶体板激光器经常面临多模式竞争如果能把两个模式刚好调到EP模式之间的耦合会重新分配增益实现单模/单偏振输出。还有拓扑光子学绕EP在参数空间走一圈两个本征矢会交换等价于一个1/2拓扑荷可以用于拓扑模式调控。所以别觉得EP只是理论概念它在COMSOL里被“仿真出来”之后很多应用层面的设计才能真正落地。2. COMSOL里搭光子晶体板几何、材料与边界设置2.1 单元胞建模与参数化先把基础几何立起来我建议即使你最终关心的是二维等效问题也先从三维单胞开始做。三维单胞更接近真实的光子晶体板纵向辐射损耗天然包含在内而EP经常恰恰出现在辐射通道和导模耦合的位置。COMSOL里选“组件 三维”用块Block做介质板用圆柱Cylinder做空气孔再用布尔操作从板里把孔挖掉。所有尺寸用全局参数定义不要写死后面扫参的底气全靠参数化。我用的初始参数表大概是这样参数数值说明晶格常数 a600 nm面内周期板厚 t220 nm板厚度孔半径 r120 nm孔半径可扫描硅折射率3.45 0.0001i微小损耗孔介质折射率1.5 - i*gammagamma为增益/损耗扫描参数空气折射率1.0上下包层及背景孔内材料折射率里的gamma是我最喜欢的扫描参数。把r和gamma都设成全局参数这样在探索EP时可以在“几何面积变化”和“材料虚部变化”两条路线之间迅速切换。COMSOL官方案例库里也有光子晶体能带相关的案例可以拿来做起点但那些案例通常算的是无损耗结构我们需要在此基础上把复折射率加上去。2.2 材料折射率实部与虚部非厄米的来源EP的出现依赖非厄米性。最直接的办法是在结构里引入增益和损耗区。比如正方晶格空气孔硅板我把孔填充材料的折射率定义为n_hole 1.5 - i * gammagamma 0表示增益gamma 0表示损耗。COMSOL里这个虚部的正负号约定要看时谐因子是e^{jωt}还是e^{-iωt}不同版本、不同物理场接口可能有差异。我做仿真时一般先用一个已知案例验证符号确认一下到底哪种设置对应衰减、哪种对应放大然后再放心扫参。更严格的PT对称结构要求折射率实部在空间上偶对称、虚部奇对称。所以如果想做PT对称光子晶体板通常是单元胞里一半孔填充增益材料-iγ_gain另一半孔填充损耗材料iγ_loss让增益和损耗互补。刚开始探索EP时不一定要一步到位构建严格的PT对称结构我一般先扫描整体gamma值在很多带交叉点附近也能观察到两个模式实部或虚部合并的迹象那往往就是EP的前兆。2.3 周期边界与PML如何把“无穷大”塞进一个单胞这是整个建模里最需要耐心的一步。面内方向也就是x和y方向使用COMSOL里的“周期条件”类型选“Floquet”输入两个方向的Bloch波矢kx、ky。波矢不能直接写死我会把它定义成扫描参数。比如研究Gamma-X方向可以令kx kx_scan * 2π/aky 0研究Gamma-M方向就让kx ky kx_scan * 2π/a。z方向由于板上下存在辐射需要在单胞上下各加一层PML。PML不是真实材料本质是一种坐标拉伸算法所以在COMSOL里把上下两个域设成“PML”域指定拉伸方向为正z和负z。PML外侧不需要再做复杂的吸收边界直接设成理想导体PEC就行。关于PML厚度我习惯至少一个工作波长第一次试算时可以从0.5λ一直扫到2λ看目标模式的虚部是否稳定。PML设置不当会在低频段冒出一堆伪模式后患无穷。2.4 网格与求解器不追求一次算到最准网格密度方面光子晶体板介质区域内建议至少保证每波长6个单元孔边界尤其是孔与板的交界要局部加密因为那里的场分布变化非常快。PML区域我常用“扫掠”网格沿z方向分5~10层不需要太多。如果全部用自由四面体PML层的网格会非常密计算量很不划算。求解器方面特征频率研究默认的求解器一般够用。我会把“所需特征数”设置为6~12个搜索基准频率设置在目标带边频率附近搜索范围不要太大否则ARPACK会找出一堆无关的高阶模。开启“模式追踪”之后参数扫描能自动跟随相近的模式这对扫EP非常有帮助。不过模式追踪也不是万能遇到模式高度简并时照样会掉线后面我专门写一节怎么处理。3. 从原理到COMSOL怎样把EP“扫”出来并验证3.1 先把能带图换成复频率图特征频率研究跑完COMSOL默认结果里会给出一个复数值freq。后处理时我建议定义两个变量freq_r real(freq) freq_i imag(freq)然后做一维绘图组横轴设为kx或gamma纵轴同时画freq_r和freq_i。这里有个亲测好用的经验EP附近的模式颜色特别容易互换所以我在同一张图里用不同颜色和线型去区分模式。比如模式1的实部用红实线、虚部用红虚线模式2的实部用蓝实线、虚部用蓝虚线。这样复频率的“接近”和“交换”能一眼看出来而不是在密密麻麻的曲线里到处找。如果你和我一样更喜欢看归一化频率可以再定义一个变量f_norm freq * a / c_const画图时纵轴用f_norm横轴用kx/(2π/a)。这样能带图直接能和论文里的图对比避免整天盯着10的14次方这种大数值看得头晕。3.2 参数扫描策略先粗后细、由简到繁想高效找到EP我总结的套路是“两步走”。第一步在无损无增益gamma 0的情况下沿高对称路径扫描kx寻找两个模式在实部频率上接近甚至交叉的位置。这个交叉点就是EP的候选区因为后续外加非厄米调制能在这里把交叉变成排斥或者在复平面上形成合并。第二步固定kx在这个候选点附近扫描gamma观察两个模式的freq_r和freq_i。如果存在某个gamma值让freq_r合并或者freq_i合并同时另一部分开始分裂那就是PT相变点也就是EP。实际操作中COMSOL的“参数化扫描”可以设置主扫描和辅助扫描。但我个人更推荐先分别做两遍单变量扫描先固定gamma0扫kx再固定kx扫gamma。因为同时扫描两个变量时如果网格不够密模式追踪很容易乱跳输出数据也很难判断对应关系。把两轮数据合并到同一张图里逻辑会非常清楚。3.3 验证EP的三种“实锤”声称找到了EP不能只靠频率曲线“看起来碰在一起”。数值仿真中两个频率差小于0.1%时肉眼根本无法区分到底是不是真的合并。我一般用三种方法交叉验证前两种最常用判据一两个模式的本征矢线性相关。在COMSOL里分别导出两个模式的电场分量归一化后计算内积。如果模长接近1说明本征矢合并如果接近0说明是正交简并。判据二参数环回。在(kx, gamma)参数平面上以EP候选点为中心取一个很小的闭环扫描一整圈把两个模式的复频率分别标成红蓝两点。绕一圈回到原点后如果红点变成了蓝点、蓝点变成了红点那就是EP的拓扑特征。Dirac点绕一圈不会有这种交换。判据三局域等效哈密顿量不可对角化这是更理论的判断。COMSOL不直接给Jordan标准型所以我基本靠前两个判据做工程确认。最容易被忽略的是本征矢判据。很多人看到一个频率交叉就兴奋结果一算内积发现两个模式还是正交的那只是Dirac点不是EP。所以这个验证步骤无论如何不能省。3.4 后处理延伸从复频率到等效电路与导纳-阻抗曲线很多初学COMSOL的朋友会问仿真谐振器的时候为什么大家总喜欢看导纳/阻抗曲线而不是复频率其实两者是相通的。光子晶体板的一个模式可以等效成RLC并联谐振回路谐振频率和品质因子分别是f0 freq_r Q |freq_r| / (2|freq_i|)如果你在COMSOL中做频域分析扫频得到导纳Y(f) G(f) iB(f)那么阻抗Z(f) 1/Y(f)。在“全局计算”里可以直接写Z_re real(1/Y) Z_im imag(1/Y)然后画Z_re和Z_im随频率的曲线阻抗虚部过零、实部出现峰的位置就是谐振频率。这个换算在微波谐振器仿真里非常常用很多朋友搜“如何从导纳曲线换算成阻抗曲线”本质就是这一步。在EP研究里导纳圆图也会出现临界耦合式的单峰收拢用电路语言描述模式行为会更直观尤其是跟射频器件交叉的时候。4. 手把手案例PT对称光子晶体板的EP扫描4.1 模型结构与参数表这一节给一个可以在普通工作站上跑完一轮扫参的演示结构。目标不是复现某篇论文而是把流程走通。结构是正方晶格空气孔硅板孔内填充增益介质上下包层均为空气。参数表如下参数数值说明晶格常数 a600 nm面内周期板厚 t220 nm板厚度孔半径 r120 nm孔半径硅折射率3.45 0.0001i微小损耗孔介质折射率1.5 - i*gammagamma为扫描参数空气折射率1.0包层/背景波矢扫描kx从0到0.5*(2π/a)沿Gamma-X路径gamma扫描0到0.05步长0.002增益强度需要说明的是这个参数组合不一定在你的目标频率范围直接出EP但流程是通用的。如果第一轮扫完没有观察到频率合并可以先改变孔半径比如从100 nm到140 nm之间扫一个粗扫描找到两个模式能量最接近的r值再回来细扫gamma。4.2 COMSOL操作步骤可以直接照着做第一步新建模型组件选三维物理场选“波动光学 电磁波、频域”研究选“特征频率”。在全局定义里把上表的参数全部填好。第二步创建单元胞。用“块”画一个长宽都为a、厚度为t的长方体作为介质板。再用“圆柱”画一个半径r、高度略大于t的柱体用“差集”把柱体从板中减掉形成空气孔。第三步分配材料。板区域使用硅的复折射率孔区域单独建一种材料折射率表达式写成1.5 - i*gamma。注意这里的gamma是一个全局参数材料属性能直接引用。第四步加边界条件。四个侧面x0、xa、y0、ya加“周期条件”类型选Floquet波矢分量设置成kx和ky。在板上方和下方各添加一个PML域PML外侧边界设为PEC。第五步划分网格。板区域用自由四面体最大单元尺寸设为a/6孔边界加一个“大小”节点把最大单元压到a/10左右。PML区域用扫掠网格沿z方向分6层。第六步设置特征频率研究。所需特征数填8搜索基准频率用 c/(a*n_eff) 估计n_eff先取2.5。打开“模式追踪”让扫参时能跟踪相近模式。第七步先固定gamma 0参数化扫描kx记录能带。然后把kx固定到能带交叉点对应的值再参数化扫描gamma。4.3 结果怎么看EP出现时的典型图谱对上面这种PT对称结构扫gamma时你会看到很典型的演变gamma很小时两个模式的freq_r保持一定频率差freq_i都接近零。随着gamma增大到临界值freq_r差逐渐减小在gamma_EP处几乎重合而freq_i开始明显分裂或合并。另一种情况是freq_r不合并、freq_i合并这取决于你站在PT相变的哪一侧。继续增大gammafreq_r重新分开但这一次两个分支的颜色很可能发生了交换——因为在EP处模式的“身份”已经互换了。在一张复频率图上这个“先合后分”的过程非常直观。我强烈建议把freq_r和freq_i放在同一张图的两个y轴上否则只看实部很容易把EP误判成普通带交叉。如果发现gamma_EP对应的频率分裂仍然比较大说明kx离真正的EP还有距离需要回到第一步重新扫kx把候选点修得更准。4.4 进一步验证绕参数环回到原点EP最迷人的地方是它的拓扑性质。为了确认刚才找到的点确实是EP而不是Dirac点我习惯在参数平面上画一个小圆kx kx_EP deltacos(theta) gamma gamma_EP deltasin(theta)这里的delta要选适当大小最好让频率分裂约为本征频率的0.1%。theta从0扫到360度步长5度。每个theta点重新求解特征值然后把所有本征频率画在复平面上。绕一整圈后你会看到模式1的复频率变成了模式2的复频率反之亦然。这就是EP的拓扑“身份交换”。COMSOL不会自动告诉你“模式颜色有没有交换”所以绕圈扫参时的模式追踪要格外小心。我用的土办法是给每个模式编号扫完后单独检查编号是否互换。绕圈半径也不能太大否则可能同时跨过多个EP太小又会淹没在数值噪音里。这个delta的选取基本要靠多试几次。5. 常见问题与避坑经验都是真金白银换来的5.1 PML导致虚部异常怎么办现象扫参时本来应该低损耗的导模虚部突然变得很大而且这些“坏点”在高频段特别集中。原因通常是PML层数不够或者PML域网格太粗导致PML本身成了伪谐振腔。对策先把PML厚度设为2λ网格在z方向至少5层然后逐步增加层数看看目标模式的虚部是否收敛。如果换了PML厚度虚部就跳变那这个模式很可能是PML伪模不是光子晶体板的真实模式。判断方法也很简单把PML厚度翻倍真正模式的频率变化应该很小伪模通常会有明显漂移。5.2 参数扫描时模式追踪“掉线”我用COMSOL跑EP扫参最头疼的就是模式追踪掉线。具体表现是频率曲线在某一步突然从低频跳到高频或者颜色曲线断成几截。原因在于“模式追踪”依赖模式之间的重叠积分而EP附近两个模式高度简并追踪器分不清谁是谁。解决办法有三个方向一是缩小扫描步长把步长控制在预期频率分裂的十分之一二是在扫描前先在一个参考点手动识别模式顺序再打开模式追踪三是把特征值搜索范围收紧只搜目标频率附近±5%。如果还是乱跳就干脆手动分批扫描一组一组导出数据离线用MATLAB或者Excel拼接。虽然笨一点但结果可控。5.3 怎样区分EP与Dirac点前面提到过这里再说一个可落地的实操方法。在COMSOL的“派生值”里定义一个积分算子对单元胞内的电场分量做归一化内积⟨E1|E2⟩ ∫ E1*·E2 dV / (sqrt(∫|E1|² dV) * sqrt(∫|E2|² dV))如果模长接近1两个本征矢重合EP如果接近0正交简并Dirac点。我第一次找到候选点后直接算这个内积结果是0.97于是确认是EP。后来又在一个看似更漂亮的交叉点上算结果是0.003才发现那只是普通简并。纸上谈兵和实际动手的区别往往就差这一步计算。5.4 计算资源与软件使用的几个实在建议三维光子晶体板、复折射率、PML、参数化扫描这几个词叠在一起计算量很快就能吃满一台普通工作站。我的经验是先用粗网格把EP大概位置摸出来再在EP附近局部加密重新扫一轮。这样一轮从两小时能压缩到二十分钟以内。如果只是验证算法先用二维有效折射率模型跑通流程再换三维完整模型也是很实用的变通思路。关于软件安装和授权我一直建议走正规渠道。学校或者单位买好license之后插件和案例库都齐全更新也稳定。不要用网上来路不明所谓“绿色版”“破解版”除了版权风险数值求解器是否被改动过也没法验证。科研数据可信度比省那点安装时间重要得多。5.5 关于单位制别被数量级绕晕COMSOL默认频率单位是Hz光子晶体板这种亚微米结构频率动辄10^14 Hz看久了确实容易眼花。我建议从一开始就定义归一化频率变量。在“定义”里写f_norm freq * a / c_const把所有能带图、EP扫描图都用这个变量做纵轴横轴用kx归一化到2π/a。这样数据和文献完全对齐后面做Q值计算也更方便Q 0.5 * abs(f_norm_r / f_norm_i)公式在COMSOL里可以直接写注意括号层级就行。这个习惯我是在吃过几次数值量级的亏之后才养成的现在强烈推荐给每个做光子晶体仿真的人。第一次跑出EP那个点的时候我对着复频率图盯了好一会儿。这玩意儿和平时看到的能带交叉完全不一样——它不是一个点上的“相遇”而是两条本征值支的黎曼面真正拧在一起绕一圈还互换身份。做仿真最大的甜头就是把抽象的非厄米拓扑变成眼前一张可以反复验证的图。如果你也在研究光子晶体板强烈建议先把这个流程跑通再去做更复杂的结构。后面我打算把EP附近的拓扑荷分布和远场辐射模式的关系再整理一下有结果了再回来更新。