二维SSH模型是研究拓扑相变的经典玩具模型。今天咱们用Matlab实现几个有意思的计算，手把手看看电子在蜂窝状晶格里的骚操作

基于紧束缚模型matlab计算二维SSH模型 包含 1.投影能带 2.原胞能带 3.x方向有限结构y方向周期的能带 x方向y方向均自由边界条件的体态场分布 5.x方向y方向均自由边界条件的边界态态场分布 6.x方向y方向均自由边界条件的本征值分布先上硬菜——构造紧束缚哈密顿量。二维SSH模型每个原胞包含四个原子跃迁强度在x/y方向交替变化function H SSH2D_Hamiltonian(kx, ky, vx, wx, vy, wy) % 基矢顺序[A,B,C,D] H zeros(4,4); % x方向跃迁 H(1,2) vx wx*exp(-1i*kx); H(2,1) conj(H(1,2)); H(3,4) vx wx*exp(1i*kx); H(4,3) conj(H(3,4)); % y方向跃迁 H(1,3) vy wy*exp(-1i*ky); H(3,1) conj(H(1,3)); H(2,4) vy wy*exp(1i*ky); H(4,2) conj(H(2,4)); end这里kx,ky是动量空间坐标v和w对应不同方向的跃迁强度。注意复指数项的正负号这关系到相邻原胞的相位积累。投影能带计算沿着高对称点画能带就像把三维地形投影到二维地图上kx linspace(-pi, pi, 50); ky linspace(-pi, pi, 50); [KX, KY] meshgrid(kx, ky); E zeros(4, 50, 50); for i 1:50 for j 1:50 H SSH2D_Hamiltonian(KX(i,j), KY(i,j), 0.5, 1.2, 0.7, 0.9); E(:,i,j) eig(H); end end figure; for n 1:4 surf(squeeze(E(n,:,:)), EdgeColor,none); hold on; end这里能看到四个能带的立体分布。注意当v基于紧束缚模型matlab计算二维SSH模型 包含 1.投影能带 2.原胞能带 3.x方向有限结构y方向周期的能带 x方向y方向均自由边界条件的体态场分布 5.x方向y方向均自由边界条件的边界态态场分布 6.x方向y方向均自由边界条件的本征值分布有限尺寸效应当x方向截断成有限链时边界态就藏不住了Nx 20; % x方向原胞数 Ny 1; % y方向保持周期 % 构造块对角哈密顿量 H zeros(4*Nx, 4*Nx); for ix 1:Nx % 原胞内耦合 block SSH2D_Hamiltonian(0, 0, 0.5, 1.2, 0.7, 0.9); H(4*(ix-1)1:4*ix, 4*(ix-1)1:4*ix) block; % 原胞间耦合 if ix Nx H(4*(ix-1)2, 4*ix1) 1.2; % x方向跃迁 H(4*ix1, 4*(ix-1)2) 1.2; end end [V,D] eig(H); edge_state V(:,2*Nx); % 取中间能隙附近态 figure; plot(abs(edge_state(1:4:end)).^2, o-); % A原子占据概率 hold on; plot(abs(edge_state(4:4:end)).^2, s-); % D原子占据概率这个画出来的态密度在两端明显增强妥妥的边界态。有趣的是A和D原子的分布不对称说明手性对称性破缺。体态与边界态同台竞技当x,y都取自由边界时得用全尺寸对角化Nx 10; Ny 10; H zeros(4*Nx*Ny, 4*Nx*Ny); for ix 1:Nx for iy 1:Ny idx (ix-1)*Ny iy; % 原胞内耦合 H(4*(idx-1)1:4*idx, 4*(idx-1)1:4*idx) ... SSH2D_Hamiltonian(0, 0, 0.5, 1.2, 0.7, 0.9); % x方向邻接 if ix Nx jdx idx Ny; H(4*(idx-1)2, 4*(jdx-1)1) 1.2; H(4*(jdx-1)1, 4*(idx-1)2) 1.2; end % y方向邻接 if iy Ny jdx idx 1; H(4*(idx-1)4, 4*(jdx-1)3) 0.9; H(4*(jdx-1)3, 4*(idx-1)4) 0.9; end end end [V,D] eig(H); eigenvalues diag(D); figure; scatter(real(eigenvalues), imag(eigenvalues), 10, filled);这里本征值分布会呈现明显的能隙而靠近零能的态对应边缘态。用scatter画出的本征值图里如果看到能隙中有孤立点那就是拓扑保护的边界态在招手。最后来个炫技操作——画三维态分布[IX, IY] meshgrid(1:Nx, 1:Ny); wave reshape(abs(V(:,200)).^2, 4, Nx, Ny); % 取第200个本征态 figure; subplot(221); scatter3(IX(:), IY(:), squeeze(wave(1,:,:)), 50, squeeze(wave(1,:,:)), filled); title(A原子分布); % 类似画出B、C、D原子分布...这种三维散点图能清晰展示波函数在空间中的局域情况。当看到波函数像镶边一样分布在样品边缘时恭喜你成功抓住了拓扑边界态的小尾巴