动手做二维SSH模型之前我把一维SSH的论文来回翻了好几遍。一维的情况很清爽胞内hopping一个参数胞间一个参数边界态、Zak phase、末端极化全都讲得明明白白。可是到了二维论文里的图变成了一张张原子排列图和能带色散图模型定义散落在附录里画图脚本通常不公开。我当时的想法很朴素能不能用Matlab把二维SSH模型从最底层的紧束缚哈密顿量开始一步步算出原胞能带再算出带表面权重的投影能带。这篇文章就是那次复现过程的完整记录包括代码骨架、参数约定、投影权重的定义方式以及我踩过的几个能让人结果错得莫名其妙的坑。适合已经会跑一维SSH、想往二维拓扑模型迈一步的人如果你只是想找个能练手的紧束缚模型把Matlab的矩阵构造和本征值求解流程摸熟这篇也够用。1. 二维SSH模型到底在算什么先有物理图景再动手写代码1.1 一维SSH的边界态与拓扑数二维推广的地基二维SSH不是凭空冒出来的模型。它的根在一维SSH链一个原胞两个格点分别记为A、B胞内hopping是v胞间hopping是w。在动量空间里哈密顿量可以写成H(k) [0, v we^{-ik}; v we^{ik}, 0]能量为E(k) ±√(v² w² 2vw cos k)。当v≠w时能谱打开一个gap。有限长度的链放在这个参数范围里如果|v/w| 1会在一端和另一端各出现一个零能的局域态这就是SSH边界态。拓扑上的表述是winding number非零Zak phase等于π。这套东西几乎所有拓扑物理的入门笔记都会讲但它真正有用的地方是给出了一个思维模板**体态的能带gap本身不包含边界态信息边界态是有限尺寸系统里才出现的额外自由度需要用开边界或者半无限边界的计算去捕捉。**二维SSH要解决的核心问题正是把这种“体态能带”和“边界态诊断”分离开来。1.2 二维SSH原胞怎么选四格点排列与hopping对称性二维SSH模型在文献里并不是唯一写法。最常用的做法是把一维的二聚化模式同时用到x和y两个方向得到一个正方形原胞里面放四个site。我采用的是这种四格点排列方式四个site分别放在原胞的左下、右下、左上、右上基矢取a1 (1,0)a2 (0,1)格点坐标为site1 (0,0)site2 (0.5,0)site3 (0,0.5)site4 (0.5,0.5)。四条边分别对应胞内hopping水平两条边上两个site之间的耦合为v_x垂直两条边的耦合为v_y。跨原胞的方向上site2向右与右侧原胞的site1耦合site4向右与右侧原胞的site3耦合振幅为w_xsite3向上与上方原胞的site1耦合site4向上与上方原胞的site2耦合振幅为w_y。为了减少参数让结果好看我建议前期先取各向同性v_x v_y vw_x w_y w。这样模型只保留两个参数v和w的比值决定相图。vw对应拓扑非平庸的区间vw是平庸区间两者中间发生能带gap闭合再打开。这个模型的漂亮之处在于它不仅能展示原胞能带和投影能带的差别还能顺带展示二维体系里比一维推广更有意思的现象不止有边缘态还可能出现四角束缚态也就是所谓的高阶拓扑角态。第5章会专门说这个。1.3 先给结论原胞能带和投影能带分别看什么原胞能带或者说bulk能带是把二维晶格当作无限周期系统玻尔兹曼变换后得到H(kx,ky)在布里渊区高对称线上扫描本征值。它回答的问题是体态是金属还是绝缘体gap在哪里色散关系长什么样。投影能带则完全不同它回答的问题是如果在这个方向切出边界gap里会不会出现边界态。做法是沿着某一个方向保留有限个原胞另一个方向继续保持周期边界然后计算能谱再按本征向量在表面层的振幅加权让局域在边界上的态在图上“变亮”。这就是标题里两个能带的实质分工。接下来我从紧束缚模型出发先把代码层面的事说清楚。2. 紧束缚哈密顿量的Matlab构造从实空间参数到H(k)2.1 先列bond清单把模型约定写清楚写代码之前我习惯先把所有hopping列成一张表避免原胞内外一堆耦合混淆。本文约定的bond清单如下bond说明振幅跨原胞方向site1 → site2原胞内水平边v_x无site3 → site4原胞内水平边v_x无site1 → site3原胞内垂直边v_y无site2 → site4原胞内垂直边v_y无site2 → 右侧原胞site1跨原胞水平w_xxsite4 → 右侧原胞site3跨原胞水平w_xxsite3 → 上方原胞site1跨原胞垂直w_yysite4 → 上方原胞site2跨原胞垂直w_yy原胞内的四个site并没有斜对角耦合site1和site4、site2和site3之间没有hopping。这一点很重要很多人第一次画二维SSH时会顺手把四方原胞的对角也连上那得到的模型就不是这个模型了拓扑性质会变。2.2 直接从H(k)矩阵入手原胞能带需要4×4哈密顿量如果只算完整的二维周期系统不需要去构造大矩阵。对每个需要计算的(kx, ky)直接写出来这个4×4的哈密顿量Hk(1,2) vx wx * exp(-1i*kx); Hk(3,4) vx wx * exp(-1i*kx); Hk(1,3) vy wy * exp(-1i*ky); Hk(2,4) vy wy * exp(-1i*ky);然后补上厄米共轭部分让矩阵变成Hermitian矩阵。用Matlab写的话核心函数可以这样组织function Hk ssh2d_Hk(kx, ky, v, w) % site order: 1(0,0), 2(0.5,0), 3(0,0.5), 4(0.5,0.5) Hk zeros(4,4); Hk(1,2) v w * exp(-1i*kx); Hk(3,4) v w * exp(-1i*kx); Hk(1,3) v w * exp(-1i*ky); Hk(2,4) v w * exp(-1i*ky); Hk Hk Hk; end这里有个相位符号的约定问题。exp里的符号取正还是取负取决于你定义的跳跃方向和布洛赫变换约定但对能带数值本身没有影响因为物理上E(k)E(-k)。真正要紧的是Hk必须是厄米的否则本征值会出现虚部。我建议所有代码都用Hk Hk Hk这种方式补全另一半这样只需要写下三角的一半不容易写错。2.3 为什么slab哈密顿量不能照搬H(k)方向上的相位叠加投影能带需要构造slab也就是沿着x方向取Nx个原胞、开边界沿着y方向保持周期边界、用ky参数化。很多人第一次写这段代码时会试图把前面的4×4矩阵按原胞数复制Nx份然后在原胞之间填w_x的耦合。这个思路本身没错但很容易漏掉一个关键细节在y方向已经是周期的前提下同一个site索引之间可能同时存在胞内hopping和跨原胞hopping的贡献必须叠加起来。以site1和site3为例。site1在(0,0)site3在(0,0.5)它们之间有一个胞内的垂直hopping v_y。同时site3向上跳半个原胞就到达上方原胞的site1这个跨原胞hopping是w_y。在固定ky之后这两个跳跃都出现在矩阵元H(site1, site3)里区别只是第二个跳跃多了一个相位因子exp(i*ky)。所以slab矩阵里这两个跳跃要写在同一个矩阵元上而不是分别写在两个不同的原胞块里。这也是滑块法和普通实空间大矩阵法的一个本质区别只要某一方向是周期边界那一方向上的跨原胞跳跃就会通过相位因子折叠回原胞内。如果漏掉这个叠加你算出来的投影能带会在能谱形状上就出错后续一切分析都没有意义。2.4 slab矩阵的Matlab代码骨架用Matlab构造x方向开边界的slab矩阵我建议用稀疏矩阵。核心代码如下function H ssh2d_slab(Nx, ky, v, w) % x方向开边界Nx个原胞y方向周期ky 属于 [0, 2*pi) N 4 * Nx; H sparse(N, N); idx (jx, s) 4 * (jx - 1) s; for jx 1:Nx % 原胞内水平边 H(idx(jx,1), idx(jx,2)) H(idx(jx,1), idx(jx,2)) v; H(idx(jx,3), idx(jx,4)) H(idx(jx,3), idx(jx,4)) v; % 原胞内垂直边同时叠加跨原胞的垂直跳跃 H(idx(jx,1), idx(jx,3)) H(idx(jx,1), idx(jx,3)) v w * exp(1i*ky); H(idx(jx,2), idx(jx,4)) H(idx(jx,2), idx(jx,4)) v w * exp(1i*ky); % x方向跨原胞跳跃开边界不乘相位 if jx Nx H(idx(jx,2), idx(jx1,1)) H(idx(jx,2), idx(jx1,1)) w; H(idx(jx,4), idx(jx1,3)) H(idx(jx,4), idx(jx1,3)) w; end end H H H; end这段代码取各向同性参数v和wsite排列与前面的bond清单一致。跑通后要验证一下固定ky0把Nx增大最低能级附近的态密度应该逐渐趋向体态能带同时slab矩阵的本征值一定是实数。如果出现复数先检查有没有漏掉H H H。3. 原胞能带沿高对称路径扫描与能带特征判读3.1 布里渊区高对称路径的选择二维正方晶格的布里渊区高对称点是Γ(0,0)、X(π,0)、M(π,π)。通常沿着Γ→X→M→Γ扫一圈就能完整看到能带的色散特征。注意这里的坐标是约化后的无量纲动量晶格常数a取1所以布里渊区边界在π而不是π/a。实际扫描时k点采样密度很重要。我一般每段路径取200个点算下来一条能带图大概600个k点。4×4的矩阵对角化实在太快Matlab跑这种规模基本是瞬间出结果所以不用吝啬采样密度。代码可以这样组织v 0.4; w 1.0; npts 200; kpath [linspace(0,0,npts) , linspace(0,pi,npts) , linspace(pi,pi,npts); linspace(0,0,npts) , linspace(0,0,npts) , linspace(0,pi,npts)]; Ebands zeros(4, 3*npts); for i 1:3*npts Hk ssh2d_Hk(kpath(1,i), kpath(2,i), v, w); Ebands(:,i) sort(eig(Hk)); end plot(1:3*npts, Ebands, k-, LineWidth, 1);3.2 能带结果怎么读gap、简并和半金属点参数取v0.4、w1.0时体态能带会打开一个有限gap。原胞能带一共有四条带因为每个原胞有四个site。在Γ点附近两条导带和两条价带分别简并简并的来源是原胞内site1/site4与site2/site3的对称性。沿着X点和M点能带色散会表现出明显的曲率变化。如果你把参数调到vw会看到某个k点处导带底和价带顶刚好碰到系统变成半金属。这个闭合点其实是拓扑相变的临界点信号扫描v从1.5降到0.5的过程中gap会先减小到零再重新打开。读能带图时有两个直觉要建立第一gap打开不代表有边界态它只是边界态存在的必要条件第二能带图的对称性E(k)E(-k)是时间反演对称的体现如果你算出来的能带图左右不对称大概率是哈密顿量构造出了问题。3.3 验证代码正确性的几个自查方法我每次写完紧束缚代码都会做三层验证强烈建议你也养成这个习惯检查厄米性。对随机(site)取一个k点计算norm(Hk - Hk)结果必须是0。slab矩阵同理。检查时间反演对称。把E(kx,ky)和E(-kx,-ky)都算出来两组本征值应当完全一致。检查极限情况。当vw1时模型退化为均匀正方格子的一部分能带会出现特定的半金属点当v0时体系变成一组完全解耦的胞间二聚化链的组合边界态特征应该非常明显且容易识别。这三层验证都通过后才能放心往下做投影能带。4. 投影能带引入表面层投影把隐藏的边界态挖出来4.1 从无限周期到slab为什么投影能带能显示拓扑边界态原胞能带对边界态是“看不见”的因为布洛赫定理默认系统无限大、没有边界。可真实的物理样品一定有边界。投影能带的思路是沿某个方向把体系截断只保留有限个原胞另一个方向仍用周期边界。这样一来哈密顿量就包含了边界信息本征态里也就可能出现局域在边界附近的新态。但光算能谱还不够因为能谱里体态的数量远多于边界态肉眼扫过去很难分辨。所以要再算一个投影权重对每个本征态计算它落在表层原胞上的概率幅平方和然后画图时用点的大小或颜色把这个权重表示出来。体态分布在全体系投影权重很小边界态集中在前几个原胞权重接近1。这样边界态就在图上一眼能看见。4.2 构造slab哈密顿量沿x开边界沿y保持周期在Matlab里就是用第2.4节的ssh2d_slab函数。对每个ky值调用一次这个函数得到维度4Nx×4Nx的矩阵然后对角化得到本征值E_i和本征向量psi_i。这里有一个性能选择的点如果Nx取60矩阵维度是240×240用eig(full(H))直接算完全没问题如果Nx取200以上矩阵变稀疏可以考虑用eigs只算特定能量范围的本征值。但我的经验是投影能带最好把所有本征值都拿到因为画全谱时体态本身也有参考价值所以前期用eig(full(H))最省心。ky扫描范围取0到π即可。因为时间反演对称让能带在负ky区间对称重复不需要扫满整个周期。4.3 投影权重的定义与表面态判据假定取Nx80表层取前3个原胞和后3个原胞。对第i个本征态投影权重P_i定义为P_i sum(|psi_i(表层索引)|.^2)表层索引就是所有属于jx1,2,3以及jxNx-2,Nx-1,Nx的site对应的行号。这样得到的P_i是一个0到1之间的数。对体态来说表层3个原胞占总原胞数的6/800.075权重大致在这个量级边界态如果完全局域在表层几个原胞内权重会接近0.5甚至更高。画图时我建议用散点图把ky、E_i作为横纵坐标点的颜色或大小映射到P_i。一个简单有效的实现如下ky_vec linspace(0, pi, 100); Nx 80; v 0.4; w 1.0; allE []; allKy []; allP []; for ky ky_vec H ssh2d_slab(Nx, ky, v, w); [psi, E] eig(full(H)); E diag(E); % 表层索引 surf_idx []; for jx [1:3, Nx-2:Nx] surf_idx [surf_idx, 4*(jx-1)1:4*(jx-1)4]; end P sum(abs(psi(surf_idx,:)).^2, 1); allE [allE; E]; allKy [allKy; ky * ones(size(E))]; allP [allP; P]; end scatter(allKy, allE, 8 30*sqrt(allP), allP, filled);这里用sqrt的作用是压缩权重分布让中等大小的投影权重也能被肉眼识别。如果数据点过大可以适当调整8和30这两个数值。颜色映射一般选parula或turbo不要用默认的jetjet在高亮端容易让人产生伪边界的感觉。4.4 收敛性检查slab厚度要取多少层才能看到干净的边界态这是个实操里很容易忽略的问题。slab厚度Nx太小两侧边界的边界态会相互耦合导致本征值偏离零能位置甚至在gap里看不到干净的边界能带。我在v0.4、w1.0参数下测试过Nx20时边界态还有点展宽因为左右表面的波函数交叠Nx60时基本收敛投影能带里gap中的亮色分支位置已经稳定Nx120时结果和Nx60几乎一样。所以保守起见日常计算取Nx80到100就够了。如果资源紧张至少不要低于40。另外投影权重对表层厚度的选择也敏感。表层只取1个原胞时边界态权重最大但体态也可能因为端点效应出现伪高权重取3到5个原胞能在边界态和体态之间拉开差距我实际用下来取3个原胞最舒服。5. 二维SSH模型真正的看点角态与高阶拓扑的数值诊断5.1 为什么二维SSH的边界态不在“能带间隙中央”而可能在角上一维SSH链的边界态是两端各一个零能态二维SSH的情况要微妙很多。普通的二维拓扑绝缘体比如Chern绝缘体或量子自旋霍尔绝缘体对应的是沿一维边缘传播的手性边缘态但二维SSH模型在这种四格点排列下如果v/w小于阈值体系并不出现完整的一维边缘能带而是在四个角上出现零能的角态。这类物态统称高阶拓扑绝缘体体态在d维gapped表面态在d-1维而拓扑保护的束缚态进一步降维到d-2维。对于二维体系高阶拓扑态就是0维的角态。这在投影能带图里会带来一个直观后果你不一定能看到像一维SSH那样贯穿gap的清晰边界色散分支。更常见的现象是在ky0或者kyπ附近出现一些能量靠近零的离散点或平带投影权重集中在表层。第一次看到这种结果的人容易怀疑代码写错了以为“gap里没有连续边界带就说明没有拓扑”。这是二维SSH最容易误判的地方。正确的诊断方式是把x和y方向都开边界直接算有限尺寸的晶块看零能附近有没有四个角态。5.2 有限尺寸全开边界的对角化直接观察角态实空间分布要验证二维SSH拓扑相最直接的方法是构造一个Nx×Ny个原胞的完全开边界矩阵。这个矩阵的构造比slab麻烦一点但原理一样把x和y方向的跨原胞hopping都显式加在矩阵里不引入任何周期方向的相位因子。核心代码骨架如下function H ssh2d_finite(Nx, Ny, v, w) N 4 * Nx * Ny; H sparse(N, N); id (jx, jy, s) 4 * ((jy-1)*Nx (jx-1)) s; for jy 1:Ny for jx 1:Nx % 胞内边 H(id(jx,jy,1), id(jx,jy,2)) H(id(jx,jy,1), id(jx,jy,2)) v; H(id(jx,jy,3), id(jx,jy,4)) H(id(jx,jy,3), id(jx,jy,4)) v; H(id(jx,jy,1), id(jx,jy,3)) H(id(jx,jy,1), id(jx,jy,3)) v; H(id(jx,jy,2), id(jx,jy,4)) H(id(jx,jy,2), id(jx,jy,4)) v; % x方向胞间 if jx Nx H(id(jx,jy,2), id(jx1,jy,1)) H(id(jx,jy,2), id(jx1,jy,1)) w; H(id(jx,jy,4), id(jx1,jy,3)) H(id(jx,jy,4), id(jx1,jy,3)) w; end % y方向胞间 if jy Ny H(id(jx,jy,3), id(jx,jy1,1)) H(id(jx,jy,3), id(jx,jy1,1)) w; H(id(jx,jy,4), id(jx,jy1,2)) H(id(jx,jy,4), id(jx,jy1,2)) w; end end end H H H; end对得到的哈密顿量求零能附近的前若干个本征态然后画实空间概率分布。我通常取NxNy30用eigs(H, 8, smallestreal)直接提取能量最靠近零的8个态。在v0.4、w1.0参数下你会看到4个简并的零能态概率密度分别集中在四个角上。把四个角的格点坐标标出来每个角态在对应的角附近有一个清晰峰。如果参数换成v1.2、w1.0零能附近的态不再局域在角上体系回归平庸绝缘体。这一步做完二维SSH的高阶拓扑性质才算真正闭环。5.3 拓展验证手段Wilson loop的思路简述角态计算是直观证据但如果你想把这套代码延伸到更复杂的模型比如加上次近邻hopping或者无序Wilson loop是更普适的判断工具。思路不复杂对每个固定的kx先沿ky积分占据态的Berry联络得到一条一维Wilson loop然后对Wilson loop矩阵取本征值得到Wannier center的位置分布。二维SSH拓扑相的特征是Wannier center流呈现特定的绕数。由于这一步不涉及投影能带的绘制这里不展开代码但如果你做完角态验证后还想再严谨一点可以从非阿贝尔Wilson loop入手这是论文里最常用的数值判据。6. 实操中的坑与调试经验让能带图真正成为论文级的图6.1 投影权重的可视化别让体态的光芒盖住边界态投影能带图最常见的问题是整张图一团黑体态的多条能带挤在一起边界态的亮色点被淹没。我的解决办法是区分两类点先用浅灰色画出所有态再用亮度较高的颜色覆盖投影权重超过0.2的态。这样即使边界态数量少也不会被体态掩盖。具体实现上可以把P_i小于阈值的点统一设置为一个很淡的颜色大于阈值的点单独画一层。还有个小技巧在scatter里用sqrt(P)作为颜色数据时要指定colormap为parula颜色条的最小值不要从0开始而是从P的最小非零值附近开始否则大部分点都落在同一色阶上。6.2 关于本征向量相位和简并态的注意事项投影权重之所有能这样简单地算是因为它只涉及概率幅的平方不涉及本征向量的U(1)相位。但如果你进一步算Wilson loop或者偶极矩就必须注意本征值的简并和本征向量的规范选择。二维SSH模型在Γ点附近有简并直接数值对角化得到的简并本征向量可以任意旋转如果程序里假设了某个固定相位顺序后续拓扑量计算就可能出错。这属于进阶问题现阶段只要记住投影能带计算对简并态是安全的因为投影权重对简并子空间内部的么正变换不变一旦开始算Wilson loop就要用光滑规范或者平行输运方法。6.3 参数扫描与相图绘制从单条能带到完整相图单条能带算通之后最自然的扩展是扫参数。把w固定为1让v在0.2到1.8之间变化对每个v值都算一次零能附近的角态能级就能画出“体态gap闭合—打开”的相图临界线。更精细的做法是同时算封边界下的角态存在条件在(v, w)二维参数空间里用颜色标记最低能级是否为零。用Matlab的contourf画出来会看到拓扑相区和平庸相区之间有一条清晰的边界线。这个过程几乎是全自动的只需要把上面几段代码嵌套到一个循环里唯一要注意的是每次对角化后要对本征值排序否则相图里会出现噪声般的跳点。6.4 常见问题快速排查表问题可能原因解决方式本征值出现虚部矩阵不是厄米矩阵检查是否补了H H检查相位符号是否统一能带图左右不对称动量路径写错或矩阵元符号不满足时间反演显示kx和ky的取值范围检查路径数据投影能带边界态不明显slab太薄或投影表层选择太少/太多Nx至少取60表层取3个原胞附近零能附近出现大量体态参数处于半金属点附近把v/w远离1比如取0.4/1.0完全开边界时角态只有两个晶块尺寸太小四角之间耦合未消失增大Nx和Ny至少30×30以上相图扫描线不光滑简并能级排序不稳定先对每个k点的本征值排序再做后续统计用Matlab做这类紧束缚计算最大的优势不是速度而是矩阵操作的直观性。从4×4的H(k)到几百乘几百的slab矩阵再到几千乘几千的有限晶块矩阵语法几乎没有变化同一套site索引逻辑可以一路复用。这也是我一直建议入门拓扑计算的人先用Matlab把这个模型跑通的原因模型足够有代表性代码量又被压缩得很少等把物理和数值方法都搞清楚之后再迁移到Python或者其他语言也不会太困难。最后分享一个小习惯每次调整参数后先把vw1的对照结果跑一遍。这个参数下体系能带图必须出现gapless的特征任何偏离都说明代码骨架被改坏了。先验基线保住再谈相图和拓扑能省去大量调试时间。