二阶光子拓扑绝缘体最迷人的地方在于几乎每个做这个方向的人第一次在仿真里看到角态场分布图时都会盯着那条局域在顶角处的亮斑看很久。角态不是凭空出现的——它在能谱里是落在带隙中的孤立频点在场分布里是集中在角落的零维束缚模式。这篇博文记录我用Comsol 5.4复现南京大学陈延峰老师组二阶光子拓扑绝缘体工作的完整过程从物理模型搭建、特征频率计算到本征场分布验证把每一步的操作思路、参数由来和坑位都摊开讲清楚。无论你是刚接触拓扑光子学的学生还是已经在用Comsol做电磁仿真但想碰一碰拓扑体系的工程师这篇文章都应该能给你省下不少摸索时间。1. 为什么要算角态二阶拓扑绝缘体的物理逻辑与仿真价值1.1 从体-边对应到体-角对应传统拓扑绝缘体最著名的特征是体-边对应一个拓扑非平庸的二维体系其体能带带有非零的陈数或Z2拓扑不变量必然导致在边界上出现单向传输的边态。而二阶拓扑绝缘体把这个对应关系往前推了一步——它不仅边界上可能有态在边界与边界的交界处也就是角上还会出现零维的角态。用一个不太严格但很好理解的说法角态是边界的边界就像一维SSH链的端态被推广到了二维的角上。我第一次看到这个概念时脑子里冒出的问题很直接这和传统的缺陷态、杂质态有什么本质区别区别在于鲁棒性。角态不是靠某个局部缺陷偶尔产生的束缚态而是受结构对称性和拓扑不变量保护的。说得再直白一点只要你没有破坏特定的离散对称性哪怕把边界切得歪一点、把两个角附近的介质柱挪一挪位置角态依然稳定存在。仿真中验证这一点最有说服力的方式就是算一个有限尺寸阵列的本征模式看带隙里有没有孤立的态再看它的场是不是压在某个角上。1.2 光子体系为什么适合做这件事光子晶体和介质谐振器阵列是实现二阶拓扑绝缘体的天然平台。和电子体系相比光子系统有几个很明显的好处第一麦克斯韦方程组本身是线性的用Comsol直接求本征频率就能得到全部模式不像电子体系要处理复杂的费米面填充问题第二微波段的介质柱阵列尺寸在厘米量级加工容易、测试方便而且场分布可以通过近场扫描直接看到第三各向同性的介质材料在微波段损耗很低仿真和实验能对上。陈延峰老师组的工作在这个方向做了很多关键推进其中一个让我印象深刻的点是他们把二阶拓扑绝缘体的角态从理论模型推到了实验测量在介质谐振器阵列里观测到了角态的本征场分布。所以用Comsol做数值复现的时候目标也很明确先把能谱算对再把场分布画出来最后对照文献里的结果确认角态频率和场形。1.3 用Comsol做这件事的合理路径Comsol在解决这类问题上有它独特的优势几何参数化方便特征频率研究配合参数扫描可以一次性把不同耦合强度下的能谱全算出来本征场直接三维可视化不需要额外导出数据再画图后处理可以做功率流、相位分析对判断角态性质很有帮助。5.4版虽然在网格剖分算法上不是最新但对于这种规模和物理场类型算力需求完全不是问题。下面我从物理模型讲起然后逐步落到Comsol的每一个设置按钮上。2. 二维SSH晶格角态模型的构造逻辑与参数依据2.1 从一维SSH到二维晶胞怎么排做二阶光子拓扑绝缘体目前最主流的可复现方案之一就是二维SSH晶格也叫三聚化格子、呼吸型格子。它的基本思路在一维SSH模型里一条链上的格点以疏-密交替排列形成强弱交替的耦合。把这个想法放到二维正方格子上每个晶胞里放4个介质柱分别靠近正方形的四个角但每个角上的柱不完全落在角点而是向内缩进一定距离。我用表格把关键参数先列出来后面所有计算都基于这套参数参数数值说明晶格常数 a8 mm相邻晶胞中心的距离介质柱半径 r1.2 mm柱截面为圆形模型为无限长柱缩进参数 δ1.2 mm典型值柱中心到晶胞角点的距离胞内柱间距约 5.6 mm晶胞内部相邻柱中心距离胞间柱间距约 2.4 mm相邻晶胞最近柱中心距离柱材料相对介电常数9.8氧化铝陶瓷微波段低损耗背景材料空气εr1柱阵列周围填充空气这里的核心是缩进参数δ。当δ a/4时胞间柱间距小于胞内柱间距柱间耦合强度呈现胞间强、胞内弱的分布此时体系处于拓扑非平庸相角态出现在四个顶角附近当δ a/4时强弱耦合翻转系统进入平庸相角态消失。δ a/4是相变点此时所有最近邻距离相等系统回到普通的正方晶格。2.2 耦合强弱为什么由间距决定很多读者第一次接触这个模型时会问为什么不用两个不同半径的柱来产生强弱耦合当然可以那叫质量修饰方案也能做。但从电磁耦合的物理直觉出发两个介质柱之间的近场耦合强度主要由它们之间的间隙决定间隙越小耦合越强。所以只用一套相同半径的柱、改变柱间距就能实现需要的耦合交替。这样做还有一个额外好处柱的材料、半径保持一致能带结构更干净不易出现因为柱尺寸不同带来的杂散模式。我用紧束缚模型的语言来描述晶胞内相邻柱间的耦合强度记为t1晶胞间相邻柱间的耦合强度记为t2。当t2 t1时对应着二聚化参数为负的SSH相体系具有非平庸的拓扑性质。在二维情况下这个条件意味着每个方向上的SSH链都是拓扑相于是边与边交汇的角上会出现一个被双重保护的零维态。2.3 有限阵列尺寸怎么选仿真中的阵列尺寸是一个需要权衡的变量。阵列太小角态和边态、体态的模式数太少能谱中能级分裂严重不好判断阵列太大特征频率求解的物理内存和时间开销明显增大。我用9×9晶胞作为标准配置——总共81个晶胞、324根介质柱在Comsol中特征频率搜索频段内大概有几百个模式既能看清带隙结构又在单台工作站上跑得动。如果你想做更精细的角态劈裂分析可以扩到15×15如果只是快速验证角态是否存在7×7也够用。2.4 为什么要用TM模式体系采用无限长介质柱近似电场只有一个纵向分量Ez磁场在横平面内这就是二维TM模式。之所以选TM而不是TE主要因为介质柱阵列在TM模式下带隙更大、耦合更容易调控而且文献中大多数光子拓扑绝缘体的实验也采用这种模式便于对比验证。Comsol中设置方法很简单在二维电磁波频域物理场中选择面外电场波即在电场分量中选择面外矢量就自动只求解Ez分量。3. Comsol 5.4建模实操几何搭建、边界处理与网格策略3.1 几何构建用数组公式批量排布晶胞千万不要手动一根一根画324根柱那是灾难。Comsol的几何序列里可以直接用阵列功能或者更灵活一点把单个介质柱定义为一个部件然后使用阵列操作指定矢量偏移和副本数。我这里推荐用参数化变量控制所有几何尺寸定义全局参数a 8[mm], rc 1.2[mm], dp 1.2[mm]缩进参数单柱中心坐标以晶胞左上角为参考点四个柱的坐标分别为 (dp, dp)、(a-dp, dp)、(dp, a-dp)、(a-dp, a-dp)单位mm用阵列操作生成9×9个晶胞偏移量为 (a, 0) 和 (0, a)这样做的最大好处是后期做参数扫描时只需要扫dp一个变量整个几何自动更新。我第一次做的时候傻乎乎地手动排了5×5个柱然后想扩成9×9时在几何里改了半小时最后发现直接改数组副本数是最省事的。3.2 材料分配与介质柱损耗处理材料很简单空气和氧化铝陶瓷。在Comsol材料库中可以直接选氧化铝如果没有可以自定义空材料填入相对介电常数9.8相对磁导率1电导率设为0。这里我要特别提醒一个细节如果完全忽略介质损耗特征频率结果将是一个纯实数但实际测量中角态的Q值会受损耗影响。如果后续想和实验数据直接对比建议在材料中加一个损耗角正切比如氧化铝在微波段的tanδ约1e-4到1e-3量级。此时特征频率变成复数虚部对应损耗可以用Q Re(f)/(2|Im(f)|)来估算品质因子。但做傅里叶能谱分析时为了模式归类方便我建议先算无损耗情况把拓扑性质看清楚后再加损耗。3.3 边界条件有限阵列需要什么这里有个关键选择是直接对有限阵列做特征频率分析四周自由空间还是在四周加完美匹配层PML来模拟开放边界。我的做法是标准情况下不用PML直接在阵列外面加一圈空气区域再在空气区域最外围设置散射边界条件SBC。原因有两点第一角态本身就是束缚态场几乎完全落在阵列内部对边界条件不敏感第二SBC比PML简单计算量小。但如果你需要计算辐射损耗或者分析角态向自由空间的泄漏场那就必须用PML了。PML厚度建议设为工作波长的半个到一个波长约15-30 mm。还有一个容易被忽略的点在2D电磁波物理场中默认的外边界是理想磁导体PMC还是电导体PEC会影响求解结果。对于TM模式Ez在PEC边界上是零在PMC边界上是非零的。如果直接用默认设置不加额外边界条件然后外边界又特别靠近阵列就相当于额外加了一个金属腔会引入腔体谐振模式把能谱搞得乱七八糟。所以要么边界设远一点至少距阵列一个波长要么显式设置散射边界条件。这一步能避免大量莫名其妙的伪模式。3.4 网格划分四步走策略介质柱阵列的网格划分是我的血泪经验输出。一开始我图省事直接用物理场控制网格默认最大单元尺寸接近1.5 mm在大约10 GHz频率下空气中波长30 mm这个精度其实够了但问题在于角态模式对介质柱间隙的网格特别敏感——只有加密间隙区域网格才能准确分辨强耦合区的场梯度。我最终采用的网格策略分四步全局最大单元尺寸设为 1.2 mm约λ/25介质柱内部区域最大单元尺寸设为 0.4 mm介质柱间隙胞间距离2.4 mm的区域做边界层网格细化最大单元尺寸 0.2 mm外圈空气区域用标准三角形最大单元尺寸放大到 3 mm这样划分下来对于9×9阵列约有30-40万自由度特征频率求解一次大约1-3分钟。这个计算量在Comsol 5.4里完全可接受。注意不要一上来就网格过度细化我试过把全局最大单元压到0.3 mm自由度直接翻到近300万单次求解半小时起步精度提升却微乎其微。4. 特征频率求解能谱计算的核心设置与模式筛选4.1 研究配置特征频率搜索基准和模式数Comsol 5.4的特征频率研究配置是我花了不少时间琢磨的地方。在研究节点添加特征频率研究然后在设置中有几个关键选项期望特征频率数我建议设50到100之间。数值太少年带隙里可能漏掉角态太多则求解时间暴增且大部分是无用的边缘模式搜索基准频率设为10 GHz。Comsol会围绕这个基准值前后搜索指定数量的特征频率特征频率搜索范围可以限到8到14 GHz之间避免求解超出物理意义的极低频模式一个经验是如果只关心带隙里的角态可以先把网格粗略地算一遍全局1.5mm找到目标频率后再用细化网格做二次验证。这个粗筛-细算流程能避免大量重复计算。4.2 从几百个模式里挑出角态算法化筛选流程算完50个特征频率后你会发现列表里密密麻麻全是数值。这时候直接从数据里找角态无异于大海捞针。我的做法是写一个简单的后处理筛选流程利用角态的空间局域性做自动化标注对每个特征频率计算全局电场能量的空间分布选取4个角区域每个角的半径设为两个晶格常数分别计算该区域内电场能量占比如果某个模式的单角度能量占比超过20%标记为候选角态检查该模式的空间场图确认场强最大值确实出现在顶角而不是体区或边区在Comsol中这可以通过派生值里的体积分功能实现。虽然没有脚本那么流畅但9×9阵列面积不大逐个模式检查也能接受。如果你熟悉Comsol的Java API也可以直接用求解器脚本批量导出所有模式的电场分量到文本文件用Matlab或Python做自动筛选。这个流程我在后续第7节详细说。4.3 能谱图怎么画横轴参数纵轴频率能谱图的正确打开方式是参数扫描。保持几何和物理场设置不变把缩进参数dp作为扫描变量从0.8 mm扫到2.6 mm步长0.2 mm每个dp都求解一次特征频率。然后把所有模式的归一化频率f·a/c画在纵轴上横轴为dp或一个无量纲参数比如δ/a。这样得到的能谱图能清晰看到体带、边带和角态在参数空间的演化。我实测的典型结果是当dp较小时强胞间耦合体能带打开一个较宽的带隙带隙中央存在一条几乎平直的角态频线随着dp增大角态频线逐渐靠近边态带或体带最后在dp接近2 mm时消失进体带。这条平直频线就是角态的标志。你看文献里的能谱图时能一眼认出它。4.4 一个必须警惕的伪模式来源特征频率计算中最容易出现的一类伪模式来自计算区域的天线腔模式——如果不加散射边界条件默认PEC边界会和介质柱阵列共同构成一个微波谐振腔产生很多不属于阵列本身的腔体模式。这些模式从频率排序上看刚好可能落在带隙区域非常容易和角态混淆。判断方法很简单看场分布是否延展到计算区域的最外侧。如果场在边界的空气区域也有明显分布那基本就是伪模式。我踩过这个坑第一次用默认边界算出来的角态频率数据点和文献对不上折腾了两天才发现是边界条件导致的计算区域谐振。5. 本征场分布怎么看角态、边态和体态的区分图谱5.1 场图步骤表面图与箭头图的组合算完特征频率在结果节点添加二维绘图组选择表面图显示Ez的实部颜色范围调到对称例如从-1到1按最大场幅值归一化这样可以清楚看到场的正负区域和节点线再添加一个箭头图显示面内的磁场Hx, Hy用来判断能量流动方向。对TM模式Ez场图最直观。一个典型的角态场分布表现为场强高度集中在阵列的某个顶角附近向体区内迅速衰减。更关键的是角态场图在角两侧的边界上也会有一小段延伸但衰减比体区慢这正是一维边界态在角上截断的表现。我第一次看到这个图时那种亮斑压在角上的视觉冲击力确实很强感觉理论一下子变成了眼前可以触摸的东西。5.2 三类模式的空间指纹对比在一张能谱图上你可能同时看到体态、边态和角态。它们的空间指纹差异很大模式类型特征频率位置空间分布对阵列尺寸的依赖体态体能带内遍布整个阵列随尺寸变化不明显边态带隙中或边带沿某一条或几条边界分布随周长变化角态带隙中央孤立频点局域在顶角附近随尺寸变化但频率几乎不变指数级收敛特别注意角态的一个独特性质当阵列尺寸变大时角态频率几乎不动只是衰减长度略有变化。这是因为角态是零维束缚态其空间衰减由带隙大小决定而非阵列尺寸。这个性质可以当作一个数值实验来验证角态身份把阵列从7×7换到13×13如果某个带隙中的态频率变化极小相对于带隙宽度就说明它是拓扑保护的角态。5.3 相位图判断模式对称性除了幅值分布相位信息对角态分类也很有用。具体操作把Ez用实部、虚部两个图分别显示或者用表面图透明度叠加显示相位色相相位亮度幅值。角态的相位分布往往具有特定的旋转对称性这取决于角点所处的对称位置和晶格对称性。如果角态是四重简并的四个角各有一个简并态那么其中任意一个本征模态都是这些简并态的线性叠加场图可能同时出现在两个或四个角上。这一点在判断时要小心如果算出的一个模式在四个角上都有分布不要急着否定它是角态先看它是不是简并组合叠加的产物。小技巧如果四个角态简并可以在Comsol中用特征频率研究时选择一个频率偏移或者通过微小的几何不对称比如把某个角的柱半径改动0.01 mm来打破简并这样单个模式下就能看到局域在一个角上的场分布。但要注意这是数值破缺不是物理破缺最终结果还是要以群论分析为准。5.4 边缘态和角态频率的分离判断在能谱图上如果边态和角态频率离得很近单看场分布可能不容易分辨。这时可以做一个沿边界的能量积分曲线取一条穿过边界区域的线沿线计算电场能量密度。边态的能量密度沿边界方向缓慢起伏角态则表现为在角点处出现一个尖锐的能量峰并且往两侧边界方向指数衰减。这个方法比单纯看2D图能更定量地说明问题。6. 参数扫描与鲁棒性验证确认角态真拓扑身份的试金石6.1 拓扑相变边界的定位方法前面提到dp a/4 2 mm是相变点。在数值上你怎么观察这个相变做法是扫描不同的dp值观察能谱中带隙闭合又打开的过程。当dp接近2 mm时带隙会逐渐变窄并最终闭合跨越相变点后带隙重新打开但此时带隙中不再存在角态模式。如果你在同一个坐标系里画出dp从1.2到2.8 mm的能谱演化这个带隙闭合-重开的画面非常漂亮也很有说服力。我特别建议在扫描时把步长在相变点附近加密到0.1mm因为带隙闭合的位置对网格精度很敏感。如果你用1.5mm粗网格可能会在dp1.8mm处就观察到带隙闭合的迹象但用0.3mm细网格后相变点会回到2.0mm附近。我第一次做扫参时没有验证网格影响导致相变点偏移了0.2mm交论文前重新细化网格才纠正过来。6.2 无序扰动测试加随机位移看角态是否健在拓扑保护的一个直观检验是给结构加随机扰动看角态是否依然存在。具体操作在全局参数里加一个随机变量delta_rand给每个柱的中心坐标加一个微小的随机偏移最大幅度设为0.1mm约晶格常数的1.3%。可以在Comsol中用均匀分布随机函数或者更简单地在几何数组中对每个晶胞分别设偏移。测试结果当扰动幅度小于一定阈值时角态频率偏移很小场分布依然局域在角上而体态和边态的频率会被扰动拉宽甚至产生局域化。这就是角态比普通缺陷态硬的地方。注意随机扰动会破坏系统的平移对称性在严格的周期结构中原本简并的四个角态会分裂成一组频率相近的态这是正常现象。6.3 与文献定量对比归一化频率的对齐方法最后一步是把数值结果与陈延峰老师组或其他文献的实验数据做对比。直接对比绝对频率意义不大因为不同实验的工作频段不同正确方式是比较归一化频率f·a/c其中a是晶格常数c是光速。同时能谱中体带、边带、角态之间的相对位置是最关键的判据——角态相对带隙中心的位置略微偏上还是偏下反映了系统的具体对称性和耦合参数这个在你调整模型参数时可以精确控制。我做对比时发现由于氧化铝介电常数在不同频段的色散和实验批次差异绝对频率偏差可能达到2%-5%但归一化能谱的形状可以高质量重合。所以别死磕绝对频率拟合能谱形状才是更科学的对比策略。7. Comsol算角态的踩坑记录收敛性、模式漏算与效率优化7.1 特征频率漏算问题需要多少模式才保险Comsol特征频率求解依据设置的期望特征频率数决定计算多少模式。但如果这个数值设置太小搜索范围内的真模式会被漏掉——不是算错是根本不在输出列表里。一个稳健的原则目标频段范围内把期望模式数设为预期模式的1.5到2倍。你可以先用粗网格算一遍看目标频段内出现了多少个模式再加倍设置。另一个经常遇到的坑是不同模式在求解器中的收敛速度差异很大。体态模式通常很快收敛角态模式因为是强束缚态网格精度不足时收敛速度会变慢。如果你发现某一次求解后角态位置的频率点有明显跳跃优先减小间隙区域网格单元尺寸而不是增加模式数。7.2 内存与求解时间9×9阵列的实用配置特征频率求解是一个广义特征值问题内存消耗取决于矩阵维度。9×9阵列、30万自由度、求100个特征频率在Comsol 5.4中大约需要8-16 GB内存。我的工作配置是16 GB内存外加启用虚拟内存但更推荐的做法是减少模式数到60个左右同时利用MUMPS求解器的多核并行。如果内存仍然吃紧可以缩小一个数量级先算7×7验证确认设置无误后改用粗网格跑9×9。7.3 改变耦合参数时几何重建的稳定性扫描缩进参数dp时Comsol每次都会重新构建几何。问题来了当dp变化时介质柱之间的间隙也随之改变极小的间隙比如dp接近2mm时某些柱间距可能小到0.5mm会导致网格生成失败提示退化单元或几何无效。遇到这个问题最简单的处理方式是限制扫描范围不要扫到柱间距小于0.3mm的参数点。另外几何重建时如果出现面的自相交错误多半是因为几何模型用了联合体而不是装配体。在介质柱阵列这种重叠检测的场景中建议几何序列用并集而不是装配体避免组装时产生内边界。缺省状态下Comsol会试图创建装配体这里要手动切换。7.4 和Comsol with Matlab联用的效率玩法如果你要扫20组参数、每组看60个特征频率那在GUI里手工点扫描太慢。Comsol 5.4支持LiveLink for MATLAB可以直接用脚本控制建模、求解和导出。我常用的脚本流程是用MATLAB设置dp值循环调用model.study(std1).run()每个dp下提取所有特征频率和对应的Ez场数据将数据保存为MAT文件后续用MATLAB做能谱聚类分析这样扫20组参数大概可以做到半小时左右无人值守完成。但注意LiveLink for MATLAB需要你本地有MATLAB环境很多人会忽略License里的这一步。如果没装MATLAB也可以直接用Comsol的批处理扫描功能效率略低但也能接受。7.5 网格收敛性验证的优先级最后一条也是最重要的一条经验任何拓扑态的频率数值结果都必须在网格收敛性验证通过后才有意义。我的验证流程是选取角态频率分别用全局最大单元1.2mm、0.9mm、0.6mm、0.45mm计算观察频率收敛趋势。如果相邻两次网格细化之间频率相对变化小于0.1%就认为网格足够。角态的收敛速度通常比体态慢因为它的场压缩在极小的空间范围内梯度大。很多论文中的数值伪影或者假角态本质都是网格不足导致的数值散射效应而不是真正的拓扑局域态。这一点请务必警醒。我个人的习惯是对每一个最终呈现在论文或报告中的角态图都会把网格细化前后的频率值和场图并排放一遍确认没有任何可疑的数值痕迹。这样做虽然繁琐但在复核阶段能省掉很多麻烦。
Comsol复现二阶光子拓扑绝缘体角态：从能谱到本征场仿真
二阶光子拓扑绝缘体最迷人的地方在于几乎每个做这个方向的人第一次在仿真里看到角态场分布图时都会盯着那条局域在顶角处的亮斑看很久。角态不是凭空出现的——它在能谱里是落在带隙中的孤立频点在场分布里是集中在角落的零维束缚模式。这篇博文记录我用Comsol 5.4复现南京大学陈延峰老师组二阶光子拓扑绝缘体工作的完整过程从物理模型搭建、特征频率计算到本征场分布验证把每一步的操作思路、参数由来和坑位都摊开讲清楚。无论你是刚接触拓扑光子学的学生还是已经在用Comsol做电磁仿真但想碰一碰拓扑体系的工程师这篇文章都应该能给你省下不少摸索时间。1. 为什么要算角态二阶拓扑绝缘体的物理逻辑与仿真价值1.1 从体-边对应到体-角对应传统拓扑绝缘体最著名的特征是体-边对应一个拓扑非平庸的二维体系其体能带带有非零的陈数或Z2拓扑不变量必然导致在边界上出现单向传输的边态。而二阶拓扑绝缘体把这个对应关系往前推了一步——它不仅边界上可能有态在边界与边界的交界处也就是角上还会出现零维的角态。用一个不太严格但很好理解的说法角态是边界的边界就像一维SSH链的端态被推广到了二维的角上。我第一次看到这个概念时脑子里冒出的问题很直接这和传统的缺陷态、杂质态有什么本质区别区别在于鲁棒性。角态不是靠某个局部缺陷偶尔产生的束缚态而是受结构对称性和拓扑不变量保护的。说得再直白一点只要你没有破坏特定的离散对称性哪怕把边界切得歪一点、把两个角附近的介质柱挪一挪位置角态依然稳定存在。仿真中验证这一点最有说服力的方式就是算一个有限尺寸阵列的本征模式看带隙里有没有孤立的态再看它的场是不是压在某个角上。1.2 光子体系为什么适合做这件事光子晶体和介质谐振器阵列是实现二阶拓扑绝缘体的天然平台。和电子体系相比光子系统有几个很明显的好处第一麦克斯韦方程组本身是线性的用Comsol直接求本征频率就能得到全部模式不像电子体系要处理复杂的费米面填充问题第二微波段的介质柱阵列尺寸在厘米量级加工容易、测试方便而且场分布可以通过近场扫描直接看到第三各向同性的介质材料在微波段损耗很低仿真和实验能对上。陈延峰老师组的工作在这个方向做了很多关键推进其中一个让我印象深刻的点是他们把二阶拓扑绝缘体的角态从理论模型推到了实验测量在介质谐振器阵列里观测到了角态的本征场分布。所以用Comsol做数值复现的时候目标也很明确先把能谱算对再把场分布画出来最后对照文献里的结果确认角态频率和场形。1.3 用Comsol做这件事的合理路径Comsol在解决这类问题上有它独特的优势几何参数化方便特征频率研究配合参数扫描可以一次性把不同耦合强度下的能谱全算出来本征场直接三维可视化不需要额外导出数据再画图后处理可以做功率流、相位分析对判断角态性质很有帮助。5.4版虽然在网格剖分算法上不是最新但对于这种规模和物理场类型算力需求完全不是问题。下面我从物理模型讲起然后逐步落到Comsol的每一个设置按钮上。2. 二维SSH晶格角态模型的构造逻辑与参数依据2.1 从一维SSH到二维晶胞怎么排做二阶光子拓扑绝缘体目前最主流的可复现方案之一就是二维SSH晶格也叫三聚化格子、呼吸型格子。它的基本思路在一维SSH模型里一条链上的格点以疏-密交替排列形成强弱交替的耦合。把这个想法放到二维正方格子上每个晶胞里放4个介质柱分别靠近正方形的四个角但每个角上的柱不完全落在角点而是向内缩进一定距离。我用表格把关键参数先列出来后面所有计算都基于这套参数参数数值说明晶格常数 a8 mm相邻晶胞中心的距离介质柱半径 r1.2 mm柱截面为圆形模型为无限长柱缩进参数 δ1.2 mm典型值柱中心到晶胞角点的距离胞内柱间距约 5.6 mm晶胞内部相邻柱中心距离胞间柱间距约 2.4 mm相邻晶胞最近柱中心距离柱材料相对介电常数9.8氧化铝陶瓷微波段低损耗背景材料空气εr1柱阵列周围填充空气这里的核心是缩进参数δ。当δ a/4时胞间柱间距小于胞内柱间距柱间耦合强度呈现胞间强、胞内弱的分布此时体系处于拓扑非平庸相角态出现在四个顶角附近当δ a/4时强弱耦合翻转系统进入平庸相角态消失。δ a/4是相变点此时所有最近邻距离相等系统回到普通的正方晶格。2.2 耦合强弱为什么由间距决定很多读者第一次接触这个模型时会问为什么不用两个不同半径的柱来产生强弱耦合当然可以那叫质量修饰方案也能做。但从电磁耦合的物理直觉出发两个介质柱之间的近场耦合强度主要由它们之间的间隙决定间隙越小耦合越强。所以只用一套相同半径的柱、改变柱间距就能实现需要的耦合交替。这样做还有一个额外好处柱的材料、半径保持一致能带结构更干净不易出现因为柱尺寸不同带来的杂散模式。我用紧束缚模型的语言来描述晶胞内相邻柱间的耦合强度记为t1晶胞间相邻柱间的耦合强度记为t2。当t2 t1时对应着二聚化参数为负的SSH相体系具有非平庸的拓扑性质。在二维情况下这个条件意味着每个方向上的SSH链都是拓扑相于是边与边交汇的角上会出现一个被双重保护的零维态。2.3 有限阵列尺寸怎么选仿真中的阵列尺寸是一个需要权衡的变量。阵列太小角态和边态、体态的模式数太少能谱中能级分裂严重不好判断阵列太大特征频率求解的物理内存和时间开销明显增大。我用9×9晶胞作为标准配置——总共81个晶胞、324根介质柱在Comsol中特征频率搜索频段内大概有几百个模式既能看清带隙结构又在单台工作站上跑得动。如果你想做更精细的角态劈裂分析可以扩到15×15如果只是快速验证角态是否存在7×7也够用。2.4 为什么要用TM模式体系采用无限长介质柱近似电场只有一个纵向分量Ez磁场在横平面内这就是二维TM模式。之所以选TM而不是TE主要因为介质柱阵列在TM模式下带隙更大、耦合更容易调控而且文献中大多数光子拓扑绝缘体的实验也采用这种模式便于对比验证。Comsol中设置方法很简单在二维电磁波频域物理场中选择面外电场波即在电场分量中选择面外矢量就自动只求解Ez分量。3. Comsol 5.4建模实操几何搭建、边界处理与网格策略3.1 几何构建用数组公式批量排布晶胞千万不要手动一根一根画324根柱那是灾难。Comsol的几何序列里可以直接用阵列功能或者更灵活一点把单个介质柱定义为一个部件然后使用阵列操作指定矢量偏移和副本数。我这里推荐用参数化变量控制所有几何尺寸定义全局参数a 8[mm], rc 1.2[mm], dp 1.2[mm]缩进参数单柱中心坐标以晶胞左上角为参考点四个柱的坐标分别为 (dp, dp)、(a-dp, dp)、(dp, a-dp)、(a-dp, a-dp)单位mm用阵列操作生成9×9个晶胞偏移量为 (a, 0) 和 (0, a)这样做的最大好处是后期做参数扫描时只需要扫dp一个变量整个几何自动更新。我第一次做的时候傻乎乎地手动排了5×5个柱然后想扩成9×9时在几何里改了半小时最后发现直接改数组副本数是最省事的。3.2 材料分配与介质柱损耗处理材料很简单空气和氧化铝陶瓷。在Comsol材料库中可以直接选氧化铝如果没有可以自定义空材料填入相对介电常数9.8相对磁导率1电导率设为0。这里我要特别提醒一个细节如果完全忽略介质损耗特征频率结果将是一个纯实数但实际测量中角态的Q值会受损耗影响。如果后续想和实验数据直接对比建议在材料中加一个损耗角正切比如氧化铝在微波段的tanδ约1e-4到1e-3量级。此时特征频率变成复数虚部对应损耗可以用Q Re(f)/(2|Im(f)|)来估算品质因子。但做傅里叶能谱分析时为了模式归类方便我建议先算无损耗情况把拓扑性质看清楚后再加损耗。3.3 边界条件有限阵列需要什么这里有个关键选择是直接对有限阵列做特征频率分析四周自由空间还是在四周加完美匹配层PML来模拟开放边界。我的做法是标准情况下不用PML直接在阵列外面加一圈空气区域再在空气区域最外围设置散射边界条件SBC。原因有两点第一角态本身就是束缚态场几乎完全落在阵列内部对边界条件不敏感第二SBC比PML简单计算量小。但如果你需要计算辐射损耗或者分析角态向自由空间的泄漏场那就必须用PML了。PML厚度建议设为工作波长的半个到一个波长约15-30 mm。还有一个容易被忽略的点在2D电磁波物理场中默认的外边界是理想磁导体PMC还是电导体PEC会影响求解结果。对于TM模式Ez在PEC边界上是零在PMC边界上是非零的。如果直接用默认设置不加额外边界条件然后外边界又特别靠近阵列就相当于额外加了一个金属腔会引入腔体谐振模式把能谱搞得乱七八糟。所以要么边界设远一点至少距阵列一个波长要么显式设置散射边界条件。这一步能避免大量莫名其妙的伪模式。3.4 网格划分四步走策略介质柱阵列的网格划分是我的血泪经验输出。一开始我图省事直接用物理场控制网格默认最大单元尺寸接近1.5 mm在大约10 GHz频率下空气中波长30 mm这个精度其实够了但问题在于角态模式对介质柱间隙的网格特别敏感——只有加密间隙区域网格才能准确分辨强耦合区的场梯度。我最终采用的网格策略分四步全局最大单元尺寸设为 1.2 mm约λ/25介质柱内部区域最大单元尺寸设为 0.4 mm介质柱间隙胞间距离2.4 mm的区域做边界层网格细化最大单元尺寸 0.2 mm外圈空气区域用标准三角形最大单元尺寸放大到 3 mm这样划分下来对于9×9阵列约有30-40万自由度特征频率求解一次大约1-3分钟。这个计算量在Comsol 5.4里完全可接受。注意不要一上来就网格过度细化我试过把全局最大单元压到0.3 mm自由度直接翻到近300万单次求解半小时起步精度提升却微乎其微。4. 特征频率求解能谱计算的核心设置与模式筛选4.1 研究配置特征频率搜索基准和模式数Comsol 5.4的特征频率研究配置是我花了不少时间琢磨的地方。在研究节点添加特征频率研究然后在设置中有几个关键选项期望特征频率数我建议设50到100之间。数值太少年带隙里可能漏掉角态太多则求解时间暴增且大部分是无用的边缘模式搜索基准频率设为10 GHz。Comsol会围绕这个基准值前后搜索指定数量的特征频率特征频率搜索范围可以限到8到14 GHz之间避免求解超出物理意义的极低频模式一个经验是如果只关心带隙里的角态可以先把网格粗略地算一遍全局1.5mm找到目标频率后再用细化网格做二次验证。这个粗筛-细算流程能避免大量重复计算。4.2 从几百个模式里挑出角态算法化筛选流程算完50个特征频率后你会发现列表里密密麻麻全是数值。这时候直接从数据里找角态无异于大海捞针。我的做法是写一个简单的后处理筛选流程利用角态的空间局域性做自动化标注对每个特征频率计算全局电场能量的空间分布选取4个角区域每个角的半径设为两个晶格常数分别计算该区域内电场能量占比如果某个模式的单角度能量占比超过20%标记为候选角态检查该模式的空间场图确认场强最大值确实出现在顶角而不是体区或边区在Comsol中这可以通过派生值里的体积分功能实现。虽然没有脚本那么流畅但9×9阵列面积不大逐个模式检查也能接受。如果你熟悉Comsol的Java API也可以直接用求解器脚本批量导出所有模式的电场分量到文本文件用Matlab或Python做自动筛选。这个流程我在后续第7节详细说。4.3 能谱图怎么画横轴参数纵轴频率能谱图的正确打开方式是参数扫描。保持几何和物理场设置不变把缩进参数dp作为扫描变量从0.8 mm扫到2.6 mm步长0.2 mm每个dp都求解一次特征频率。然后把所有模式的归一化频率f·a/c画在纵轴上横轴为dp或一个无量纲参数比如δ/a。这样得到的能谱图能清晰看到体带、边带和角态在参数空间的演化。我实测的典型结果是当dp较小时强胞间耦合体能带打开一个较宽的带隙带隙中央存在一条几乎平直的角态频线随着dp增大角态频线逐渐靠近边态带或体带最后在dp接近2 mm时消失进体带。这条平直频线就是角态的标志。你看文献里的能谱图时能一眼认出它。4.4 一个必须警惕的伪模式来源特征频率计算中最容易出现的一类伪模式来自计算区域的天线腔模式——如果不加散射边界条件默认PEC边界会和介质柱阵列共同构成一个微波谐振腔产生很多不属于阵列本身的腔体模式。这些模式从频率排序上看刚好可能落在带隙区域非常容易和角态混淆。判断方法很简单看场分布是否延展到计算区域的最外侧。如果场在边界的空气区域也有明显分布那基本就是伪模式。我踩过这个坑第一次用默认边界算出来的角态频率数据点和文献对不上折腾了两天才发现是边界条件导致的计算区域谐振。5. 本征场分布怎么看角态、边态和体态的区分图谱5.1 场图步骤表面图与箭头图的组合算完特征频率在结果节点添加二维绘图组选择表面图显示Ez的实部颜色范围调到对称例如从-1到1按最大场幅值归一化这样可以清楚看到场的正负区域和节点线再添加一个箭头图显示面内的磁场Hx, Hy用来判断能量流动方向。对TM模式Ez场图最直观。一个典型的角态场分布表现为场强高度集中在阵列的某个顶角附近向体区内迅速衰减。更关键的是角态场图在角两侧的边界上也会有一小段延伸但衰减比体区慢这正是一维边界态在角上截断的表现。我第一次看到这个图时那种亮斑压在角上的视觉冲击力确实很强感觉理论一下子变成了眼前可以触摸的东西。5.2 三类模式的空间指纹对比在一张能谱图上你可能同时看到体态、边态和角态。它们的空间指纹差异很大模式类型特征频率位置空间分布对阵列尺寸的依赖体态体能带内遍布整个阵列随尺寸变化不明显边态带隙中或边带沿某一条或几条边界分布随周长变化角态带隙中央孤立频点局域在顶角附近随尺寸变化但频率几乎不变指数级收敛特别注意角态的一个独特性质当阵列尺寸变大时角态频率几乎不动只是衰减长度略有变化。这是因为角态是零维束缚态其空间衰减由带隙大小决定而非阵列尺寸。这个性质可以当作一个数值实验来验证角态身份把阵列从7×7换到13×13如果某个带隙中的态频率变化极小相对于带隙宽度就说明它是拓扑保护的角态。5.3 相位图判断模式对称性除了幅值分布相位信息对角态分类也很有用。具体操作把Ez用实部、虚部两个图分别显示或者用表面图透明度叠加显示相位色相相位亮度幅值。角态的相位分布往往具有特定的旋转对称性这取决于角点所处的对称位置和晶格对称性。如果角态是四重简并的四个角各有一个简并态那么其中任意一个本征模态都是这些简并态的线性叠加场图可能同时出现在两个或四个角上。这一点在判断时要小心如果算出的一个模式在四个角上都有分布不要急着否定它是角态先看它是不是简并组合叠加的产物。小技巧如果四个角态简并可以在Comsol中用特征频率研究时选择一个频率偏移或者通过微小的几何不对称比如把某个角的柱半径改动0.01 mm来打破简并这样单个模式下就能看到局域在一个角上的场分布。但要注意这是数值破缺不是物理破缺最终结果还是要以群论分析为准。5.4 边缘态和角态频率的分离判断在能谱图上如果边态和角态频率离得很近单看场分布可能不容易分辨。这时可以做一个沿边界的能量积分曲线取一条穿过边界区域的线沿线计算电场能量密度。边态的能量密度沿边界方向缓慢起伏角态则表现为在角点处出现一个尖锐的能量峰并且往两侧边界方向指数衰减。这个方法比单纯看2D图能更定量地说明问题。6. 参数扫描与鲁棒性验证确认角态真拓扑身份的试金石6.1 拓扑相变边界的定位方法前面提到dp a/4 2 mm是相变点。在数值上你怎么观察这个相变做法是扫描不同的dp值观察能谱中带隙闭合又打开的过程。当dp接近2 mm时带隙会逐渐变窄并最终闭合跨越相变点后带隙重新打开但此时带隙中不再存在角态模式。如果你在同一个坐标系里画出dp从1.2到2.8 mm的能谱演化这个带隙闭合-重开的画面非常漂亮也很有说服力。我特别建议在扫描时把步长在相变点附近加密到0.1mm因为带隙闭合的位置对网格精度很敏感。如果你用1.5mm粗网格可能会在dp1.8mm处就观察到带隙闭合的迹象但用0.3mm细网格后相变点会回到2.0mm附近。我第一次做扫参时没有验证网格影响导致相变点偏移了0.2mm交论文前重新细化网格才纠正过来。6.2 无序扰动测试加随机位移看角态是否健在拓扑保护的一个直观检验是给结构加随机扰动看角态是否依然存在。具体操作在全局参数里加一个随机变量delta_rand给每个柱的中心坐标加一个微小的随机偏移最大幅度设为0.1mm约晶格常数的1.3%。可以在Comsol中用均匀分布随机函数或者更简单地在几何数组中对每个晶胞分别设偏移。测试结果当扰动幅度小于一定阈值时角态频率偏移很小场分布依然局域在角上而体态和边态的频率会被扰动拉宽甚至产生局域化。这就是角态比普通缺陷态硬的地方。注意随机扰动会破坏系统的平移对称性在严格的周期结构中原本简并的四个角态会分裂成一组频率相近的态这是正常现象。6.3 与文献定量对比归一化频率的对齐方法最后一步是把数值结果与陈延峰老师组或其他文献的实验数据做对比。直接对比绝对频率意义不大因为不同实验的工作频段不同正确方式是比较归一化频率f·a/c其中a是晶格常数c是光速。同时能谱中体带、边带、角态之间的相对位置是最关键的判据——角态相对带隙中心的位置略微偏上还是偏下反映了系统的具体对称性和耦合参数这个在你调整模型参数时可以精确控制。我做对比时发现由于氧化铝介电常数在不同频段的色散和实验批次差异绝对频率偏差可能达到2%-5%但归一化能谱的形状可以高质量重合。所以别死磕绝对频率拟合能谱形状才是更科学的对比策略。7. Comsol算角态的踩坑记录收敛性、模式漏算与效率优化7.1 特征频率漏算问题需要多少模式才保险Comsol特征频率求解依据设置的期望特征频率数决定计算多少模式。但如果这个数值设置太小搜索范围内的真模式会被漏掉——不是算错是根本不在输出列表里。一个稳健的原则目标频段范围内把期望模式数设为预期模式的1.5到2倍。你可以先用粗网格算一遍看目标频段内出现了多少个模式再加倍设置。另一个经常遇到的坑是不同模式在求解器中的收敛速度差异很大。体态模式通常很快收敛角态模式因为是强束缚态网格精度不足时收敛速度会变慢。如果你发现某一次求解后角态位置的频率点有明显跳跃优先减小间隙区域网格单元尺寸而不是增加模式数。7.2 内存与求解时间9×9阵列的实用配置特征频率求解是一个广义特征值问题内存消耗取决于矩阵维度。9×9阵列、30万自由度、求100个特征频率在Comsol 5.4中大约需要8-16 GB内存。我的工作配置是16 GB内存外加启用虚拟内存但更推荐的做法是减少模式数到60个左右同时利用MUMPS求解器的多核并行。如果内存仍然吃紧可以缩小一个数量级先算7×7验证确认设置无误后改用粗网格跑9×9。7.3 改变耦合参数时几何重建的稳定性扫描缩进参数dp时Comsol每次都会重新构建几何。问题来了当dp变化时介质柱之间的间隙也随之改变极小的间隙比如dp接近2mm时某些柱间距可能小到0.5mm会导致网格生成失败提示退化单元或几何无效。遇到这个问题最简单的处理方式是限制扫描范围不要扫到柱间距小于0.3mm的参数点。另外几何重建时如果出现面的自相交错误多半是因为几何模型用了联合体而不是装配体。在介质柱阵列这种重叠检测的场景中建议几何序列用并集而不是装配体避免组装时产生内边界。缺省状态下Comsol会试图创建装配体这里要手动切换。7.4 和Comsol with Matlab联用的效率玩法如果你要扫20组参数、每组看60个特征频率那在GUI里手工点扫描太慢。Comsol 5.4支持LiveLink for MATLAB可以直接用脚本控制建模、求解和导出。我常用的脚本流程是用MATLAB设置dp值循环调用model.study(std1).run()每个dp下提取所有特征频率和对应的Ez场数据将数据保存为MAT文件后续用MATLAB做能谱聚类分析这样扫20组参数大概可以做到半小时左右无人值守完成。但注意LiveLink for MATLAB需要你本地有MATLAB环境很多人会忽略License里的这一步。如果没装MATLAB也可以直接用Comsol的批处理扫描功能效率略低但也能接受。7.5 网格收敛性验证的优先级最后一条也是最重要的一条经验任何拓扑态的频率数值结果都必须在网格收敛性验证通过后才有意义。我的验证流程是选取角态频率分别用全局最大单元1.2mm、0.9mm、0.6mm、0.45mm计算观察频率收敛趋势。如果相邻两次网格细化之间频率相对变化小于0.1%就认为网格足够。角态的收敛速度通常比体态慢因为它的场压缩在极小的空间范围内梯度大。很多论文中的数值伪影或者假角态本质都是网格不足导致的数值散射效应而不是真正的拓扑局域态。这一点请务必警醒。我个人的习惯是对每一个最终呈现在论文或报告中的角态图都会把网格细化前后的频率值和场图并排放一遍确认没有任何可疑的数值痕迹。这样做虽然繁琐但在复核阶段能省掉很多麻烦。
