二维SSH模型Matlab复现：紧束缚能带、投影能带与高阶拓扑角态

动手做二维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_x；site3向上与上方原胞的site1耦合，site4向上与上方原胞的site2耦合，振幅为w_y。

为了减少参数让结果好看，我建议前期先取各向同性：v_x = v_y = v，w_x = w_y = w。这样模型只保留两个参数，v和w的比值决定相图。v<w对应拓扑非平庸的区间，v>w是平庸区间，两者中间发生能带gap闭合再打开。这个模型的漂亮之处在于它不仅能展示原胞能带和投影能带的差别，还能顺带展示二维体系里比一维推广更有意思的现象：不止有边缘态，还可能出现四角束缚态，也就是所谓的高阶拓扑角态。第5章会专门说这个。

1.3 先给结论：原胞能带和投影能带分别看什么

原胞能带，或者说bulk能带，是把二维晶格当作无限周期系统，玻尔兹曼变换后得到H(kx,ky)，在布里渊区高对称线上扫描本征值。它回答的问题是：体态是金属还是绝缘体，gap在哪里，色散关系长什么样。投影能带则完全不同，它回答的问题是：如果在这个方向切出边界，gap里会不会出现边界态。做法是沿着某一个方向保留有限个原胞，另一个方向继续保持周期边界，然后计算能谱，再按本征向量在表面层的振幅加权，让局域在边界上的态在图上“变亮”。这就是标题里两个能带的实质分工。接下来我从紧束缚模型出发，先把代码层面的事说清楚。

2. 紧束缚哈密顿量的Matlab构造：从实空间参数到H(k)

2.1 先列bond清单：把模型约定写清楚

写代码之前，我习惯先把所有hopping列成一张表，避免原胞内外一堆耦合混淆。本文约定的bond清单如下：

原胞内的四个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(jx+1,1)) = H(idx(jx,2), idx(jx+1,1)) + w; H(idx(jx,4), idx(jx+1,3)) = H(idx(jx,4), idx(jx+1,3)) + w; end end H = H + H'; end

这段代码取各向同性参数v和w，site排列与前面的bond清单一致。跑通后要验证一下：固定ky=0，把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、简并和半金属点

参数取v=0.4、w=1.0时，体态能带会打开一个有限gap。原胞能带一共有四条带，因为每个原胞有四个site。在Γ点附近，两条导带和两条价带分别简并，简并的来源是原胞内site1/site4与site2/site3的对称性。沿着X点和M点，能带色散会表现出明显的曲率变化。如果你把参数调到v=w，会看到某个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)都算出来，两组本征值应当完全一致。
• 检查极限情况。当v=w=1时，模型退化为均匀正方格子的一部分，能带会出现特定的半金属点；当v=0时，体系变成一组完全解耦的胞间二聚化链的组合，边界态特征应该非常明显且容易识别。

这三层验证都通过后，才能放心往下做投影能带。

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 投影权重的定义与表面态判据

假定取Nx=80，表层取前3个原胞和后3个原胞。对第i个本征态，投影权重P_i定义为：

P_i = sum(|psi_i(表层索引)|.^2)

表层索引就是所有属于jx=1,2,3以及jx=Nx-2,Nx-1,Nx的site对应的行号。这样得到的P_i是一个0到1之间的数。对体态来说，表层3个原胞占总原胞数的6/80=0.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，不要用默认的jet，jet在高亮端容易让人产生伪边界的感觉。

4.4 收敛性检查：slab厚度要取多少层才能看到干净的边界态

这是个实操里很容易忽略的问题。slab厚度Nx太小，两侧边界的边界态会相互耦合，导致本征值偏离零能位置，甚至在gap里看不到干净的边界能带。我在v=0.4、w=1.0参数下测试过：Nx=20时，边界态还有点展宽，因为左右表面的波函数交叠；Nx=60时基本收敛，投影能带里gap中的亮色分支位置已经稳定；Nx=120时结果和Nx=60几乎一样。所以保守起见，日常计算取Nx=80到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的清晰边界色散分支。更常见的现象是，在ky=0或者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(jx+1,jy,1)) = H(id(jx,jy,2), id(jx+1,jy,1)) + w; H(id(jx,jy,4), id(jx+1,jy,3)) = H(id(jx,jy,4), id(jx+1,jy,3)) + w; end % y方向胞间 if jy < Ny H(id(jx,jy,3), id(jx,jy+1,1)) = H(id(jx,jy,3), id(jx,jy+1,1)) + w; H(id(jx,jy,4), id(jx,jy+1,2)) = H(id(jx,jy,4), id(jx,jy+1,2)) + w; end end end H = H + H'; end

对得到的哈密顿量求零能附近的前若干个本征态，然后画实空间概率分布。我通常取Nx=Ny=30，用eigs(H, 8, 'smallestreal')直接提取能量最靠近零的8个态。在v=0.4、w=1.0参数下，你会看到4个简并的零能态，概率密度分别集中在四个角上。把四个角的格点坐标标出来，每个角态在对应的角附近有一个清晰峰。如果参数换成v=1.2、w=1.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 常见问题快速排查表

用Matlab做这类紧束缚计算，最大的优势不是速度，而是矩阵操作的直观性。从4×4的H(k)到几百乘几百的slab矩阵，再到几千乘几千的有限晶块矩阵，语法几乎没有变化，同一套site索引逻辑可以一路复用。这也是我一直建议入门拓扑计算的人先用Matlab把这个模型跑通的原因：模型足够有代表性，代码量又被压缩得很少，等把物理和数值方法都搞清楚之后，再迁移到Python或者其他语言也不会太困难。

最后分享一个小习惯：每次调整参数后，先把v=w=1的对照结果跑一遍。这个参数下体系能带图必须出现gapless的特征，任何偏离都说明代码骨架被改坏了。先验基线保住，再谈相图和拓扑，能省去大量调试时间。