来源论文: https://arxiv.org/abs/2607.05853v1 生成时间: Jul 08, 2026 07:22

基于数值原子轨道与双 k 网格策略的高效 Bethe-Salpeter 方程（BSE）计算：ABACUS+LibRPA 理论、算法与深度复现指南

0. 执行摘要

在凝聚态物理与材料化学领域，精确预测半导体、绝缘体及低维材料的激发态性质与光学响应函数是一项核心而艰巨的任务。基于一粒子格林函数理论的 $GW$ 加上 Bethe-Salpeter 方程（$GW$+BSE）方法，因其能够系统地处理光致激发过程中电子与空穴之间的强关联库仑相互作用（激子效应），被公认为预测中性激发能和光学吸收光谱的“金标准”。然而，传统 $GW$+BSE 方法在处理周期性固体体系时面临极其严峻的计算瓶颈：构建激子波函数和精确捕捉吸收光谱的峰值和振子强度，往往需要在倒空间（布里渊区）采用极其密集的 $k$ 网格采样。若对 $GW$ 准粒子自能和 BSE 相互作用内核均采用等同的密集网格计算，其极高的计算复杂度（通常随系统尺度呈 $O(N^4)$ 甚至更高阶标度）将使得实际应用难以为继。

针对这一关键瓶颈，中国科学院物理研究所任新国研究员团队及其合作者发表了题为 “Efficient Bethe-Salpeter Equation Calculations Based on Numerical Atomic Orbitals and Norm-Conserving Pseudopotentials: Dual-k-mesh Strategy” 的重要工作。该工作基于自主开发的国产第一性原理软件 ABACUS 以及后 DFT 多体格林函数计算框架 LibRPA，巧妙地结合了数值原子轨道（Numerical Atomic Orbitals, NAOs）、局部恒等分辨技术（Localized Resolution-of-Identity, LRI）以及双 $k$ 网格（Dual-$k$-Mesh）插值策略。通过将筛选库仑相互作用 $W$ 投影至实空间局部辅助基组（ABFs）中，利用实空间矩阵元 $W_{\mu\nu}(\mathbf{R})$ 的天然局域性与短程衰减特性，实现了相互作用内核从粗 $k$ 网格到任意密集 $k$ 网格的高效傅里叶插值。这一双网格工作流能够在基本不损失计算精度的前提下，将周期性固体 BSE 计算的整体计算效率提升数个数量级。

本博客将对该项工作进行深度的理论与工程化技术解析，涵盖核心科学原理、LRI 局部化展开方法、双 $k$ 网格傅里叶插值算法流程、关键体系的收敛性行为（如 NAO 基组、辅助基组及 $k$ 点采样等）、分子与周期性固体（Si, MgO, GaN 等）的基准测试、代码实现与复现配置方案，并最终客观评述该方法的科学局限性与未来发展前沿。

1. 核心科学问题、理论基础、技术难点与方法细节

1.1 核心科学问题：BSE 内核的非局域性与密集 $k$ 网格瓶颈

在中性激发过程中，体系的光学吸收主要由电子-空穴对的协同运动决定。BSE 方程在数学上可以看作是两粒子传播子的 Dyson 方程。其相互作用内核（Kernel）$K$ 包含两项：

1. 直接库仑项（Direct Term）：描述电子与空穴之间被介质筛选的库仑吸引作用，依赖于筛选库仑相互作用 $W$；
2. 交换库仑项（Exchange Term）：描述电子与空穴之间的无筛选（裸）库仑排斥相互作用，体现了激子的自旋多重度特征。

对于三维周期性晶体，由于光致激发形成的激子（特别是 Wannier-Mott 弱束缚激子）在实空间具有较大的空间伸展（激子玻尔半径可达数十到数百埃），在倒空间中其波动特征则表现为极度局域化。这意味着，若要精确捕捉激子能级和振子强度，必须在倒空间中进行极其细致的布里渊区（BZ）采样（例如 $20 \times 20 \times 20$ 甚至更高的 $k$ 网格）。然而，筛选库仑相互作用 $W$ 的计算涉及极其昂贵的非相互作用响应函数 $\chi^0$ 的构建与矩阵求逆，这在密集的 $k$ 网格下是几乎无法承受的。传统的应对方案往往通过直接插值准粒子能级，并在较粗的 $k$ 网格上近似计算 BSE 内核，但这会引入显著的插值误差并导致吸收光谱形状的严重畸变。

