这个题目我盯着看了好一会儿，第一个反应是“这哥们儿肯定也是被魔角石墨烯的图洗过脑的人”。两层结构叠一起转个角度，能带就开始不老实，这种从凝聚态物理跨界到光子学的研究思路，近两年在顶刊上确实反复出现。但真正自己动手在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 几何结构参数与材料
结构参数我直接给了固定值（这是快速验证的思路，确定物理后你再参数化扫尺寸）：
| 参数 | 值 | 说明 |
|---|---|---|
| 晶格常数 a | 500 nm | 决定带隙中心频率大致在通信波段附近 |
| 介质柱半径 r | 0.18a | 太粗带隙宽但平带效应弱，太细带隙小，0.18a比较均衡 |
| 介质柱折射率 | 3.46 | 对应高阻硅，无损耗 |
| 背景折射率 | 1.0 | 空气 |
| 层间距 d | 0.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: ky
- kx、ky作为参数，具体数值通过扫描全局参数输入
然后特征频率研究求解器会自动按你给的kx和ky来解对应的本征模式频率。这样我只需要在全局参数里定义一系列k点路径，就能依次求解能带图上的每个点。
3.2 扫描布里渊区路径：从Γ到M再到K
能带图不是全布里渊区都画，只需沿着高对称点路径跑就行。方形晶格的第一布里渊区，三个典型点分别为：
- Γ点（0,0）
- X点（π/a, 0）
- M点（π/a, π/a）
但二维光子在纯二维模型中的色散通常只沿边界上的路径来看比较有代表性，即Γ→X→M→Γ。很多文章里对正方晶格的TM模带隙，主要观察Γ到M的方向，因为带隙边缘在这个方向上更容易出现。磁性结构还要额外关注X到M的走向。我这里按照Γ → X → M → Γ的顺序扫描。
具体操作是在“研究”的“参数化扫描”里添加一个辅助函数，例如定义一个“s”来线性映射高对称点：
- s=0对应Γ点
- s=1对应X点
- s=1到2之间对应X到M
- s=2到3之间对应M到Γ
每个s区间内，kx、ky的取值给出线性路径。例如：
- Γ到X：kx=3.14159*s/Lx，ky=0（单位为1/m，但在COMSOL里直接用s乘以基础波矢）
- X到M：kx=3.14159/Lx，ky=3.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做数值验证是完全可行的。扭转光子晶体这个方向，玩的就是几何自由度，只要你愿意把角度扫描做细，很多意想不到的模式都会从能带图里冒出来。
