干货：声学拓扑-基于COMSOL计算拓扑不变量(1)

来源：两江科技评论

本文导读：

1，  欢迎加入香港浸会大学马冠聪老师课题组；

2，  猫哥的拓扑之路；

3，  一维模型的Zak phase计算；

4，  二维模型的陈数计算；

5，  二维valley type光子晶体的Berry Curvature；

一、欢迎加入香港浸会大学马冠聪老师课题组

模数哥按：作为COMSOL杂货店最资深的特约编辑，猫哥贡献过一些超级硬核的帖子，如声学拓扑- 2018 Nature 外尔晶体中声表面波的拓扑负折射，也指导模数哥写过：声学拓扑-2016 Nat. Phys. 基于能带反转的声拓扑绝缘体(专为猫哥写的广告帖)。在开始猫哥的雄文之前，模数哥要隆重介绍猫哥正在博后的马冠聪老师课题组，马老师不仅是科研达人，也是运动健将，现在是香港浸会大学物理系助理教授，研究兴趣包括声波超材料，声子晶体，声拓扑，非厄米系统等相关的经典波现象，欢迎各位有志青年前来学习、工作和交流。

二、猫哥的拓扑之路

要说当今物理里面最火的，还是凝聚态；凝聚态里最火的，就是拓扑。我在南大读研那会儿，是靠光子晶体、声子晶体起家，那时拓扑在凝聚态领域已经大热，基于光子晶体的拓扑也有一阵风了。当时组里面有个同学就做光子晶体拓扑，在国内已经算是很早的了。每次听他报告我都是云里雾里，诚然他做报告的水平已经是公认的很好了，奈何当时自身水平太低，也无心学业，枉费老师师兄们营造的良好环境，哈哈~也怪校园处于闹市，周边花花世界太精彩。后来放飞自我的留学生涯中某天又重拾某位物理系女神给的计算陈数的小代码，把玩一番后三月不知肉味，开始正式走上了这条不归路。

拓扑离不开边界态和拓扑不变量。如果谈可操作性的话，可以通过计算超原胞看出边界态来，但是拓扑不变量的计算就不是那么容易了。大部分情况下，人们提出的做拓扑光声的系统是对应于某个紧束缚模型，或者其简并打开这个过程可以用紧束缚模型去描述，而基于紧束缚近似可以很好地去计算拓扑不变量。另外也可以利用高对称点的态的对称性去辅助判断。所以一直以来，猫哥都很少用COMSOL直接计算，因为这确实是个大坑，难点主要在于算一遍时间久，debug难度大。但是最近越来越多的模型很难直接match到紧束缚，而其中的物理内涵也需要通过直接对COMSOL波函数进行诸如算wannier center的骚操作去揭示，所以用COMSOL算各种拓扑不变量还是很有必要的。

三、一维模型的Zak phase计算

    拓扑不变量多如牛毛。最常见的是一维有Zak phase，二维有Chern number。其算法有一定相似性，本质上都沿着闭合曲线或者曲面的积分。积分内核长这个样子：

〈ψ_n (k)|i∂_k |ψ_n(k)〉

这是个矢量函数，一维情况下的通常就是沿着某个闭合曲线做线积分。二维情况下就是∇k×作用一下得到标量函数再面积分。如果用COMSOL来计算这个的时候，就是涉及到对不同求解参数下的本征值之间的交叠积分，这里还是有很多细节的：

1．理论上波函数要归一化，声学系统光学系统里面归一化的方式有些许的差异。但是实际上我发现归一化或者不归一化结果也没影响。当然下面的计算还是归一化了的。

2． COMSOL的波函数和基于紧束缚近似写的哈密顿量h(k)的波函数有差别。

3，为了方便调用COMSOL的波函数，个人推荐comsol link with matlab，在matlab里面用COMSOL自带的函数调用模型的解。这个方法的难点在于让系统准确地识别符号语言与数值的区别，好处是报错比较快。 

其中以一维为例，武汉大学的肖孟Meng Xiao 教授2014年在港科大求学时曾有一篇PRX文章，Surface Impedance and Bulk Band Geometric Phases in One-Dimensional Systems，里面提到过一个光子晶体模型，其中第三条带的Zak phase为1，用这个模型来debug，再合适不过了。猫哥试了一把，是真π(原文Fig.1b):

该一维光子晶体的第三条带的Zak phase 为π.

如果仅有一个”π”字，在文章图里即使标得再大，也不令人信服。因为有很多种方法可以推断出这个带的Zak phase是π，如果你直接跟审稿人说是用COMSOL算的，审稿人估计也不太信，所以我们得掏出点更可视化的数据。

四、二维模型的陈数计算

下面再讨论一下Chern 绝缘体的陈数计算。其实Chern 绝缘体的陈数，可以有好几种算法，比如用类似于前面一维案例的操作。首先我们定义两个可以覆盖布里渊区的矢量ka, kb，对ka上的每个点，算沿着kb的积分可以得到一个数，用这个数对ka画图，就会出来很有物理味道的结果，把它放在文章里面，就可以实锤地说明咱们确实算过，没有偷懒。Zheng Wang教授和Y. D. Chong教授2008年发表在PRL上工作Reflection-Free One-Way Edge Modes in a Gyromagnetic Photonic Crystal，是光学拓扑的milestone [模数哥按：咱们最近也有重复Y. D. Chong教授早年工作：光学拓扑-2013 PRL 基于环阵列耦合的拓扑绝缘体]，以这篇文章的二维光子晶体为例，给出如下结果：

原文Fig.2二维YIG光子晶体，其第二条带与第三条带的Wilson loops.  

这样一来，内行人就很清楚了，我们可以计算其每个固定k的“Zak phase”,从结果来看斜率和陈数正负有关系，绕圈多少（突变）则代表了陈数是多少。这个算法似乎可以总结为研究Berry curvature的线积分。根据定义，陈数也可以通过Berry curvature的面积分计算。离散计算Berry curvature的方法也不是那么直观，有篇日本人写的参考文献被引用了500多次。再一番小操作。可以得到这两个图，积分积一下，可以发现结果恰好是1和-2. 笔者用的离散点是20乘以20，最后得到的结果是1.0000与-2.0000，精度惊人，说明这个离散计算还是收敛很快的。

第二条带与第三条带的Berry Curvature.其积分分别为1和-2.

五、二维valley type光子晶体的Berry Curvature

   二维系统中也有很多陈数为0的拓扑相，比如valley。Valley的总陈数为0，但是在k点或者K‘点有非零的berry curvature分布，而且这两个分布恰好相反，如果将k点或者k‘点附近的Berry curvature积分，结果就是1/2或者-1/2. 下图就是一个简单的通过打破空间反演对称性实现的valley type的光子晶体。

简单的二维石墨烯结构光子晶体，两个sub lattice具有不同半径的圆柱，其第二条带Berry Curvature在K与K‘点分别具有非0的Berry curvature.