1.2 理论基础与多体微扰框架

1.2.1 一粒子格林函数与 $G^0W^0$ 准粒子能级

多体微扰理论的起点是基于 Kohn-Sham (KS) 轨道构建的一粒子非相互作用格林函数 $G^0$。在实空间和虚时间表示下：

其中 $\phi_\mu$ 代表数值原子轨道，$\mathbf{R}$ 为晶格矢量，$\boldsymbol{\tau}_\mu$ 为胞内原子位置。基于 $G^0$，可顺次计算非相互作用极化函数 $\chi^0$、介电函数 $\varepsilon$、筛选库仑相互作用 $W$ 以及自能算符 $\Sigma$。在经典的 $G^0W^0$ 级次下，准粒子能级 $E^{\text{GW}}_{n\mathbf{k}}$ 通过对 KS 本征能级 $\varepsilon_{n\mathbf{k}}$ 进行一阶修正得到：

1.2.2 Bethe-Salpeter 方程的矩阵形式

在电子-空穴对表象（通常选择占据态 $v$ 和未占态 $c$ 的乘积基组 $\psi_{v\mathbf{k}} \psi^*_{c\mathbf{k}}$）下，静态近似下的完整的 BSE 矩阵方程写为：

其中 $\Omega_S$ 为激子激发能，$\begin{pmatrix} X^S \\ Y^S \end{pmatrix}$ 分别为激子态 $S$ 的正向与反向激发振幅（对应激子波函数）。其中矩阵元 $A$ 与 $B$ 的具体公式为：

这里 $\alpha$ 为自旋多重度因子（对于自旋单态 $\alpha = 2$，自旋三态 $\alpha = 0$）；$v$ 项代表裸库仑项，贡献激子自能的交换作用；$W$ 项代表静态筛选库仑项，贡献电子-空穴间的直接吸引作用。积分表达式如下：

在 Tamm-Dancoff 近似（TDA） 下，耦合矩阵 $B$ 被完全忽略，BSE 演变为标准的厄米矩阵特征值问题：

1.3 技术难点与应对策略

1.3.1 难点一：四指标库仑积分的高效展开与计算

由于 NAO 并非平面对称的基组，计算非局部轨道乘积对应的四指标积分极其复杂。传统方法不仅计算量达到 $O(N^4)$，且内存开销巨大。

应对方案（LRI 技术）：通过引入局域辅助基组（Auxiliary Basis Functions, ABFs）$P_\mu$，将空间中两个同胞或临近胞的 NAO 轨道乘积进行二体展开：

该展开系数通过在库仑度规下极小化拟合误差精确确定，从而将极其困难的四中心积分简化为辅助基组下的二中心库仑矩阵元：

1.3.2 难点二：倒空间裸库仑相互作用的 $q \to 0$ 奇点

在周期性体系中，裸库仑算符在长波极限下存在 $1/q^2$ 的发散。在 ABFs 表象中，这种奇异性表现为“头部（Head）”和“非对角翼部（Wing）”项的发散，若不加处理将严重破坏电中性体系的收敛性。

应对方案：在计算对称化介电函数和筛选库仑矩阵时，该工作通过平面对抗基投影（利用 bare Coulomb 矩阵的最大本征值特征向量作为平面的 $G=0$ 对应分量）进行解析奇点消减；同时采用 Spencer-Alavi 截断库仑势方法，在 Born-von Kármán 超胞内严格截断长程尾巴，确保了数值计算的极高稳定度。

1.4 双 $k$ 网格（Dual-$k$-Mesh）算法方法细节

