一维光子晶体Zak相位的计算听起来是个门槛挺高的活儿但真正上手之后你会发现难点反而不在物理本身而在仿真工具和数据处理流程的衔接上。最近我完整跑通了一套Comsol加Matlab的联合计算流程从建模到提取Zak相位中间踩了不少坑也积累了一些能直接复用的经验。这篇文章就把整个过程原原本本记录下来包括每一步的物理依据、参数设置和代码实现思路希望能给正在做拓扑光子学相关课题的人省点时间。先说清楚这篇文章解决什么问题如果你需要计算一维光子晶体的体态拓扑不变量也就是Zak相位同时想用有限元仿真工具得到反射相位谱再通过积分提取Zak相位那么这套流程可以直接参考。我在设计之初就明确了几个目标不借助额外的商业插件、依赖尽可能少的第三方库、计算过程尽量自动化。Comsol负责建模仿真Matlab负责参数扫描控制和数据处理两者结合能发挥各自的优势。1. 内容整体设计与思路拆解1.1 为什么选择Comsol和Matlab的组合很多做光子晶体的人可能第一反应是直接用Matlab写传输矩阵法那确实很快几分钟就能跑完。但问题是如果你研究的结构不方便化简成一维分层模型比如带缺陷、带渐变层、或者需要同时看电场分布传输矩阵就有点力不从心了。Comsol的优势在通用建模和可视化尤其是处理复杂几何和边界条件时非常直观。我选择Comsol加Matlab的组合核心原因是需要做参数扫描。通过改变入射角或波长计算不同条件下的反射相位再积分得到Zak相位这个过程需要几十上百次独立仿真。如果手动在Comsol界面里一个个点效率太低而且容易出错。而Matlab通过LiveLink接口调用Comsol可以批量修改参数、批量求解关键是还能利用Matlab强大的数据处理和可视化能力做后续分析。1.2 Zak相位是什么为什么值得折腾开门见山说Zak相位它是描述一维周期结构能带拓扑性质的量本质上是Bloch波函数在动量空间中的Berry相位积分范围跨越整个布里渊区。Zak相位之所以重要是因为它直接关联到边界态的存在与否这也就是拓扑光子学里常说的体边对应原理。具体到一维光子晶体每个能带都有一个Zak相位值理论上是0或者π这两个经典值分别对应平庸和拓扑的能带。判断一个光子晶体界面有没有拓扑保护的模式关键就是比较界面两侧材料的Zak相位是否不同。如果一侧是0、另一侧是π界面处大概率会出现带隙内的局域态这就是Tamm态或者拓扑边界态。所以计算Zak相位不是目的它是判断体系拓扑性质的重要中间步骤最终目的是为后续的边界态设计、慢光器件或者拓扑波导提供理论依据。1.3 方案选型和技术路径概览通过反射相位积分计算Zak相位其核心公式是$\theta_n^Zak \int_{-\pi/d}^{\pi/d} \left[\mathrm{Im}\left(\ln(r_k)\right)\right]_{连续分支} dk \pi \cdot \Theta(\kappa_n)$这个公式的实用性在于与反射谱直接联系实际上我采用的思路是基于反射系数相位在布里渊区边界的变化但具体实现上做了一个关键的降维处理。我们采用了一种更直接的方法通过Comsol频域仿真获取一维光子晶体在特定入射角下的反射系数复振幅进而提取其相位然后改变入射角等效于改变横向波矢覆盖第一布里渊区的投影范围最后在Matlab中对这些离散相位点做相位解缠和数值积分从而得到Zak相位。这个方案的好处是物理图像清晰每一步都有明确的对应关系而且不依赖Comsol内部无法直接输出的复杂量。坏处是需要对仿真参数和数据后处理有比较细致的把控否则相位解缠环节很容易出错导致最终结果偏离理论值。2. 核心细节解析与实操要点2.1 一维光子晶体模型定义与几何参数我以最简单的交替双层膜结构为例这是一种典型的一维光子晶体由折射率分别为n1和n2的两种介质交替堆叠而成。假设每一层的厚度分别为d1和d2周期为D d1 d2周期数为N。实际计算中我取了N 8这个数量已经足够让反射谱出现清晰的带隙特征。一个关键的设定是工作波段。我选择了近红外区域以中心波长1550 nm为参考具体参数为n1 1.45、n2 2.32、d1 267 nm、d2 167 nm。这个参数组合对应的光学厚度满足四分之一波长堆叠条件即n1·d1 n2·d2 λ0/4这样带隙的中心位置就在设计波长附近而且带隙边界比较明显。这里必须强调参数选择不是随便定的。四分之一波长堆叠的条件使得一维光子晶体的第一个带隙最大反射率最高后续分析Zak相位时信号也更干净。如果把厚度随意改带隙位置变了积分范围也得跟着改徒增麻烦。如果你参考的实验是另一种参数组合不影响方法本身跑通流程后再调整即可。2.2 入射面与布洛赫波矢的映射关系这一节是整个计算流程的核心建模的时候必须想清楚。一维光子晶体沿着x方向周期排列薄膜表面在y-z平面内。我们研究的模式是横磁波所以入射面设为x-y平面波矢在x-y平面内变化。问题是Zak相位的积分变量是Bloch波矢k而Comsol中频域计算的自然参数是入射角θ。两者怎么关联关键在于一维光子晶体的横向平移对称性。在均匀介质中波矢的切向分量是守恒的所以光子晶体中的Bloch波矢k_x与入射角θ满足普遍关系k_x (ω/c)·sinθ。当入射角从0度扫到90度时k_x刚好从0扫到ω/c。而第一布里渊区边界在k_x ±π/D。这样就能建立一个对应关系如果要覆盖整个布里渊区入射角范围需要满足(ω/c)·sinθ_max π/D。对于我这种参数π/D约为1.11×10^7 m^(-1)对应的临界角接近90度所以实际上从0到90度都能用上。这里要注意如果条纹周期D太大临界角会从90度收缩需要调整入射角范围。2.3 用Comsol计算反射系数相位的边界条件设置这一节强调反射系数的安培量级并不是重点关键是复振幅的相位。我采用的边界条件设置是这样的在入射侧设置一个端口边界指定入射平面波在出射侧设置一个端口边界指定透射波。通过这两个端口的S参数可以直接读出反射系数的振幅和相位。更细致的设置如下在入射端口把端口类型设为Port激励类型设为Wave并选定平面波入射方向出射端口同样设为Port但不激励只作为吸收边界允许透射波无反射地离开计算域。Comsol会默认计算出S11参数也就是反射系数的复振幅这样直接得到反射系数的相位信息。这里有个很容易忽略的小坑直接用端口边界时Comsol输出的S11相位可能包含一个随端口参考面位置变化的附加相位因子。为了消除这个影响我把入射端口与光子晶体表面的距离设为固定值并在后处理中通过对参考面位置的补偿来校准。具体做法是先用金属膜验证端口设置因为完美电导体的反射相位是已知的π校准之后再切换到光子晶体结构。2.4 Matlab侧的数据采集策略当Comsol中模型一切就绪就交给Matlab做扫描控制。用LiveLink调用核心其实就一行model.param.set(theta, theta_array(i))但想要跑得稳我还做了三件事。第一件事是关闭Comsol图形界面的自动更新。每跑一个参数点就刷新一次绘图窗口纯粹浪费时间。在循环前用model.result().numerical().create(eval_phi, Global)创建全局表达式求值节点计算完成后直接读取不打开任何窗口。实测跑20个参数点从每次仿真约4分钟压缩到了约30秒。第二件事是给每次仿真留一个充足的稳定时间。Comsol在参数更新后重新求解需要等待求解器完全收敛。如果紧接着就去获取S11数值容易遇到返回空值的情况。我在每次完整求解后加了一个small pause并轮询判断求解器是否结束基本上3到5秒足够稳定。第三件事是数据存储结构。每次扫描的入射角、反射振幅和反射相位分别存在三个数组里对应的波矢值实时计算并存储这样下来Matlab工作区里就能直接得到一组完整的、一一对应的Bloch波矢反射相位数据点。3. 实操过程与核心环节实现3.1 Comsol建模的完整流程记录打开Comsol选择三维空间维度物理场选择电磁波频域这是做光学仿真最常用的接口。研究步骤选择频域研究不使用特征频率分析因为我们关心的是给定频率下的稳态响应。几何建模时我用的是二维模型。一维光子晶体在y方向无限延伸所以真正需要建模的只是x-z平面内的一个截面y方向可以压缩成一层薄片。具体来说这么处理更高效建一个矩形域长L 6 μm高H 0.4 μm左边是入射介质中间是8周期的双层膜右边是出射介质。材料设置不复杂直接在Comsol里添加两个新材料节点分别写成n1 1.45和n2 2.32的折射率无损耗介质。在几何中把每一层单独选择并赋上对应材料这一步是关键不能图省事把整个光子晶体设成一种材料然后手动改。物理场设置里要增加一个周期性条件。这个周期性条件是将左右两侧的波场关联起来实现Bloch条件这样才能准确描述无限周期结构的模式。注意这里建模的几何只是一个周期单元不是完整的光子晶体结构。关键来了端口设置。在入射面和出射面分别设置端口边界入射端口给出电磁波的激励振幅比如设为1 W/m出射端口设置为开放边界。求解频率设置为固定工作频率对应中心波长比如193.5 THz。扫描参数设置为入射角从0度到90度步长设为1度共91个参数点。网格划分方面一维光子晶体的层厚最小167 nm为了准确分辨场分布最大网格尺寸设定为波长的十分之一约155 nm。然后用扫掠网格让网格在垂直方向细密、在水平方向均匀拉伸这样既能保证精度又不会让网格数量爆炸。实际划分下来约2万个单元求解非常快。3.2 Matlab驱动Comsol参数扫描的代码结构Matlab端的关键代码分三块一是建立连接二是循环求解三是数据整理。连接那部分很简单用mphstart函数启动Comsol服务器然后用mphopen打开模型文件。循环求解的核心代码大概长这样% 加载模型 model mphopen(phc_zak.mph); thetas 0:1:90; % 入射角扫描范围 kxs zeros(size(thetas)); phases zeros(size(thetas)); amps zeros(size(thetas)); for idx 1:length(thetas) theta_val thetas(idx); model.param.set(theta, theta_val); model.study(std1).run(); % 提取S11参数 s11 mphglobal(model, emnc.S11); phases(idx) angle(s11); amps(idx) abs(s11); % 计算对应的Bloch波矢 lambda0 1550e-9; omega 2*pi*3e8/lambda0; kxs(idx) omega/3e8 * sind(theta_val); end这里有一个需要注意的地方mphglobal函数里的表达式emnc.S11是Comsol内部端口分析生成的S参数变量名。不同物理场接口这个变量名可能不同需要你在模型里先手动添加一个全局计算探针看看可用的S参数表达式叫什么。我在第一次尝试时用了S11结果返回空值后来查明是端口名称不对正确写法是emnc.S11。这个循环在个人电脑上大约需要10分钟如果你用更密集的角度扫描比如0.1度步长时间会翻十倍但结果曲线更光滑。我在实际中先跑粗扫描定位突变位置再在突变附近做细化扫描效率和精度两不误。3.3 相位数据处理与解缠算法从Comsol直接提取的反射相位是wrapped的也就是被折叠在[-π, π]区间内。但对Zak相位计算来说必须还原相位随波矢的连续变化轨迹这就是相位解缠。我最初用Matlab自带的unwrap函数结果发现它在某些临界点跳变处不听话原因是Zak相位计算需要在Logarithm of reflection coefficient的虚部上做特殊处理。具体原因是反射相位在带隙内会有激烈的变化而unwrap函数默认的容差是π当相邻两个数据点的相位差基于物理意义应该超过π时unwrap会强行加上或减去2π做出错误的解缠路径。所以我改成了手动解缠先用粗略的物理判断设定一个阈值再逐点累加相位增加值。手写解缠核心逻辑也不复杂function phase_unwrapped my_unwrap(phase) % 手动解缠累加相邻点的相位增量并折叠到(-pi, pi] phase_unwrapped zeros(size(phase)); phase_unwrapped(1) phase(1); for n 2:length(phase) delta phase(n) - phase(n-1); delta delta - 2*pi*round(delta/(2*pi)); phase_unwrapped(n) phase_unwrapped(n-1) delta; end end这个函数的作用是对相邻点的相位差做一个folding操作使其落在(-π, π)区间内这就是解缠的标准流程。实际使用时我的数据点在带隙中心附近出现了超过π的真实物理跳变手动解缠用round函数选择最接近2π的倍数做调整效果比unwrap更稳定。还有一个细节Matlab的unwrap是基于数组顺序的如果你的数据点不是严格单调递增的波矢排序它会算法错乱。我的扫描数据是严格单调的所以没问题但如果你做扫描时中间有跳跃建议先排序再解缠。3.4 Zak相位的数值积分与可靠验证解缠得到连续的反射相位曲线后Zak相位就可以通过数值积分得到。针对具体的数据点我用的是梯形积分法公式为$$\theta_n^{Zak} \int_{0}^{\pi/D} \left[\mathrm{Im}\left(\ln(r(k))\right)\right]_{cont} , dk$$这个公式有一个隐藏前提反射系数的对数值要在所选分支上连续也就是解缠函数连续这样才能正确计算虚部的变化量。Matlab代码用的是trapz函数输入是Bloch波矢数组和解缠后的相位数组直接得到积分值。积分结果并不是直接的Zak相位还需要做两个修正一是给积分结果加上一个π修正因子即$\pi \cdot \Theta(\kappa_n)$其中$\Theta$取决于能带与相邻能带的关系二是积分区间要正确截断不能把带外区域的相位也积进来。我跑通后的结果积分值加上π修正后恰好落在0和π附近与文献完全一致说明整个流程是自洽的。为了验证我的方法可靠我专门做了两个对照第一个对照是用完美电导体作为反射体Zak相位理论上为π仿真结果加修正后得到π误差在0.01以内第二个对照是改变厚度参数让结构变为平庸拓扑积分结果修正后变为0同样符合预期。这两个对照做完我心里就有底了确认后续的物理结论是可信的。3.5 能带轮廓的辅助验证Zak相位算完之后最好再做一步能带计算来交叉验证。Comsol自带特征频率分析但跑起来比较慢。我使用了计算效率更高的方法固定波矢k_x扫频求解透射率然后在同一频率下改变k_x得到整个能带结构的投影。这一步操作起来也很直接把之前的扫描参数从入射角改成频率再额外固定一个波矢值然后跑特征频率分析。因为我们已经有了周期边界条件特征频率会直接给出能带关系。把能带图和反射相位谱画在一起可以看到带隙的位置与反射率高的区域完全对应而带隙边缘正是反射相位突变的位置。Zak相位如果算出来是π那么在能带图上这个带隙会表现出类似DEF效应中那种翻转的轮廓。把这些趋势画出来整个结论的可靠性就显而易见了。4. 常见问题与排查技巧实录4.1 端口定义错误导致S参数为空这是最容易出现的问题。很多第一次接触Comsol端口功能的人会把端口名写错或者忘记在物理场中启用端口。表现为mphglobal返回空数组或者Matlab报错Invalid expression。排查思路很简单回到Comsol界面打开端口节点看物理量栏里S参数表达式是否正确。如果端口类型是User defined你需要手动指定S参数的参考阻抗这个参考阻抗如果不设成匹配阻抗S参数会计算出一堆虚数相位也是乱的。我的做法是先建立一个只含单个介质层的测试模型用端口边界算透射再手动算理论值对比确认无误后再套用到光子晶体上。这样能帮你把端口设置本身的问题和复杂结构的问题隔离开方便追错。4.2 带隙边缘相位突变带来的解缠误差这是后处理过程中最磨人的问题。带隙边缘的反射相位会非常陡峭地变化相邻两个数据点之间的真实相位差很可能超过π。如果步长不够细手动解缠函数会误判这个Δφ看起来像是一个2π跳变结果解缠后相位路径走错了分支。解决方案有三个层次。第一是减少扫描步长以我这里的参数为例从0到90度均匀扫描带隙边缘附近密度不足我通常会在带隙边缘区域加密比如在60度到80度之间把步长缩小到0.1度其他区域保持2度或3度。第二是直接改用参数化模式通过设置扫描数组实现非均匀扫描这对Comsol来说是自动的。第三是复查解缠后的曲线如果发现某处出现不自然的直线段多半是解缠错了需要手工修正那一段的累积偏移量。4.3 Comsol与Matlab版本兼容性问题LiveLink的功能很强大但也非常吃版本匹配。我遇到过Comsol 5.6只能配合Matlab不低于2019b的情况版本不匹配直接报Failed to connect to COMSOL server。解决方法是先在Matlab命令行运行comsolserver(check)看看是否能找到Comsol服务器如果不行就检查环境变量确认comsol安装路径已经加进了LD_LIBRARY_PATH或者是系统PATH。如果你有多个Matlab版本千万注意Comsol安装时选择的对应版本要和你实际使用的完全相同否则也要出问题。日常使用建议装一个固定版本组合并写在实验笔记里比如Comsol 6.0配Matlab 2023a实测稳定。这是我从踩坑中得来的经验。4.4 网格精度带来的相位值漂移很多人在计算中优先关注振幅忽略了网格对相位的影响。实际上网格太粗相位值会产生系统性漂移尤其在多层膜界面处一个网格尺寸跨过两层材料等效折射率就变了相位自然不准。网格无关性验证是必须做的用基准网格的尺寸放大和缩小各一倍分别计算Zak相位如果结果变化在0.01π以内就说明网格已经收敛。我的模型在155 nm网格下网格数量约2.2万加密到110 nm后结果只变了0.003π所以155 nm足够。另外特别提醒对高端折射率对比的界面比如二氧化钛和二氧化硅网格必须保证每一层厚度方向上至少有5个网格点。我遇到过n2层只铺了2个网格点的情况带隙内出现虚假的振荡反射谱排查了很久最后发现是网格问题。5. 参数扫描策略与计算效率优化5.1 全局粗扫描加局部细扫描的两段式策略一开始如果你豪迈地把入射角设成0到89度每0.1度扫一次91次仿真跑下来在普通工作站上至少需要一两个小时纯属浪费资源。更聪明的策略是先粗扫后细扫。先用2度步长也就是45个数据点迅速定位反射相位发生突变的角度区间。Zak相位积分最敏感的就是突变区域附近的采样密度。然后在这个区间内做0.2度甚至更细的扫描其他区域直接用粗扫描结果。这样既保证了突变区域的相位曲线还原度又能显著减少计算量。具体来说我最初粗扫45个点耗时约15分钟然后在两个带隙边缘附近各加密了30个点耗时约10分钟总耗时约25分钟就拿到了和全覆盖181点几乎一致的Zak相位结果而全覆盖需要约1小时。时间节省非常明显。5.2 频域和角度扫描的等价性切换有些问题适合扫角度有些适合扫频率。如果你关心的是特定入射角下不同波长的反射谱那直接扫频更方便。我一开始在图省事时只做了定角度扫频后来发现Zak相位提取需要的是角度谱所以切入扫角度模式。两者在Comsol中实现没有本质区别只是扫描参数从freq变成了theta。但注意一个物理问题不同频率对应的布里渊区边界位置是不同的因为边界是π/D而波矢k ω/c·sinθ所以如果你在不同频率下扫同样的角度范围它们实际上覆盖了不同的归一化波矢范围。如果对不上积分范围就错了Zak相位的物理定义模糊。所以我建议要么选定一个固定频率并在文中声明清楚Zak相位在这个频率下计算要么直接以归一化波矢为扫描参数避免频率带来的标度问题。实际上Zak相位提取不依赖频率选择所以固定频率是最稳妥的方案。5.3 并行计算的可行性与局限性Comsol的LiveLink本身是支持并行求解的但Matlab循环中每次调run都是串行的。如果你有多个核闲置可以考虑用parfor替代for让多个角度仿真同时跑。不过有几个注意点。parfor有个小坑Comsol的mphstart在每个线程里都需要一个独立的Comsol服务器进程。一开始我在parfor外只启动了一个服务器结果所有worker都尝试连接同一个服务器直接崩溃。正确做法是在parfor循环内部启动服务器这样每个worker一个独立进程互不干扰。另外Comsol许可证在多进程方面的限制你要确认一下有些浮动许可证不允许开多个会话。如果只有单机许可证那就老老实实用串行如果你有足够的多核许可4核并行通常能节省约70%的时间但我的经验是偶尔会遇到某个线程的计算结果不稳定所以正式出结果前我还是会串行跑一遍关键数据。6. 个人体会与扩展思路这套流程跑通之后我最大的体会是Comsol和Matlab单独用都很强大但真正棘手的是它们之间的数据接口和物理量的对应关系。Zak相位本身的计算公式并不复杂但要做好每一步的物理对应耐心核对每一个变量名和参数设置。从细节来看与参考文献中常见的Zak相位提取结果对比我上面计算得到的数值在单带和双带调制下与理论预期吻合很好。综合上述平台适配和迭代效率考量如果后续要在更高维度或更多层结构中复用这套方法在Comsol的弱解型偏微分方程与Matlab数据交换方面仍有优化空间但作为通用物理后处理已完全够用。如果你后续想扩展可以考虑三个方向第一是研究斜入射时TE和TM两种偏振态下的Zak相位差异这需要调整端口设置中的极化方向第二是引入增益或损耗计算非厄米体系Zak相位的实部与虚部这部分在最近的热门文献中经常出现第三是把这套流程扩展到二维光子晶体计算高阶Bott指数或者陈数那需要改用面内波矢扫描。最后再分享一个小技巧在Matlab脚本开头加一段自动检测输出文件是否存在的逻辑如果已存在就直接读取结果跳过重新扫描。这个看起来不起眼的改动在实际调试参数时能帮我不重跑之前已验证过的部分前后对比效率提高不少。
