4. 周期晶格傅里叶变换与动量空间哈密顿量

在实空间的紧束缚模型中,我们面对的是一个 N×N 的哈密顿量矩阵(N 为格点数),直接对角化虽然可行,但计算量随 N 增长迅速,而且对角化后得到的本征态是各格点上的振幅分布,物理图像不够直观。周期晶格的关键优势在于平移对称性——正是这一性质让我们可以通过傅里叶变换把实空间哈密顿量转换到动量空间,使矩阵维度从 N×N 降为 m×m(m 为每个原胞内的轨道/子格数),从而以极低的代价获得能带色散关系 E(k)。本节围绕三个核心环节展开:布洛赫基的构造、傅里叶变换的具体实现、以及动量空间矩阵元的计算。

一、布洛赫基的构造

设晶格具有平移不变性,原胞由格矢 R 标记,每个原胞内有 m 个轨道(或子格),轨道指标记为 α = 1, 2, …, m。实空间基矢为 |R, α⟩,表示电子局域在第 R 个原胞的第 α 个轨道上。布洛赫定理告诉我们,在周期势场中,哈密顿量与平移算符对易,因此可以将基矢取为平移算符的本征态。布洛赫基定义为:

|ψα(k)⟩ = (1/√Nc) ΣR eik·R |R, α⟩

其中 Nc 为原胞总数,k 取第一布里渊区内的离散值(对有限周期晶格)或连续值(热力学极限)。这个变换的本质是把“格点位置”这一自由度通过平面波因子 eik·R 编码为相位,而保留原胞内的轨道指标 α 作为矩阵的索引。换句话说,布洛赫基把 N×N 的大矩阵块对角化为 Nc 个 m×m 的小矩阵,每个 k 点对应一个。

在数值实现中,我们通常不需要显式构造布洛赫基向量,而是直接利用傅里叶变换关系写出动量空间哈密顿量的矩阵元。但理解布洛赫基的结构对于正确实现变换至关重要——尤其是当原胞内包含多个轨道时,相位因子只与格矢 R 有关,与原胞内的轨道指标 α 无关。

二、傅里叶变换的实现

实空间紧束缚哈密顿量的矩阵元为 Hαβ(R) = ⟨0, α|H|R, β⟩,表示原胞 0 的 α 轨道与原胞 R 的 β 轨道之间的跃迁积分。动量空间哈密顿量由傅里叶变换给出:

Hαβ(k) = ΣR Hαβ(R) eik·R

注意这里的求和遍历所有使 Hαβ(R) 非零的格矢 R。对于厄米哈密顿量,需满足 Hβα(−R) = Hαβ(R)*,这一关系在实现中应作为正确性检查。以一维链为例,若只有最近邻跃迁 −t,则 H(k) = −t(eika + e−ika) = −2t cos(ka),这是最经典的色散关系。对于更复杂的模型(如包含次近邻跃迁或复数跃迁积分),只需在求和中逐项累加即可。

数值实现时,通常按以下流程操作:首先从实空间模型定义中提取所有非零的 Hαβ(R) 项,存入一个列表,每项包含 (α, β, R, 复数值);然后对每个目标 k 点,初始化一个 m×m 的零矩阵,遍历列表,将 Hαβ(R)·eik·R 累加到矩阵的 (α, β) 位置;最后对矩阵做厄米化处理(若数值误差导致微小非厄米性)并对角化。下面给出一个可参考的 Python 实现片段:

import numpy as np

def hamiltonian_k(k, hoppings, m):
    """hoppings: list of (alpha, beta, R_vector, amplitude)"""
    Hk = np.zeros((m, m), dtype=complex)
    for alpha, beta, R, amp in hoppings:
        phase = np.exp(1j * np.dot(k, R))
        Hk[alpha, beta] += amp * phase
    # 厄米化,消除数值误差
    Hk = 0.5 * (Hk + Hk.conj().T)
    return Hk

