这个题目我盯着看了好一会儿第一个反应是“这哥们儿肯定也是被魔角石墨烯的图洗过脑的人”。两层结构叠一起转个角度能带就开始不老实这种从凝聚态物理跨界到光子学的研究思路近两年在顶刊上确实反复出现。但真正自己动手在COMSOL里复现一遍你会发现从几何建模到能带提取每一步都有暗坑和看论文图的感觉完全不是一回事。我前阵子也手贱搞了一把把整个流程从零到尾盘了一遍这篇就把我的操作路径和绕坑心得写出来给想试试手但还没敢动的人一个参考。1. 为什么“扭转”会改变光子晶体能带先搞懂物理再动手1.1 魔角石墨烯给我们的启发魔角石墨烯的事情大家多少都听过两层石墨烯叠在一起相对旋转一个特定的小角度约1.1°居然就出现了超导态。这个现象的核心物理并不神秘就是两层周期结构叠加后形成了一种更大周期的莫尔条纹也就是莫尔超晶格。这个超晶格周期远大于原始晶格常数导致能带在动量空间被折叠布里渊区变小原本活泼的电子被束缚住费米速度变慢于是产生了几乎色散为平的能带。平带意味着动能几乎为零电子的相互作用效应就被放大了超导和强关联就是从这里来的。光子晶体就是光学里的“石墨烯”。它是介电常数周期性排列的结构有光子带隙能操控光。既然电子在莫尔超晶格里会产生迷你平带光子为什么不行这个问题很自然就蹦出来了。答案是确实行。两个光子晶体平板叠在一起转一个角度后同样会形成莫尔条纹同样会折叠能带、产生平带。只不过这平带对应的不是超导而是慢光、局域态、高品质因子光学模式这些好玩的东西。1.2 莫尔超晶格与能带折叠要理解能带计算里的变化得先明白一个关键操作能带折叠。周期性结构转到一定角度后两层结构的复合周期不再是原来的晶格常数而是莫尔周期。在倒空间里这意味着原来的布里渊区缩小成了迷你布里渊区原来第一布里渊区里的能带会被“折”进这个迷你布里渊区。折叠之后原本分散在不同动量位置的模式被搬到了同一个位置如果这些模式之间有耦合就会出现能带排斥打开带隙如果耦合弱两条能带就会靠近甚至出现几乎平坦的色散。我做的模型是两层二维光子晶体板上层相对下层旋转。从直观上看转角会破坏两层之间严格的周期性匹配所以必须构建一个更大的超胞来描述复合结构。在超胞的布里渊区里做能带计算原本的能带就会因为超胞的“盖子”太小而产生大量折叠分支这些折叠分支的平坦程度直接反映了层间耦合强度——耦合越强平带越平也越容易从透射谱上看到对应的窄峰。1.3 看能带之前先想清楚你要复现什么这个项目最忌讳的就是一上来直接建双层三维平板然后闷头跑特征频率。我自己犯过这个错误结果算了三天文件几十个G啥也没看清。正确做法是先明确你要回答什么问题是想看转角从0到30度变化时平带怎么移动还是想看特定角度下的透射谱有没有对应的慢光峰目标是能带就专注能带计算目标是透射就专注端口分析。两者虽然可以共用几何但网格策略、边界条件、扫描参数完全不同混在一起做只会让计算量爆炸。我最终的策略是用二维模型先跑通物理逻辑再用对称性和参数化来筛选角度最后再考虑三维平板。二维模型虽然忽略了一些真实的垂直方向效应但对于“转角-能带-平带”这种基本物理机制来说已经完全足够用来验证思路了。透射谱部分也可以用二维做算出来的S参数趋势在定性上是可信的。2. COMSOL几何与物理场设置模型框架怎么搭2.1 选择物理场接口和研究类型COMSOL里做周期性光子晶体最常用的是RF模块下的“电磁波频域”接口。二维模型里需要选择极化方式TE模式是电场沿面内TM模式是磁场沿面内。光子晶体领域最经典的介质柱方形晶格结构带隙通常在TM模式下比较明显如果做的是孔板结构TE模式更容易出带隙。我用的介质柱方形晶格所以选TM模式也就是求解面内的电场分量。研究类型要按需求拆成两类能带计算用“特征频率”研究。扫描倒空间波矢k求解本征模式频率。透射谱用“频域”研究。设置入射端口和出射端口扫频得到S参数。这两类研究可以放在同一个模型文件里但建议分开否则每次切换研究类型时网格和求解器配置都会互相干扰反而降低效率。2.2 几何结构参数与材料结构参数我直接给了固定值这是快速验证的思路确定物理后你再参数化扫尺寸参数值说明晶格常数 a500 nm决定带隙中心频率大致在通信波段附近介质柱半径 r0.18a太粗带隙宽但平带效应弱太细带隙小0.18a比较均衡介质柱折射率3.46对应高阻硅无损耗背景折射率1.0空气层间距 d0.6a层间耦合强度靠这个调越小耦合越强介质柱的折射率和层间距是最影响结果的。旋转角度改变了柱子的相对位置层间距则直接决定两层之间光学模式的交叠程度。d太小模式强烈排斥能带分裂过大d太大两层各自为政几乎不耦合平带只是简单的能带折叠谈不上“扭转效应”。我算下来d在0.5a到1.0a之间比较合适太小了物理上接近一个单层的东西太大了就变成两个孤立层。几何构建首先要建单层晶格。先画一个晶格常数为a的正方形单元里面放一根圆截面柱体作为介质柱这是基础单元。双击“几何”节点用“阵列”功能把单元扩展成例如5×5或7×7的超胞。注意COMSOL二次开发时直接用参数化几何比较方便当扫描角度改变时可以用参数控制几何更新。我刚开时遇到一个麻烦如何在旋转时保持两层各自独立。解决方案是给每层单独设置几何域上层作为一个整体组件包含介质柱和空气背景这个组件整体绑定到“旋转”节点里通过旋转节点的角度参数来控制视角大小。这样上层是完整的一块旋转的是一个层而不是一根根柱子单独转。2.3 超胞处理公度角与非公度角的建模方案这一步是整个建模最核心也最容易让人迷糊的地方。旋转后的两个晶格要能在一个有限尺寸的模型里用周期边界条件模拟就要求整个结构严格周期性。这意味着选择的转角必须让复合结构在某个更大的周期内重复。比如转角为某特定值使得两层晶格在一个超胞内恰好重合这样的角度就是“公度角”比如30°、26.5°、21.8°等可以通过简单的几何关系构造超胞这也是常规做法。非公度角的情况更复杂无法用有限超胞严格描述。实际处理往往用大超胞近似逼近但超胞尺寸越大计算量呈指数增长而目前大部分关于魔角的研究也都是取大超胞做近似取几十纳米周期的超胞已经算很奢侈的尝试了。我这里先用了公度角快速验证选了30°它有一个相对简单的超胞构造方式而且算出来确实能看到明显的平带痕迹。如果想扫描从0到45度的连续变化可以选择在几何里设置“旋转”参数让上层旋转角度从0开始以2°步长扫描同时要在几何更新中保证每次计算都重新生成几何。假如当前超胞设定是5×5基础单元对于公度角30°在超胞适应下未必严格匹配因此我实际采用了10×10的超胞来减少边界近似带来的误差。非公度角度比如10°以内的小角度就需要更大的超胞比如30×30甚至更大这时内存消耗会非常高二维先跑通即可满足大部分物理验证需求。3. 能带计算Floquet边界、波矢扫描与模式提取3.1 Floquet周期边界条件的设置细节能带计算的基础是周期性结构所以边界条件要用Floquet周期边界在COMSOL里叫“周期”边界条件类型选“Floquet周期”。这相当于把单个超胞在空间里无限重复模拟一个真实无限的周期性结构。设置方式并不复杂把超胞的一个边界设为周期边界对面边界单独指定“Floquet周期”并添加上对应的波矢分量。关键点是波矢的传递方式——COMSOL需要你输入一个“周期矢量”默认对应kx和ky的偏移量。在二维方形晶格里如果沿x和y方向周期都是晶格常数则k矢量直接定义为一个二维分量边界两端的相位差就是exp(-ikx×a)和exp(-iky×a)这个关系在帮助文档里写明。实际设置时在“Floquet周期”功能里直接设定“波矢”参数一般来说x、y分量是固定的。我的参数设置是周期矢量x: kx周期矢量y: kykx、ky作为参数具体数值通过扫描全局参数输入然后特征频率研究求解器会自动按你给的kx和ky来解对应的本征模式频率。这样我只需要在全局参数里定义一系列k点路径就能依次求解能带图上的每个点。3.2 扫描布里渊区路径从Γ到M再到K能带图不是全布里渊区都画只需沿着高对称点路径跑就行。方形晶格的第一布里渊区三个典型点分别为Γ点0,0X点π/a, 0M点π/a, π/a但二维光子在纯二维模型中的色散通常只沿边界上的路径来看比较有代表性即Γ→X→M→Γ。很多文章里对正方晶格的TM模带隙主要观察Γ到M的方向因为带隙边缘在这个方向上更容易出现。磁性结构还要额外关注X到M的走向。我这里按照Γ → X → M → Γ的顺序扫描。具体操作是在“研究”的“参数化扫描”里添加一个辅助函数例如定义一个“s”来线性映射高对称点s0对应Γ点s1对应X点s1到2之间对应X到Ms2到3之间对应M到Γ每个s区间内kx、ky的取值给出线性路径。例如Γ到Xkx3.14159*s/Lxky0单位为1/m但在COMSOL里直接用s乘以基础波矢X到Mkx3.14159/Lxky3.14159*(s-1)/Ly这样一条连续路径就完成了。扫描时把kx和ky作为全局参数特征频率研究在每个步长下计算出若干个本征值之后再把频率值与k点对应起来画曲线。3.3 特征频率计算参数选择和模式筛选特征频率计算默认找的是最靠近某个“基准频率”的模式但这个基准频率如果不设置好求解器会抓一堆无关的低频或高阶模浪费迭代时间。我的设置是在“特征频率”研究中把所求特征频率的搜索基准设在2πc/(a×0.5)左右即大约对应晶格常数一半波长附近的频率这个区域一般就是第一光子带隙所在位置。搜索数量设为12个足够覆盖前几条能带又不至于数量太多导致计算时间指数增长。求解器选“迭代”和“直接”默认都能跑但如果模型网格大建议直接用PARDISO直接求解器内存稍微吃点但稳得多不会在特征值算法里反复跳模式。关于“丢失本征值”的一个细节COMSOL在计算特征频率时对所有本征值都有全局排序但如果在某些k点突然多出一个模式是因为在扫描范围内模式条带发生交叉或弯曲。结果里你会看到能带线之间有跳动这属于正常的模式重排不需要慌。真正要注意的是由于边界设置错误导致的假模式这种情况我会在第五节详细讲。3.4 能带图上找“平带”跑完扫描之后把特征频率画成二维图k点为横轴频率为纵轴就是能带图。我第一次看到自己算出来的有平带的能带图时第一反应是怀疑是否算错了——因为很平色散几乎为零了。这种带宽极窄的模式在色散曲线上就是一条几乎水平的直线而且经常穿插在其他陡峭能带的中间乍一看不像什么特殊模式。判断是不是物理平带有个标准方法在COMSOL里打开这个模式的电场分布图看能量是否局域在莫尔位点附近即两层柱体交叠或错位形成的“势点”区域。真正的扭转诱导平带其模场分布往往会形成清晰的周期性局域化周期对应莫尔周期而不是原始晶格周期。如果模场分布就是普通的光子晶体模式上下两层重复那就说明这只是折叠能带没用。我的经验是当转角接近某个公度角时模式场分布在超胞里会出现明显的畴结构能量集中在扭曲区域。这种模场形态和平带是互相印证的数据拿出来自己才有底气。4. 透射谱计算端口、PML与S参数提取4.1 给模型加入射端口和透射端口透射谱和能带是两套仿真思路。能带对应无限周期结构里的模式而透射谱对应有限厚度结构在外部光束照射下的响应。要做透射谱必须给模型在入射端加一个端口边界条件出射端加另一个端口边界模拟光从一端输入、从另一端输出。二维TM模式下端口设置要指定“模式”即入射光的横向空间分布。最简单的做法是用“周期性端口”功能类型选“衍射级”入射端口指定衍射级为0即垂直入射平面波如果要做斜入射对应的能带图上的不同k点可以修改衍射级对应的波矢量。我只做了垂直入射对应能带图中的Γ点附近。端口1在超胞下方端口2在超胞上方两者的参考阻抗自动设置为平面波本征阻抗。这里的关键参数是端口模式的数量不要偷懒只加一个衍射级。在超胞较大的时候高阶衍射模式在某些频点会变成传播模式对应能带里的高阶带如果不加S矩阵可能会漏掉能量导致透射反射不是1。4.2 完美匹配层与边界设置透射计算要求模拟一个“开放”的区域因此上下两端不能直接加端口然后堵上还需要在端口外侧加一段完美匹配层PML吸收向外传播的波避免边界反射干扰透射谱。PML设置有几个常见错误几何尺寸太薄导致吸收不彻底材料类型设置错误PML域必须单独指定为“PML”特征不能直接用普通材料域。在COMSOL里添加“完美匹配层”节点选择“频域”然后手动框选PML域。厚度我取约等于中心频率波长的1.5倍安全不反射。水平方向x方向依然用周期边界条件因为我们是超胞模型水平方向要模拟无限周期。对透射模型本质上还是一个周期结构的单胞做散射计算水平和垂直方向边界条件不同水平方向周期垂直方向入射/出射PML。4.3 从S参数到透射谱频域扫描完成之后COMSOL会在“全局计算”里列出S11反射和S21透射参数。但它们频率单位一般是dB也可以用abs()取线性值。要画透射率直接用20*log10(abs(S21))画dB值。我碰到的第一个反常现象是透射谱里在能带图带隙位置的透射率确实掉到-30dB以下但带通区域内也不是光滑的。扭转结构里边有大量折叠能带分支光的能量在层间来回散射谱线会出现一堆附加的峰谷结构——这本身也是物种内部的真实信息不要当bug。重点看带隙边缘附近的陡降和某些频点位置的窄峰它们对应平带和局域模式透射峰越高越细说明光在这个能量上强烈耦合进平带然后被局域。把能带和透射谱放在一起看是最说明问题的把你算出来的平带横坐标范围对应到频率轴上然后在透射谱里查这个频率点。正常情况下应能观察到透射谱中的特征峰或谷和能带的模式位置对得上。这就是扭转结构最漂亮的证据链。5. 实操踩坑记录这些坑我基本都踩过5.1 特征频率“丢根”和数值假模最烦人的问题就是能带算到一半某些k点突然少了一条带导致能带图画出来断断续续。原因是特征频率扫描时COMSOL本质上是找矩阵特征值如果两个模式在某个k点非常接近迭代求解器容易漏掉一个。解决的办法有几个我逐一试过增加搜索模式的个数比如从8个增加到15个能有效减少丢根的可能性但注意把计算时间压回来把搜索基准频率移到目标区域正中央不要把基准调到带隙中间否则特征值收敛难换用“所有特征值”搜索法而不是“围绕中心频率”搜索。还有一种更隐蔽的情况几何上出现无意的平直平移对称性比如某些高阶模式的简并度没有被COMSOL一起算出来这时候能带图上会在某个k点出现交叉。判断真假模的方法可以看模场分布真的物理模式空间分布规则、按模式序号排列清晰假模往往集中在边角或内部非物理边界处。5.2 网格密度对结果的巨大影响光子晶体仿真里网格是最影响结果的因素。我开始图省事用物理场控制的默认网格“超细化”直接开算结果能带图上平带的位置上下漂浮不定透射谱的峰形也粗糙至极。问题在于默认网格不会针对介质柱边界和层间间隙做充分加密莫尔周期本身很大而关键的光学模式局域在很薄的层间区域差一层网格结果差很多。我的最终网格策略是介质柱内部和柱体周围的空气区域用“边界层”网格加密网格最大单元尺寸设为a/40其他区域设置为a/20。对二维模型来说这个密度计算的代价不算大但精度提升是肉眼可见的。更关键的是角度变化后模型几何发生变化如果网格固定不重新生成结果直接失真。所以参数化扫描的每一步都要勾选“重新生成几何并重新划分网格”。5.3 参数化扫描转角时几何不更新的坑这是我在做角度扫描时最大的一次卡壳。直接在“参数化扫描”里给上层旋转角度φ设了一串值跑起来后频繁报错或者不报错但结果竟然完全没变。检查后发现罪魁祸首是“移动网格”和“几何更新”的不匹配我用了“移动网格”接口去旋转上层但移动网格只负责变形网格并不重新生成几何当旋转角度较大时网格完全畸变计算结果当然就错了。COMSOL里几种处理方式差异很大如果是通过“旋转”特征直接修改几何在参数化扫描里更新几何是稳妥的每次计算前COMSOL会根据新的参数重新构建几何并重新生成网格如果要用“变形几何”那性质就变了它适合小变形像转角扫描动辄十几度变形几何的网格很难跟上结果会非常不稳定。实际上对于转角扫描我推荐直接用COMSOL的全局参数控制几何的旋转角度并配合“重新划分网格”的设置而不是用移动网格功能。移动网格更适合流体、结构力学那种连续变形场景放在周期性旋转上就是自找麻烦。5.4 计算资源的平衡二维起步是理智选择直接三维起步会极其痛苦。三维的双层平板光子晶体超胞稍微大一点网格数量轻松突破千万特征频率扫描一次要好几个小时甚至直接内存耗尽。我的血泪教训是二维模型定参数、定路径、验证物理图像所有机制确认后再考虑三维做单个角度的精细化结果而不是一上来就全角度扫描。二维模型能降低计算量之外还有一个额外好处能带图中每个模式都明确区分TE或TM不会因为三维模式混合而搞不清模式对应关系。三维里TM和TE之间往往存在耦合模式不纯能带分析就复杂得多新手很容易迷失。5.5 我的调试顺序从简单到复杂给一个我验证过的心态和方法流程正好帮你避开很多弯路先做单层光子晶体算能带确认带隙位置在预期频率附近网格策略、Floquet条件全部跑通做双层0度转角对齐对比区分出“普通耦合带隙”和“扭转诱导带隙”的特征旋转一个小角度比如5度仍然用较小超胞近似确认能带里有新的折叠分支出现切换到大角度公度角比如30度精确构造超胞找平带平带确定后再拷贝一个模型副本改频域研究算透射谱验证对应频率的特征。这套顺序最大的好处是每一步的结果都能作为下一步的参照。如果你第一步就发现单层能带和你预期完全对不上千万别硬着头皮做双层回头检查参数和几何才是节省时间的方式。关于转角结构我目前的复现还停留在固定超胞近似和二维简化上后续我打算继续做的方向是把介质柱换成六角晶格结构、引入第三层对比扩展调制模式、或者干脆加增益介质看这种莫尔平带能不能把增益往带边赶。这些方向在前沿文献里都有对应的实验报道用COMSOL做数值验证是完全可行的。扭转光子晶体这个方向玩的就是几何自由度只要你愿意把角度扫描做细很多意想不到的模式都会从能带图里冒出来。