基于紧束缚模型，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*ix+1) = 1.2; % x方向跃迁
        H(4*ix+1, 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原子分布...

这种三维散点图能清晰展示波函数在空间中的局域情况。当看到波函数像镶边一样分布在样品边缘时，恭喜你成功抓住了拓扑边界态的小尾巴！