简介:SSH模型是固体物理中描述一维链状体系电子与声子行为的重要理论模型,也是计算物理课程中理解能带形成的典型范例。该压缩包提供一个MATLAB求解脚本,用于构建一维弹簧质点链并数值求解其能带结构,特别适合固体物理、计算物理方向的本科生和科研人员,既可用于课堂演示,也可作为快速仿真工具。资源以RAR格式打包,共1个文件,即一个.m脚本,整体大小仅674B,代码精炼、依赖少,易于阅读和修改参数。目前已有139人学习浏览。运行脚本后,用户可以直观观察耦合强度变化如何影响能带形状,理清电子能带与声子能带的对应关系,并通过调整模型参数进一步探索拓扑相变特性,为学习拓扑绝缘体、声子晶体等前沿专题奠定数值模拟基础。 这段时间一直在整理计算物理相关的教学项目,翻到一个很有意思的案例:压缩包名叫“SSH_Model_for_1D_Spring_Mass_Chain”,里面是一套用一维弹簧质量链来演示SSH模型(Su-Schrieffer-Heeger)的代码和笔记。乍一看这个名字挺绕,但拆开之后你会发现,它横跨了固体物理、能带理论和计算物理三个核心板块,而且用的是非常直观的经典力学系统。对于正在学凝聚态物理或者刚入坑拓扑物态的人来说,这是一个特别好的切入点。这篇内容我就围绕这个项目,把背后的物理模型、能带计算思路、数值实现过程以及我调试时踩过的坑一并整理出来。
1. 这个项目到底在解决什么问题
1.1 标题信息拆解:SSH模型、一维链、能带与声子
先把这个压缩包名字拆开看。“SSH_Model”指的是Su-Schrieffer-Heeger模型,这是1979年三个人提出的一维紧束缚模型,最初用来描述聚乙炔(polyacetylene)的金属-绝缘体转变和孤子激发。“1D_Spring_Mass_Chain”代表一维弹簧质量链,也就是一堆质量块由弹簧串联的经典体系。两者放在一起,本质是在做一件事——用弹簧的松紧交替来模拟SSH模型里电子跳跃积分的大小交替。
后面跟着的“能带”和“能带 声”这两个关键词,更是直接点出了项目的核心输出:计算这个一维链在周期性结构下的色散关系,也就是能带;而“声”对应的其实是声子色散,{\displaystyle \omega }{\displaystyle k}关系。因为在经典弹簧链里,格波色散就是声子的能带,和电子在周期性势场中的能带结构在数学上同构。算完声子色散,再去对比电子SSH模型的能带,两者给出的物理图像完全一致。
“计算物理”则是整个项目的方法论标签:不是手推公式,而是借助数值对角化、周期边界条件、开放边界条件等手段,把色散关系、本征态、边界态这些抽象概念用可视化的方式算出来、画出来。
1.2 为什么用经典弹簧链来学拓扑能带
很多人在初学紧束缚模型时,都会被“跃迁积分”“胞内/胞间耦合”这些术语劝退。可一旦换成弹簧质量链,一切突然变得特别好理解:把质量块当成原子,把弹簧当成电子跃迁或者原子间相互作用,弹簧软硬交替就对应SSH模型里的交变跃迁强度。
这种映射的好处有三层:
- 一维链的动力学方程直接就是牛顿第二定律,不需要先啃二次量子化和 Bloch 定理,门槛低。
- 色散关系通过求解动力学矩阵的本征值就能得到,不需要求复数波函数,用 MATLAB、Python 都能很轻松地实现。
- 边界态(零能模)在弹簧链里可以直接看到——如果链两端连接的是软弹簧,就会出现一种只在边界附近振荡、频率为零的局域模,对应电子体系里的拓扑边界态,这种“看得见摸得着”的直观性非常难得。
所以在我看来,这个项目不只是一段代码,更是一套兼顾了理论理解与数值实操的完整学习方案。适合固体物理课程配套实验、计算物理作业、或者想搞懂拓扑绝缘体入门概念的自学者。
2. 核心物理:弹簧质量链如何长出能带和拓扑相
2.1 从运动方程到色散关系的基础推导
我们先从最简单的情况说起:原子质量为m,弹簧劲度系数为k,相邻原子间距为a。对第n个原子写牛顿方程:
[ m\ddot{u}n = k(u{n+1} + u_{n-1} - 2u_n) ]
采用平面波试探解 (u_n = A e^{i(qna - \omega t)}),代入后得到:
[ -m\omega^2 = k(e^{iqa} + e^{-iqa} - 2) = 2k(\cos(qa)-1) ]
所以色散关系为:
[ \omega(q) = 2\sqrt{\frac{k}{m}}\left|\sin\left(\frac{qa}{2}\right)\right| ]
这就是单原子链的声学支色散。在长波极限 (q\to 0) 下,(\omega \approx \sqrt{k/m}, qa),对应的声速是 (v = a\sqrt{k/m}),和连续介质中的弹性波完全一致。
这一步是整个项目的物理起点。它告诉我们:色散关系本质上是振动本征模式的频率与波矢之间的关系,而这个关系完全由动力学矩阵决定。接下来只要把弹簧分成“强-弱-强-弱”交替排列,立刻就会用到双原子链的结论。
2.2 二聚化与能带折叠
SSH模型的关键是“二聚化”。假设每个原胞内有两个原子A、B,原子间有两种弹簧:胞内弹簧 (k_1),胞间弹簧 (k_2),并且 (k_1 \neq k_2)。此时原胞长度为 (2a),但实际原子间距还是a。我们可以写出两个耦合的运动方程:
[ m\ddot{u}{n,A} = -k_1(u{n,A}-u_{n,B}) - k_2(u_{n,A} - u_{n-1,B}) ] [ m\ddot{u}{n,B} = -k_1(u{n,B}-u_{n,A}) - k_2(u_{n,B} - u_{n+1,A}) ]
设解为Bloch形式 (u_{n,A}=A e^{iq(2na)}),(u_{n,B}=B e^{iq(2na+a)}),代入后得到2x2动力学矩阵:
[ D(q) = \frac{1}{m}\begin{pmatrix} k_1+k_2 & -k_1 - k_2 e^{-i2qa} \ -k_1 - k_2 e^{i2qa} & k_1+k_2 \end{pmatrix} ]
对角化后会出现两支色散:声学支和光学支。在布里渊区边界 (q=\pi/(2a)) 处,两支之间的能隙大小为 (|k_1-k_2|)。如果 (k_1=k_2),能隙闭合,系统回到单原子链的折叠能带——就像把一段直线按两倍周期折叠回来,两条能带在边界相接。
这里可以看到能带理论的第一个深刻结论:周期性势场会让能谱在倒空间里折叠,打开的能隙大小由交替势的强度决定。对电子SSH模型来说,这个交替势就是交变的跳跃积分 (t_1) 和 (t_2),而对弹簧链来说,则是 (k_1) 与 (k_2)。
2.3 SSH模型的紧束缚哈密顿量对应关系
在电子体系里,SSH模型的动量空间哈密顿量可以写成:
[ H(q) = \begin{pmatrix} 0 & t_1 + t_2 e^{-iqa} \ t_1 + t_2 e^{iqa} & 0 \end{pmatrix} ]
这个哈密顿量的本征值给出两条电子能带:
[ E_\pm(q) = \pm \sqrt{t_1^2 + t_2^2 + 2t_1 t_2 \cos(qa)} ]
对比弹簧链的动力学矩阵,你会发现数学结构几乎完全一样,只是矩阵元里的弹簧劲度系数 (k_1, k_2) 换成了跃迁积分 (t_1, t_2)。这就是两个看似风马牛不相及的体系之间存在“同构关系”的绝佳例子。做完弹簧链的数值模拟,再去推导电子能带,很多原本抽象的东西一下子就有了画面。
2.4 拓扑相变判据:能隙闭合与边界态
一个被反复讨论的问题是:对SSH模型,什么时候系统是拓扑非平凡的?关键指标有两个。
第一个是体态能隙是否闭合。当 (k_1 = k_2)(或者说 (t_1=t_2))时,能带在布里渊区边界闭合,系统发生拓扑相变。在 (k_1 < k_2) 和 (k_1 > k_2) 两侧,能带结构看起来非常相似——都有能隙,但拓扑性质完全不同。这就引出了第二个指标:缠绕数(winding number)或者Zak相位。弹簧链里没法直接测电子波函数的相位,但通过观察开放边界条件下的局域态,就能非常直观地判断拓扑性:当 (k_2 > k_1) 时,链两端出现零能边界模;而当 (k_1 > k_2) 时,没有边界态。这个差异在电子体系里对应的是“有/无拓扑保护端点态”。弹簧链的边界模可以非常直观地通过振动振幅分布看到——端点振动幅度特别大,中间呈指数衰减,这就是实际物理中的“边缘态”。
3. 数值实现:从哈密顿量到能带图与边界态
3.1 参数设定与原胞划分
我自己用Python复现了一遍这个项目,参数选得很常规:原子数 (N=200),原胞数 (N_c=100),每个原胞两个原子。弹簧劲度系数设为 (k_1=1.0)、(k_2=3.0),质量全部取 (m=1.0)。这个比例对应电子SSH模型里 (t_2>t_1) 的拓扑非平凡区。实际调试时可以先用这个参数组合,因为它给出的能隙适中,边界态也足够局域,图像特别清楚。
需要提前说明的是:直接使用开放边界条件时,动力学矩阵是 ((2N_c)\times(2N_c)) 的实对称矩阵,数值对角化非常快,100个原胞也就零点几秒。但如果你想要解析色散关系曲线,还是应该先对无限周期体系取Bloch波矢,再对每个q点对角化2x2矩阵。
3.2 周期性边界条件下的能带计算
周期性边界条件下的核心就是扫描第一布里渊区中的波矢q。因为原胞长度为 (2a),布里渊区范围是 (q \in [-\pi/(2a), \pi/(2a)])。代码思路大致如下:
import numpy as np import matplotlib.pyplot as plt m = 1.0 k1, k2 = 1.0, 3.0 Nq = 400 qs = np.linspace(-np.pi/2, np.pi/2, Nq) bands = [] for q in qs: D = np.zeros((2,2), dtype=complex) D[0,0] = (k1 + k2) / m D[1,1] = (k1 + k2) / m D[0,1] = -(k1 + k2 * np.exp(-2j*q)) / m D[1,0] = -(k1 + k2 * np.exp(2j*q)) / m evals = np.linalg.eigvalsh(D) bands.append(np.sort(evals)) bands = np.array(bands) plt.plot(qs, bands[:,0], 'b') plt.plot(qs, bands[:,1], 'r') plt.xlabel('q') plt.ylabel(r'$\omega^2$') plt.show()注意这里我画的是 (\omega^2) 而不是 (\omega)。原因是动力学矩阵的本征值直接给出的是频率平方,取平方根后画 (\omega(q)) 也可以,但有些情况下 (\omega^2) 曲线更平滑,能隙打开的位置更清楚。我建议两个都画一下,对比会有收获。
真实输出结果中,两个支在 (q=0) 处时,声学支频率为0,光学支频率为 (\sqrt{2(k_1+k_2)/m});在 (q=\pi/2) 处,两支之间的频率间隔正好是 (2|k_1-k_2|/m) 的平方根相关量。这个“间隙随 (|k_1-k_2|) 线性增大”的规律,在图上非常直观。
3.3 开放边界条件与拓扑边界态判定
周期性边界条件只能看到体态能带,边界态必须在开放边界条件下才会出现。构建实空间动力学矩阵时,每个原子与最近邻之间的弹簧要么是 (k_1) 要么是 (k_2),取决于它处于原胞内还是原胞间。我用一个循环来组装 (2N_c\times 2N_c) 三对角矩阵,然后调用 numpy 的 eigh 直接对角化。
Nc = 100 Ntot = 2 * Nc M = np.zeros((Ntot, Ntot)) # 交替填充 k1 (胞内) 和 k2 (胞间) for i in range(Ntot - 1): if i % 2 == 0: M[i, i+1] = M[i+1, i] = -k1 / m else: M[i, i+1] = M[i+1, i] = -k2 / m # 对角项 for i in range(Ntot): M[i, i] = (k1 + k2) / m if i == 0 or i == Ntot - 1: M[i, i] = (k1 + k2) / m # 两端各少一个弹簧,更精确处理见下文 evals, evecs = np.linalg.eigh(M)这里有一个容易踩坑的地方:链两端各只有一个弹簧,所以端点原子的对角项应当是 (k_1/m) 或 (k_2/m),而不是 ((k_1+k_2)/m)。你如果偷懒不对角项修正,边界态的零频特征会被一个小的虚偏移污染,图上看不明显,但数值上会变得不严格为零。
正确做法是用“交错数列”生成弹簧类型后,逐点统计每个原子求和的弹簧劲度系数之和,再填入对角元。这样得到的本征谱里会清楚地看到:除了两条体态能带,还会在零频率(对应零能)附近出现一个额外的态。画出这个本征向量,空间分布会在两端同时出现两个峰值,中间指数衰减。这,就是拓扑边界态。
3.4 边界态空间分布分析与可视化
判定边界态最直接的方法是画本征向量分布。找本征值最接近0的那个本征向量,然后对它所有分量取模平方:
idx = np.argmin(np.abs(evals)) psi = evecs[:, idx] plt.plot(np.abs(psi)**2) plt.xlabel('site index') plt.ylabel(r'$|\psi|^2$')如果用的是 (k_1=1.0, k_2=3.0),会得到两个峰值分别落在第0号和第N-1号原胞位置,中间距离相隔100个原胞,每个峰值内部按 (|\psi_n|^2 \sim e^{-n/\xi}) 衰减。局域长度大致由 (\xi \approx 1/\ln(k_2/k_1)) 给出。对这个参数,(\ln 3 \approx 1.098),所以 (\xi \approx 0.91) 个原胞,边界模衰减非常快。
设置 (k_1=3.0, k_2=1.0) 之后重新对角化,靠近零频的边界态消失,所有本征态都是体态。这就构成了一组漂亮的对比图:同样的能隙、同样的体态,但拓扑性质完全不同。很多教材用“Zak相位”来区分这两种情况,但在数值实验里,边界态的存在与否就是最直观的拓扑判据。
4. 常见问题与调试心得
4.1 能带图出现“假交叉”或能隙不明显
这是最常碰到的问题。原因多半是波矢扫描范围没取对。SSH模型原胞长度是两倍原子间距,所以第一布里渊区边界在 (q=\pi/(2a)),而不是 (\pi/a)。如果你按单原子链习惯扫描到 (\pi/a),能看到四支能带镜像重复,好像能带发生了交叉,其实只是重复了。另一个原因是没有把 (q=0) 到边界、再回到负边界的完整路径扫完,建议一次性扫完整个BZ再画图。
4.2 频率出现负值
这里要特别提醒:动力学矩阵本征值对应 (\omega^2),确定是非负数。如果出现负本征值,说明组装矩阵时对角项写错了,或者是把相互作用势的符号搞反了。我曾把交叉项写成正号,结果色散关系在边界附近出现虚频,图上直接报Warning。检查方法:对一个均匀链,解析解是 (\omega^2 = (2k/m)(1-\cos q)),用数值代码去对均匀链跑一遍,和解析解对比就能定位问题。
4.3 边界态太“胖”或衰减太慢
如果边界态波及到链中段,说明边界态局域长度超过了链长的一半。解决办法很简单——把原子数加到大几百甚至上千,(\xi) 不变,链变长之后边界态自然就分离开。还有一种情况是你用了周期性边界条件来观察边界态,那肯定看不到边界,因为边界根本不存在。所以计算边界态时一定要用开放边界条件。
下面我整理了一份快速排查表,按实际经验排了优先级:
| 现象 | 可能原因 | 处理方式 |
|---|---|---|
| 能带在边界闭合/交叉 | 波矢范围扫到 ( \pi/a) 而不是 ( \pi/(2a)) | 改成 ( q \in [-\pi/(2a), \pi/(2a)] ) |
| 出现负频率 | 动力学矩阵对角项或交叉项符号错误 | 均匀链解析解对比自检 |
| 看不到边界态 | 用了周期性边界条件 | 切换为开放边界条件 |
| 边界态不局域 | 链长太短或 ( k_2/k_1 ) 接近1 | 增大 N 或调大参数比 |
| 能隙关闭位置不对 | 原胞边界定义反了 | 检查 ( k_1, k_2 ) 顺序与 q 相位因子 |
4.4 关于“能带声”这个标注意义
标题里的“能带 声”我理解为“声子能带”,这在项目里占了很重要的比重。弹簧质量链算出来的色散关系本质上是声子色散,只是通过参数重新映射后能和电子SSH模型的电子能带完全对上。做固体物理课程设计时,如果能同时给出电子能带和声子色散两张图,对理解电子-声子相互作用也很有帮助——链振动的量子化对应声子,电子在链上的跳跃对应电子能带,两者在同一个几何结构上耦合,这正是凝聚态物理最迷人的地方。
我之前试过在这个模型的基础上把每个格点质量也做二聚化,也就是 (m_A\neq m_B),色散关系会变得更复杂,打开声学支与光学支之间的带隙。这个改造方向很有价值,它对应的是更一般的双原子链问题,也是离子晶体中光学声子产生介电响应的物理基础。如果你对声子晶体或者超材料感兴趣,可以把弹簧链模型延伸到二维和三维,用有限元软件模拟出完整的声子带隙结构,那个视觉冲击力更强。
5. 从模型到真实物理:这个项目的延伸价值
5.1 拓扑不变量与体边对应原理
前面提到边界态的出现对应 (k_2>k_1),但如果没有拓扑不变量这个概念,你很难解释为什么边界态这么“结实”。在实际的SSH模型里,系统绕参数空间一圈的Zak相位为 (\pi) 时对应拓扑非平凡相,边界态的存在受拓扑保护,一定出现在体能隙中的零能位置。弹簧链不能直接测电子相位,但体边对应原理的数值体验是完全一致的——不改变体能隙结构,只改变连接方式(比如把 (k_1) 和 (k_2) 的排列顺序对调),边界态就会消失或出现,这正是拓扑相变的特征。
为了更清楚地展示这个性质,建议你额外做一个实验:在链的中间插入一个“缺陷”——把中间某两个原胞之间的弹簧从 (k_2) 换成 (k_1)。此时系统相当于“体态绝缘,但内部出现一个界面”,在缺陷处同样会产生零能局域态,这就是域壁(domain wall)上的边界模。图一出来,你就明白为什么SSH模型能用来解释聚乙炔中的孤子激发。
5.2 从弹簧链到真实拓扑材料
表面上看,一维弹簧和拓扑绝缘体八竿子打不着,但其实很多拓扑物态的核心就在于“交变耦合”这个简单机制。比如声子晶体里的拓扑边界态、电路网络里的SSH模型、光子晶体里的拓扑保护边缘态,全都在沿用这个物理思想。作为计算物理的学习项目,从弹簧链出发去理解拓扑不变量和边界态的对应关系,之后再切换到实际材料的紧束缚模型或者第一性原理计算软件,思路就会顺畅得多。
做计算物理还有一个很实在的好处:你会逐渐意识到,物理直觉和数值验证是互相成就的。模型是先于计算去预测,计算再反向修正直觉。SSH这个项目很适合充当这种“训练场”,因为它的解析解是已知的,任何数值错误都能被立即察觉,特别适合用来建立学习信心和调试能力。
6. 写在最后的实操感悟
如果让我给后人一个建议,那就是:不要只跑参数出图就结束,而是亲手把链的所有本征矢都画出来,逐个去看它们长什么样。低频声学支的本征矢在实空间是一段平滑的正弦波,高频光学支则是相邻原子反向运动,边界态的分布更是直接揭示了拓扑保护的意义——这种把数值结果和物理图像联系起来的观察,才是计算物理最有价值的收获。
我用这个项目带过好几位低年级的同学,他们都能在半天内跑通代码,但真正花时间仔细看本征矢形状、改变参数观察相变的人,后续学拓扑能带理论时明显扎实很多。补一个我在玩这个模型时发现的小技巧:把边界原子的质量逐步减少到接近零,边界态的零频会往复数方向漂移,这在实验上对应的是端头质量负载对边界态频率的调控。如果你以后做声子晶体实验或者超材料设计,这个技巧在设计边界态传感器或波导时非常有用。
本文还有配套的精品资源,点击获取