双 $k$ 网格的核心物理图像建立在以下事实上：虽然准粒子激发态和激子波函数需要高精度的倒空间采样（Dense $k$-mesh），但筛选相互作用 $W(\mathbf{r}, \mathbf{r}', \omega=0)$ 在实空间具有极强的空间定域性（在半导体和绝缘体中，由于静电屏蔽作用，它表现为典型的短程性质衰减）。 这一空间定域性意味着，我们仅需在较粗的网格（Coarse $k$-mesh）上计算其在实空间的矩阵元，再通过傅里叶变换即可插值到任意密集的倒空间 $k$ 点上。

算法工作流（对应论文中的图 2）包含以下核心步骤：

[Step 1: DFT 计算 (ABACUS)]
  │ └─ 计算粗网格与细网格下的 KS 能级 ε_{nk} 与本征函数 c_{sn}^k
  │ └─ 构建实空间 LRI 展开系数 C^μ_{st}(R) 和辅助基表象库仑矩阵 V_{μν}(R)
  ▼
[Step 2: G0W0 计算 (LibRPA)]
  │ └─ 基于粗网格 (Coarse q-mesh) 计算静态筛选库仑矩阵 W_{μν}(q, iω ≈ 0)
  │ └─ 逆傅里叶变换至实空间短程胞元 W_{μν}(R)
  │ └─ 计算粗网格准粒子自能并解析延拓，插值获得细网格下的 E^{GW}_{nk}
  ▼
[Step 3: BSE 矩阵构建与插值 (ABACUS)]
  │ └─ 在细网格 (Dense k-mesh) 上，将 W_{μν}(R) 进行傅里叶插值至任意 W_{μν}(q_dense)
  │ └─ 组合细网格下的轨道本征系数 c_dense，构建细网格下的直接项与交换项内核矩阵
  ▼
[Step 4: 对角化与光谱构建 (ABACUS)]
    └─ 求解大尺寸 BSE 矩阵 A（或含 B 的全矩阵）
    └─ 获取激子能级 Ω_S 及其振幅 X_S，计算宏观介电函数 ε_2(ω)

2. 关键 Benchmark 体系、计算所得数据与性能展示

为了系统评估该理论及算法实现的精确度与效率，研究团队对分子体系（标准 Thiel 分子集）和典型周期性固体（共价半导体 Si，宽禁带离子绝缘体 MgO，典型三维极性材料 GaN）进行了多维度的基准测试。

2.1 分子测试集：Thiel’s Set 的系统对标

Thiel 分子测试集包含 28 个有机分子，是评估激发态方法精确度的经典标准。作者将基于 LRI 和 ABACUS+LibRPA 的计算结果，与同基组下的全电子高精度第一性原理软件 FHI-aims 进行了直接对比，计算包含了每个分子的前 5 个激发态（Singlet 与 Triplet，TDA 与 Full BSE 共计 140 个激发态）。

基组配置：Tier2（等同于高精度三 Zeta 极化基组 TZDP）。

表 1：ABACUS+LibRPA 与 FHI-aims 在 Thiel 分子集上的统计偏差分析（单位：meV）

关键数据结论与物理分析：

1. 极佳的吻合度：ABACUS 作为一个基于赝势的程序，其与全电子程序 FHI-aims 的平均绝对偏差（MAD）控制在 30 - 44 meV 之间，这充分验证了多体微扰矩阵元构建与 LRI 积分在数值实现上的高度一致性与严谨性。
2. 共轭长度相关性：数据分析表明，在直链聚烯烃系列（乙烯 $C_2 \to$ 丁二烯 $C_4 \to$ 己三烯 $C_6 \to$ 八四烯 $C_8$）中，随着分子共轭链长度的增加，两种代码之间的激发能偏差呈现单调温和增长趋势（例如，在三态 Full BSE 路径下，从乙烯的 $7.6 \text{ meV}$ 增大至八四烯的 $208.0 \text{ meV}$）。这反映出大离域体系中，激子波函数具有更广的空间分布，放大了不同频域展开网格（LibRPA 使用极高精度的 Minimax 时间-频率网格，而 FHI-aims 使用 Gauss-Legendre 网格）所带来的微小积分差异。

2.2 周期性体系收敛性行为评估

2.2.1 NAO 基组收敛性与“速度-长度”规范对等诊断

在理想的完备基组下，BSE 光学吸收光谱的**速度表象（Velocity Gauge）和长度表象（Length Gauge）**是严格等价的：

由于计算中使用的 NAO 基组是不完备的，这两者之间会产生偏差（即基组不完备误差 BSIE）。

研究团队利用苯分子的吸收光谱，对比了 DZP（双 Zeta 双极化）、TZDP（三 Zeta 三极化）、QZTP（四 Zeta 四极化）三种不同精度的原子轨道基组下的计算结果（如图 3 所示）。计算表明，随着基组由 DZP 拓宽到 QZTP，速度与长度规范之间的光谱曲线重合度显著提升，特别是在 $12 - 16 \text{ eV}$ 的高激发能区，QZTP 几乎消除了规范不对等性，这强有力地提供了一种内部不确定度评估手段。

2.2.2 辅助基组（ABF）与 PCA 压缩收敛性

LRI 方法的一大技术参数是主成分分析（Principal Component Analysis, PCA）的奇异值筛选阈值。阈值越小，保留的实空间辅助基（ABF）数量越多，理论精度越高，但计算负荷随之增加。

作者在 wurtzite 结构的 GaN 体系中（如图 4 所示），测试了不同 PCA 阈值（$2\times 10^{-3}, 10^{-3}, 10^{-4}, 10^{-5}$）下的吸收光谱，这些阈值对应在单个晶胞中分别产生 988, 1178, 1704, 2348 个辅助基函数。结果显示，即使在最粗糙的 $2\times10^{-3}$ 阈值下，GaN 的吸收光谱与极高精度的 $10^{-5}$ 谱线之间的绝对差值也保持在 $0.3 \text{ a.u.}$ 以下。这证明了相对于 $GW$ 准粒子计算（需要大量的极高能级非占据态因而对 ABF 要求极高），BSE 局域激子核的计算对辅助基的完整度敏感性显著降低，使得我们可以大胆采用较大的 PCA 截断阈值来极大加速计算。

2.2.3 $k$ 点采样收敛加速：Offset Averaging（偏移平均法）

周期性固体最棘手的问题是极其缓慢的 $k$ 点收敛。如图 5 所示，Si 在传统的等间距 $8\times8\times8$ 到 $14\times14\times14$ 网格下，在 $4.5 - 5.5 \text{ eV}$ 范围内的峰值极不稳定。尤其是在传统的无偏移网格中，高对称性 $k$ 点的非物理贡献会导致出现虚假的强吸收峰。

为了消除此伪影，作者实现了 Offset Averaging 方法（通过在布里渊区中采用 4 个相互独立的非 $\Gamma$ 对称性偏置网格进行独立 BSE 求解并取权重加权平均，这四个偏移量分别为：$(0,0,0)$ 权重 1/27，$(1/3,0,0)$ 权重 8/27，$(1/3,1/3,0)$ 权重 6/27，及 $(2/3,1/3,0)$ 权重 12/27）。如图 6 所示，采用偏移平均后的 $14\times14\times14$ 谱线与极度昂贵的均匀 $21\times21\times21$ 采样谱线完全重合，成功平滑了吸收曲线，展示出极强的实用收敛加速效能。

2.3 典型材料吸收光谱与激子结合能（Exciton Binding Energy）

利用该方法，作者对经典共价半导体硅（Si）与强离子型晶体氧化镁（MgO）的光谱进行了基准计算，并与经典的平面对抗程序 BerkeleyGW 及 VASP 的结果进行了仔细比对（分别见论文的图 7 与图 8）。结果表明，ABACUS+LibRPA 完美地复现了 Si 的 $E_1$ ($3.4 \text{ eV}$) 与 $E_2$ ($4.3 \text{ eV}$) 吸收双峰，且精确给出了 MgO 位于 $7.6 \text{ eV}$ 的强激子吸收主峰及其相对振子强度分布。

激子结合能 $E_B$ 的精确预测是检验 BSE 计算的最关键指标，它定义为体系最小直接准粒子带隙与最低激子本征激发能之差：

表 2：多种典型周期性半导体与绝缘体的激子结合能 $E_B$ 基准测试数据（单位：meV）

核心结论与深入分析：

1. 极强的网格依赖性与物理机理：激子结合能 $E_B$ 的收敛极为缓慢。以 GaN 为例，随着细网格从 $11\times11\times7$ 逐步加密到 $21\times21\times14$，$E_B$ 从 $74 \text{ meV}$ 急剧单调下降到 $13 \text{ meV}$（如图 9 所示）。在物理上，随着网格密度增加，倒空间分辨率提高，这允许体系更好地描述大空间尺度、小动量动能扩展的激子波函数（其物理图像上属于大玻尔半径的弱束缚 Wannier 激子）。相比之下，强激子束缚体系（如 MgO，$E_B \approx 284 \text{ meV}$）由于在实空间极度局域，在倒空间本身就大范围弥散，其网格敏感度则要低得多。
2. 与实验值的优异契合度：本工作计算得到的 Si ($17 \text{ meV}$), CdS ($29 \text{ meV}$), GaN ($13 \text{ meV}$) 的 $E_B$ 与实验测得的物理值在极窄的公差内完全契合。这一极高吻合度不仅由于双网格下插值内核消除了数值振荡，也得益于基于 LRI 的高精度静态筛选库仑矩阵构建。

3. 代码实现细节、复现指南与开源链接

3.1 软件架构与依赖链

该高效双网格 $GW$+BSE 工作流通过三方国产及开源计算生态的深度整合完成：

1. ABACUS (Atomic Band Unfolding and Computations with United Structure): 国产开源第一性原理电子结构计算软件，负责完成 DFT 自洽计算，输出多体微扰计算所需的基本实空间积分矩阵元、Kohn-Sham 本征轨道系数与能级。
2. LibRPA: 作为多体计算的核心动力学引擎，负责完成时间-空间格林函数求解、自能 $\Sigma$ 矩阵计算、准粒子修正 $G^0W^0$ 能级输出以及静态实空间筛选库仑作用元 $W_{\mu\nu}(\mathbf{R})$ 的打包。
3. LibRI: 高效并行局部 LRI 积分收缩库，提供了分布式内存下（MPI/OpenMP 混合）大规模三中心与双中心 LRI 积分的高性能计算支持。

3.2 深度复现配置方案与关键参数设置

要启动一个周期性材料（如 Si）的双网格 $GW$+BSE 计算，需要依次执行三个阶段：

第一阶段：ABACUS 下完成 Ground-state DFT 基础计算

编写 INPUT 控制文件：

INPUT_PARAMETERS
# 核心电子结构参数
system_name       Si
symmetric_force   1
calculate_force   1

# 数值原子轨道与赝势设定
basis_type        lcao
pseudo_dir        ./
orbital_dir       ./

# 倒空间网格（此时可以设定较细的 target k-mesh，用于提前准备密集的 KS 本征态）
kpoint_type       monkhorst-pack
grid_speed        1

# 必须开启的高级关联泛函相关接口配置
exx_rmesh_times   2.0      # 确定实空间辅助基 LRI 的积分半径截断，推荐设置为2.0倍的LCAO半径

设置 KPT 网格文件（作为 target 密网格）：

K_POINTS
0
Monkhorst-Pack
14 14 14 0 0 0

第二阶段：LibRPA 下计算粗网格自能 $G^0W^0$ 与静态 $W$

此时配置 LibRPA 运行环境，并读入 ABACUS 生成的 LRI 二/三中心积分。LibRPA 的输入模板中需要指定：

{
  "calculation_type": "G0W0",
  "coarse_k_mesh": [6, 6, 6],
  "dense_k_mesh": [14, 14, 14],
  "frequency_grid": {
    "type": "minimax",
    "points": 30
  },
  "read_lri_matrix": true,
  "output_static_W": true
}

此处 coarse_k_mesh 表示我们仅在 $6\times6\times6$ 级别上求解昂贵的多体极化函数。LibRPA 运算结束后，将在主目录下写入二进制的静态筛选库仑矩阵文件 W_mu_nu_R.bin，以及准粒子修正后的密网格能级数据 QP_energy_dense.dat。

第三阶段：ABACUS 端重载并执行双 $k$ 网格 BSE 求解

在进行 BSE 的主计算时，ABACUS 读取上一步的能级和实空间 $W_{\mu\nu}(\mathbf{R})$ 矩阵，通过内置的傅里叶插值器直接在 $14\times14\times14$ 网格上完成激子 Hamiltonian 的对角化：

INPUT_PARAMETERS
bse_calculation   true        # 开启 BSE 流程
bse_variant       tda         # TDA：仅对角化 A 矩阵；full：解非厄米完整 BSE
bse_spin_channel  singlet     # 计算单态激发光谱

# 导入多体修正输入
read_qp_energy    true
qp_energy_file    QP_energy_dense.dat

# 实空间 W 傅里叶插值配置
read_static_W     true
static_W_file     W_mu_nu_R.bin

# 吸收光谱计算参数
lor_broadening    0.1         # 罗伦兹展宽因子 (eV)
energy_max        15.0        # 光谱最大输出能量 (eV)

3.3 开源项目仓库链接

• ABACUS 开发主分支项目： https://github.com/abacusmodeling/abacus-develop

• LibRPA 核心计算框架： https://github.com/Srlive1201/LibRPA

• LibRI LRI 二中心与三中心积分处理引擎： https://github.com/abacusmodeling/LibRI

4. 关键引用文献与学术局限性点评

4.1 关键引用文献及科学承袭关系

1. 多体微扰与极化内核经典文献：
   • Strinati, G. (1988), “Application of the Green’s functions method to the study of the optical properties of semiconductors,” Riv. Nuovo Cimento 11, 1. (奠定了两粒子格林函数求解 BSE 的理论框架)。
   • Rohlfing, M. and Louie, S. G. (2000), “Electron-hole excitations and optical spectra from first principles,” Phys. Rev. B 62, 4927. (确立了第一性原理平面对抗 $GW$-BSE 计算的标准计算模式)。
2. 数值原子轨道基组下的 BSE 先驱工作：
   • Zhou, R., Yao, Y., Blum, V., Ren, X., and Kanai, Y. (2025), J. Chem. Theory Comput. 21, 291. (本工作的理论直系前身，首次在 FHI-aims 全电子 NAO 基组下实现了周期性体系的 BSE 方程)。
3. 不均匀 $k$ 点采样与双网格策略文献：
   • Alvertis, A. M. et al. (2023), Phys. Rev. B 108, 235117. (对比了基于非均匀布里渊区采样的激子能级收敛加速，是极具代表性的收敛性对标文献)。
   • Spencer, J. and Alavi, A. (2008), Phys. Rev. B 77, 193110. (提供了实空间截断库仑相互作用以消除 $q \to 0$ 奇点的具体数学形式)。

4.2 本工作局限性批判与科学反思

虽然本项工作在理论算法和数值重构效率上取得了突破性的进展，但作为多体微扰计算前沿方法，依然存在以下不可忽视的学术局限性：

1. 局限一：静态筛选库仑相互作用近似（Static Screening Approximation, $\omega = 0$） 本工作及常规 BSE 大多采用了传统的静态筛选近似，即假定电子-空穴相互作用过程中库仑屏蔽的响应是瞬间完成的，这忽略了高频能区动态筛选效应（Dynamical Screening）的调制作用。对于一些具有超高频声子、极性强关联特征的过渡金属氧化物，静态近似会系统性地低估准粒子能带边折叠能，并导致吸收峰位置出现轻微的高能或低能漂移。

2. 局限二：电子-声子耦合效应（Electron-Phonon Coupling）的物理缺失 理论计算的光谱与表 2 中的实验值对比仍有一些微小系统性偏差。这是因为本计算完全基于固定晶格近似（Born-Oppenheimer 近似下的 $0 \text{ K}$ 刚性结构）。在实际实验中，晶格的热振动和电声相互作用会引起准粒子带隙的温升收缩，同时声子辅助激发会给光学吸收带边带来显著的红移和不对称的谱线拓宽。未纳入电声相互作用是导致 MgO 激子带隙低估、谱线过窄的物理本因。

3. 局限三：双 $k$ 网格的“伪不均匀性”与对目标密网格的内存限制 尽管对筛选相互作用 $W$ 实现了粗网格计算，但在对角化激子 Hamiltonian $A$ 时，BSE 矩阵的总维度为 $N_v \times N_c \times N_k$。在最密集的 $21\times21\times21$ $k$ 网格下，矩阵维度呈爆发式增长，极其消耗内存，且传统的直接对角化算法的计算复杂度高达 $O(M^3)$。目前的双网格插值虽然加速了自能内核构建，但并没有减小最后一步求解大对角化矩阵的数学复杂性。若能结合真正自适应非均匀布里渊区“小斑块（non-uniform local patch）”细密采样，才能从根本上攻克高维激子 Hamiltonian 的对角化屏障。

5. 补充理论推导与进阶前沿探讨

5.1 激子 Hamiltonian 块结构形式及数学物理推导

为了展示从两粒子 Dyson 方程到本征值方程的平滑数学过渡，我们有必要详细推导激子激发算符的特征问题形式。我们在第 3 页看到了关联函数 Dyson 方程的算符表达（公式 6）：

我们令基组为 $\chi_{v c \mathbf{k}} = \psi_{v\mathbf{k}}(\mathbf{r}_1)\psi^*_{c\mathbf{k}}(\mathbf{r}_2)$。在这个两粒子乘积空间内，我们将 Dyson 方程投影到此基组上，通过使用恒等算符展开：

在无温绝热激发极限下，费米分布函数差值 $(f_p - f_o)$ 在 $v \to c$（占据态到未占态）激发时值为 $1$；在 $c \to v$（未占态到占据态）去激发时值为 $-1$。

我们将指标分类定义，让 $(i, j)$ 代表占据态价带，$(a, b)$ 代表未占态导带：

1. 当对应物理过程为：正向激发（即 $v \to c$ 对，记为物理指标 $ia$）： 其对角能量项为 $E_a - E_i > 0$。
2. 当对应物理过程为：反向去激发（即 $c \to v$ 对，记为物理指标 $ai$）： 其对角能量项为 $E_i - E_a < 0$。

此时，Dyson 矩阵算符的左侧括号项可完美拆分为四个块：

• 左上块 $A$ (激发-激发空间)： $$A_{ia, jb} = (E_a - E_i) \delta_{ij}\delta_{ab} + K_{ia, jb}$$
• 右上块 $B$ (激发-去激发空间)： $$B_{ia, jb} = K_{ia, bj}$$
• 左下块 $-B^*$ (去激发-激发空间)： 由于费米差因子为 $-1$，矩阵元带负号。
• 右下块 $-A^*$ (去激发-去激发空间)： 对角项为 $E_i - E_a$，整体矩阵元带负号。

这就完美推导出了公式 19 所呈现的非厄米系统本征值矩阵：

该形式完美体现了由于去激发通道的存在，两粒子体系在数学上与随时间变化的含时密度泛函理论（TDDFT）具有高度一致的本征值代数结构（即著名的 Casimir-Casimir 关联结构）。

5.2 物理图像探讨：为什么激子结合能对 $k$ 网格极其敏感？

我们可以利用量子力学不确定性原理，定性理解激子结合能对 $k$ 网格的收敛规律。Wannier 激子表现为弱束缚电子-空穴对，两者在实空间具有极其开阔的相对运动空间分布。其波函数包络满足氢原子模型：

其中 $a_B^*$ 为激子的有效玻尔半径。在实空间，激子波函数的特征伸展度为 $a_B^*$。根据傅里叶变换的量子力学不确定性关系：

在倒空间中，该激子波函数所对应的布里渊区动量动能展开范围为 $\Delta k \sim 2\pi / a_B^*$。

1. 当体系激子结合能很低时（例如单晶硅 Si，$E_B \approx 15 \text{ meV}$）： 激子玻尔半径 $a_B^*$ 极大（通常超过数十纳米，跨越数十个晶胞）。此时，其在倒空间中的动量宽度 $\Delta k$ 变得极小，所有物理权重被高度压缩、局域在 $\Gamma$ 点附近微小的倒空间区域。如果倒空间离散 $k$ 网格不够细密，则在此极窄带宽范围内的物理状态无法被网格“点样积分”捕捉，积分被严重低估，从而导致巨大的数值收敛不稳定性。
2. 当体系激子结合能很高时（例如 MgO，$E_B \approx 80 - 145 \text{ meV}$）： 激子在实空间极小，在倒空间极为弥散（$\Delta k$ 很大）。此时采用较粗的 $k$ 网格就足以精确采样其动量谱特征，因此对网格加密不敏感。

双 $k$ 网格插值方案正是这一物理本质的绝佳工程解。它在粗网格上计算短程的筛选相互作用（该作用随 $r$ 衰减极快，对应的 $k$ 空间谱线极平滑，易用粗网格捕捉），而把长程、极窄的动量转换激子包络完全留在细网格的傅里叶插值阶段解析处理。这一优美而统一的物理图像正是本篇论文的精髓所在。

5.3 进阶前沿应用：有限动量激发与动量解析光谱

在目前的主流 BSE 实现中，人们通常只关注光子的长波极限（$q \to 0$），因为可见光的动量相对电子的布里渊区尺度极小，属于垂直光学跃迁（如图 1 所示）。然而，前沿实验物理学（如非弹性 X 射线散射 (IXS)、电子能量损失谱 (EELS)）能够实现高动量转移，探测处于有限动量状态（Finite-momentum $q \neq 0$）的非垂直激子跃迁。

在双 $k$ 网格框架下，向 $q \neq 0$ 激发谱扩展的理论路径是十分清晰的：仅需将插值转换算符中的静态筛选库仑矩阵 $W_{\mu\nu}(\mathbf{R})$ 乘以空间相位因子 $e^{i\mathbf{q}\cdot\mathbf{R}}$，即可自然重构带有限动量位移的非局域 BSE 内核。这一方向将是材料科学预测声子辅助间接带隙激子效应及动量解析多体跃迁谱的一大利器。