3. 一维双原子链与能隙计算

在紧束缚框架下，能隙最简洁的物理来源是原胞内存在两个不等价的格点。一维双原子链（如 A-B-A-B…）正是这一机制的最小模型：每个原胞含两个原子，分别记为 A 与 B，原胞内跃迁记为 t₁，原胞间跃迁记为 t₂，两格点的在位能分别为 ε_A 与 ε_B。当 t₁ ≠ t₂ 或 ε_A ≠ ε_B 时，布里渊区边界处打开能隙。本节围绕原胞选取、两带哈密顿量、能隙形成机制以及数值对角化与能带绘图展开。

原胞选取与基矢约定

建模第一步是明确原胞。取晶格常数 a 为原胞长度，第 n 个原胞内 A 位坐标为 x_A = n·a，B 位坐标为 x_B = n·a + d（0 < d < a）。原胞间最近邻有两类：A(n) 与 B(n) 之间距离 d，对应跃迁 t₁；B(n) 与 A(n+1) 之间距离 a−d，对应跃迁 t₂。若取 d = a/2，则 t₁ 与 t₂ 对应等距但不同化学环境，仍可不等。全文统一采用此坐标约定，正方向沿链向右。

原胞选取直接决定布洛赫相位因子。以 A 位为原胞参考点，B 位相对位移为 d，则倒空间哈密顿量的非对角元含相位 e^{−ikd} 与 e^{ik(a−d)}。若错误地把 B 位当作参考点，相位因子会整体平移，虽不改变本征值，但会改变波函数相位与绘图中的能带形状细节，因此必须在代码中固定约定。

两带紧束缚哈密顿量

在基矢 (|A,k⟩, |B,k⟩) 下，哈密顿量为 2×2 矩阵：

H(k) = [ ε_A            −t₁ − t₂·e^{−ika} ]
       [ −t₁ − t₂·e^{ika}      ε_B          ]

其中 t₁、t₂ 取正实数表示跃迁幅度，负号来自紧束缚符号约定。对角化得到两条能带：

E±(k) = (ε_A + ε_B)/2 ± sqrt( ((ε_A − ε_B)/2)² + |t₁ + t₂·e^{ika}|² )

关键量是 |t₁ + t₂·e^{ika}|² = t₁² + t₂² + 2t₁t₂·cos(ka)。当 k = ±π/a 时 cos(ka) = −1，该量取最小值 (t₁ − t₂)²；当 k = 0 时取最大值 (t₁ + t₂)²。因此能隙出现在布里渊区边界 k = ±π/a，大小为：

ΔE = 2·sqrt( ((ε_A − ε_B)/2)² + (t₁ − t₂)² )

该式清楚显示能隙的两个独立来源：格点能差 ε_A − ε_B（离子性来源）与跃迁不对称 t₁ − t₂（二聚化来源）。二者平方和开方，说明它们对能隙的贡献不能简单线性叠加。

能隙形成机制

从物理图像看，若 t₁ = t₂ 且 ε_A = ε_B，则双原子链退化为单原子链，两带在 k = ±π/a 处简并，能隙为零。当 t₁ ≠ t₂ 时，二聚化使原胞内与胞间耦合不等，布里渊区边界处两支布洛赫波发生相干散射，形成驻波，其能量因耦合不对称而劈裂，从而打开能隙。当 ε_A ≠ ε_B 时，即使 t₁ = t₂，两个子晶格能级本身错开，边界处简并被解除，同样产生能隙。实际材料中两种机制常同时存在。

数值对角化与能带绘图

下面给出可直接运行的 Python 示例，采用 NumPy 对角化 2×2 矩阵，并绘制能带。示例参数为 ε_A = 0.0、ε_B = 0.5、t₁ = 1.0、t₂ = 0.6、a = 1.0，属于合理假设的教学示例，非实测数据。

import numpy as np
import matplotlib.pyplot as plt

a = 1.0
epsA, epsB = 0.0, 0.5
t1, t2 = 1.0, 0.6

Nk = 400
k = np.linspace(-np.pi/a, np.pi/a, Nk)
E = np.zeros((Nk, 2))

for i, kk in enumerate(k):
    H = np.array([[epsA, -t1 - t2*np.exp(-1j*kk*a)],
                  [-t1 - t2*np.exp(1j*kk*a), epsB]], dtype=complex)
    E[i] = np.linalg.eigvalsh(H)

gap = E[:,1].min() - E[:,0].max()
print("数值能隙 =", gap)
print("解析能隙 =", 2*np.sqrt(((epsA-epsB)/2)**2 + (t1-t2)**2))

plt.plot(k, E[:,0], label="下带 E-")
plt.plot(k, E[:,1], label="上带 E+")
plt.axvline(np.pi/a, ls="--", c="gray")
plt.axvline(-np.pi/a, ls="--", c="gray")
plt.xlabel("k"); plt.ylabel("E(k)")
plt.legend(); plt.tight_layout(); plt.show()

运行后数值能隙应与解析式一致，约为 2·sqrt(0.25² + 0.4²) ≈ 0.943。绘图时应把 k 轴范围设为 [−π/a, π/a]，并在 ±π/a 处标注虚线，便于观察能隙位置。若发现两带在边界处交叉而非分开，应检查相位因子符号是否与基矢约定一致，或 t₁、t₂ 是否被误设为相等。

从图中可读出三点：其一，能隙位于布里渊区边界 ±π/a，而非 k = 0；其二，上带在边界处取极小、下带取极大，二者之差即能隙；其三，能带关于 k = 0 对称，这是时间反演对称性的体现。若把 ε_B 调回 0 且令 t₁ = t₂，图中两条曲线将在边界处相接，能隙闭合，可作为代码正确性的自检用例。

参数扫描与结果判读

为理解两个机制各自的贡献，可固定 t₁ = 1.0、t₂ = 0.6，扫描 ε_B 从 0 到 1.0，记录数值能隙并与解析式对比。预期能隙随 |ε_A − ε_B| 单调增大，在 ε_B = 0 时能隙为 2|t₁ − t₂| = 0.8，在 ε_B = 1.0 时约为 2·sqrt(0.25 + 0.16) ≈ 1.28。若数值结果偏离解析值超过 1e−10，应优先检查 k 网格是否覆盖边界点、矩阵是否厄米、相位因子是否与基矢一致。

完成上述步骤后，读者应能独立完成从原胞约定、哈密顿量构造、解析能隙推导到数值验证与能带绘图的完整链路，并能通过参数扫描区分二聚化与离子性两种能隙来源。