# 示例:一维两轨道模型,原胞内两个子格 A/B
a = 1.0
hoppings = [
    (0, 1, np.array([0.0]), -1.0),   # A(0)-B(0) 胞内跃迁
    (0, 1, np.array([-1.0]), -1.0),  # A(0)-B(-1) 胞间跃迁
]
m = 2
k_points = np.linspace(-np.pi/a, np.pi/a, 201)
bands = np.array([np.linalg.eigvalsh(hamiltonian_k(k, hoppings, m))
                  for k in k_points])

这段代码的核心逻辑是:对每个 k 点,遍历跃迁列表,用相位因子 eik·R 调制实空间跃迁积分,累加得到 m×m 的 H(k),再对角化得到 m 条能带。注意示例中 A(0)-B(-1) 的跃迁对应 R = −a,相位因子为 e−ika,这与 A(0)-B(0) 的胞内跃迁共同构成两轨道模型。实际使用时,跃迁列表应直接从实空间模型的参数定义中生成,避免手工遗漏。

三、动量空间矩阵元的计算与验证

得到 H(k) 后,对角化即可获得本征值 En(k) 和本征矢 |un(k)⟩(周期部分)。本征矢的分量 un,α(k) 描述了第 n 条能带在 k 点处各轨道的权重分布,是后续计算贝里相位、有效质量、光吸收矩阵元等物理量的基础。数值上应关注以下几点:

  • 厄米性检查:H(k) 必须满足 Hαβ(k) = Hβα(k)*。若偏差超过 10−12 量级,说明跃迁列表的对称性录入有误。
  • 周期性检查:H(k + G) = H(k),其中 G 为倒格矢。可在布里渊区边界处验证。
  • 时间反演对称性:若无磁场且无磁性,En(k) = En(−k),可作为额外校验。
  • 与实空间对角化对比:对有限尺寸(如 20 个原胞)的实空间哈密顿量做数值对角化,其本征值应与动量空间在对应离散 k 点上的本征值一致(允许边界效应导致的微小偏差)。

下图展示了从实空间跃迁参数到动量空间能带的完整数据流,以及关键验证节点:

实空间 → 动量空间变换流程与验证节点 实空间跃迁参数 H_αβ(R) 列表 (α, β, R, 复数值) 例:A(0)-B(0): −t 傅里叶变换 H(k) 矩阵构造 H_αβ(k) = Σ H_αβ(R)·e^{ik·R} m×m 复矩阵 逐 k 点独立计算 对角化 能带与本征矢 E_n(k), |u_n(k)⟩ m 条能带 轨道权重分布 关键验证节点 厄米性检查 H(k) = H(k)† 偏差 < 1e−12 周期性检查 H(k+G) = H(k) 布里渊区边界验证 时间反演对称 E_n(k) = E_n(−k) 无磁场时成立 实空间对比 有限尺寸对角化 本征值一致性 示例:一维两轨道模型,最近邻跃迁 −t,原胞内 A/B 子格 H(k) = [[0, −t(1+e^{−ika})], [−t(1+e^{ika}), 0]] E_±(k) = ±t|1+e^{−ika}| = ±2t|cos(ka/2)|

上图中四个验证节点分别对应变换流程的不同阶段:厄米性检查针对 H(k) 矩阵构造的正确性,周期性检查验证傅里叶求和的完备性,时间反演对称性检查物理模型的合理性,实空间对比则是最直接的端到端校验。在实际仿真中,建议将这四个检查写成自动化断言,每次修改跃迁参数后自动运行,避免因手工录入错误导致后续物理量计算全部失效。

当原胞内轨道数 m 较大时(如 m ≥ 4),建议将跃迁列表按 (α, β) 分组存储,每组内再按 R 排序,这样在遍历 k 点时可以利用向量化操作一次性计算所有相位因子,显著提升效率。对于二维或三维晶格,k 点通常取布里渊区内的均匀网格,网格密度需保证能带色散曲线的光滑性,一般每个方向不少于 50 个点。

总结来说,布洛赫基构造提供了从实空间到动量空间的理论桥梁,傅里叶变换给出了矩阵元的具体计算公式,而数值实现的核心在于正确维护跃迁列表、逐 k 点累加相位因子、并对结果做系统性验证。掌握这一流程后,后续的能带绘制、态密度计算、拓扑不变量求解都可以在此基础上直接展开。