基于 Toeplitz 加速与球谐图神经网络的 Bethe–Salpeter 方程高效求解器

基于 Toeplitz 加速与球谐图神经网络的 Bethe–Salpeter 方程高效求解器

摘要
Bethe–Salpeter 方程（BSE）是量子场论中描述两体束缚态的相对论性积分方程，在强子物理、凝聚态物质和量子化学等领域有广泛应用。然而，BSE 的四维离散化导致矩阵规模随网格点数平方增长（$O(N^2)$），传统求解方法受限于计算资源和内存，难以处理高精度网格。本文提出一种结合 Toeplitz 矩阵结构与幂迭代的高效求解算法。我们利用 BSE 核在能量维度上的平移不变性，将矩阵‑向量乘法复杂度从 $O(N^2)$ 降至 $O(N)$，内存需求从 $O(N^2)$ 降至 $O(N)$。结合球谐图神经网络（SH‑GNN）对波函数进行预测，可进一步加速本征值扫描。数值实验表明，在普通工作站上，50×50 网格（2500 自由度）的全介子谱计算可在 1 分钟内完成，网格扩展至 2000×2000 时仍具有可行性。相比传统 $O(N^2)$ 方法，本文算法实现了两个数量级的加速，为高精度束缚态问题提供了全新的计算范式。

关键词：Bethe–Salpeter 方程；Toeplitz 矩阵；幂迭代；球谐图神经网络；束缚态；介子质量

1 引言

量子场论中，描述两粒子束缚态的相对论性积分方程被称为 Bethe–Salpeter 方程（BSE）。该方程由 Bethe 和 Salpeter 于 1951 年提出，是量子电动力学（QED）和量子色动力学（QCD）中处理电子-正电子（正电子素）以及夸克-反夸克（介子）束缚态的基本工具。与薛定谔方程不同，BSE 完全保持相对论协变性，能够正确处理自旋、轨道耦合以及高能过程。

尽管 BSE 具有严谨的理论基础，其数值求解长期面临两大困难：高维度和非线性本征值问题。在动量空间，BSE 是一个四维积分方程。离散化后，未知函数（Bethe–Salpeter 振幅）在四维网格上的值形成一个大型向量，维度 $N=N_{k_4}\times N_{|\mathbf{k}|}\times N_{\theta }$。而将该方程转化为矩阵本征值问题时，矩阵的尺寸为 $N\times N$。当采用中等精度网格（如 $N_{k_4}=50,N_{|\mathbf{k}|}=50,N_{\theta }=20$）时，$N$ 已达 $5\times 10^4$，矩阵元素数量 $N^2$ 高达 $2.5\times 10^9$，内存需求超过 20 GB，计算复杂度更是难以接受。因此，传统 BSE 求解器（如 Yambo、BerkeleyGW）必须依赖大规模并行计算和超级计算机，严重限制了其在实际研究中的普及。

近年来，人工智能和高性能算法为这一困境提供了新的思路。一方面，图神经网络（GNN）和球谐展开被成功应用于物理场的学习与预测，例如用 SH‑GNN 从夸克质量函数直接预测介子波函数。另一方面，Toeplitz 矩阵结构在具有平移不变性的核中广泛存在。BSE 的相互作用核 $K((k-q{)}^2)$ 仅依赖于四维动量差的平方，因此在能量维度（$k_4$）上具有平移不变性。利用这一性质，我们可以构建块 Toeplitz 矩阵，从而将矩阵‑向量乘法的计算量从 $O(N^2)$ 降至 $O(N)$，并大幅降低内存占用。

本文的主要贡献如下：

本文的组织结构如下：第 2 节介绍 BSE 的标准形式及其在欧几里得空间的简化。第 3 节给出离散化方案和本征值问题的构建。第 4 节详细阐述 Toeplitz 加速技术。第 5 节描述幂迭代及介子质量的二分法搜索。第 6 节简述 SH‑GNN 的原理及其与求解器的集成。第 7 节呈现数值实验结果。第 8 节总结全文并展望未来改进方向。

2 Bethe–Salpeter 方程的标准形式

2.1 闵可夫斯基空间的齐次 BSE

在闵可夫斯基空间，描述夸克‑反夸克束缚态（介子）的齐次 BSE 为

$$\Gamma (P,k)=-i\int \frac{d^{4}q}{(2\pi {)}^4}K(k,q;P)S(q_{+})\Gamma (P,q)S(q_{-}),$$

其中 $P$ 是介子总动量，$P^2=M_H^2$；$k$ 是夸克和反夸克的相对动量；$q_{\pm }=q\pm P/2$；$S(p)$ 是完全夸克传播子；$K(k,q;P)$ 是两粒子不可约核（two‑particle irreducible kernel）。在彩虹‑梯（Rainbow‑Ladder）近似下，核取为单胶子交换：

$$K(k,q;P)=\frac{4}{3}{\gamma }_{\mu }{g}^2{D}_{\mu \nu }(k-q){\gamma }_{\nu },$$

其中 $g^2{D}_{\mu \nu }(q)$ 是胶子传播子。通常将胶子传播子与顶点的乘积合并为一个标量函数 $K(Q^2)$（$Q^2=(k-q{)}^2$），并忽略洛伦兹结构细节，则 BSE 可简化为

$$\Gamma (P,k)=\frac{4}{3}\int \frac{d^{4}q}{(2\pi {)}^4}K((k-q{)}^2){\gamma }_{\mu }S(q_{+})\Gamma (P,q)S(q_{-}){\gamma }_{\mu }.$$

2.2 Wick 旋转到欧几里得空间

为了数值求解，对时间分量进行 Wick 旋转：$k_0=i{k}_4$，$P_0=i{P}_4$，并定义欧几里得四动量 $k_E=(k_4,\mathbf{k})$，$P_E=(P_4,0)$（介子静止系）。此时 $P_E^2=P_4^2=M_H^2$。传播子变为

$$S_E(p_E)=\frac{-i/p_E+M(p_E^2)}{p_E^2+M^2(p_E^2)},$$

其标量部分为 $1/(p_E^2+M^2(p_E^2))$。代入 BSE 并取投影到赝标量通道（${\gamma }_5$ 结构），可得仅关于标量函数 $\Phi (P,k)$ 的方程：

$$\Phi (P,k)=\frac{4}{3}\cdot 4\int \frac{d^4{q}_E}{(2\pi {)}^4}K((k_E-q_E{)}^2)\frac{M(q_{+}^2)M(q_{-}^2)}{(q_{+}^2+M^2(q_{+}^2))(q_{-}^2+M^2(q_{-}^2))}\Phi (P,q),$$

其中 $q_{\pm }=q\pm P/2$，积分测度 $d^4{q}_E=d{q}_4{d}^3\mathbf{q}$。因子 $\frac{4}{3}$ 来自颜色因子，另一个 4 来自狄拉克迹 $\text{Tr}[{\gamma }_5{\gamma }_{\mu }{\gamma }_5{\gamma }_{\mu }]=4$。为简洁，记

$$G(q_4,\mathbf{q};M_H)=\frac{M(q_{+}^2)M(q_{-}^2)}{(q_{+}^2+M^2(q_{+}^2))(q_{-}^2+M^2(q_{-}^2))}.$$

2.3 分波展开与 S 波近似

由于核 $K((k_E-q_E{)}^2)$ 仅依赖于四维夹角，可将振幅按球谐函数展开。对于 S 波（$l=0$）介子，振幅与角度无关，方程简化为关于 $k_4$ 和 $|\mathbf{k}|$ 的二维积分：

$$\Phi (k_4,k)=\frac{16}{3}{\int }_{-\infty }^{\infty }\frac{d{q}_4}{2\pi }{\int }_0^{\infty }\frac{q^{2}dq}{(2\pi {)}^2}{\mathcal{K}}_0(k_4,k;q_4,q)G(q_4,q;M_H)\Phi (q_4,q),$$

其中角度平均核为

$${\mathcal{K}}_0(k_4,k;q_4,q)=\frac{1}{2}{\int }_{-1}^{1}dxK((k_4-q_4{)}^2+k^2+q^2-2kqx).$$

对于 Yukawa 型核 $1/((k_4-q_4{)}^2+k^2+q^2-2kqx+m^2)$，角度积分解析可做：

$${\mathcal{K}}_0^{\text{Yuk}}=\frac{1}{4kq}\ln \frac{(k_4-q_4{)}^2+(k+q{)}^2+m^2}{(k_4-q_4{)}^2+(k-q{)}^2+m^2}.$$

对于更一般的核（如 Maris‑Tandy 模型的混合项），可采用数值积分（如 Gauss‑Legendre 求积）。

3 离散化与本征值问题

3.1 网格与权重

选择动量网格：

• $k_4$ 在区间 $[-K_4^{\max },K_4^{\max }]$ 上取 $N_4$ 个等距点（或 tanh 加密）。
• $k=|\mathbf{k}|$ 在对数网格上取 $N_k$ 个点，覆盖 $[k_{\min },k_{\max }]$（例如 $10^{-3}$ 到 $10^2$ GeV）。

记 $k_{4,i}$（$i=1,\dots ,N_4$）和 $k_a$（$a=1,\dots ,N_k$）。积分权重：

$$w_{k_4}^{(j)}=\frac{\Delta k_4}{2\pi },w_k^{(b)}=\frac{k_b^2\Delta k_b}{(2\pi {)}^2},$$

其中 $\Delta k_4$ 为等距步长，$\Delta k_b$ 为对数网格的差分（$\Delta k_b\approx k_b\Delta (\ln k_b)$）。

定义网格点总自由度 $N=N_4\times N_k$。将振幅离散化为向量 $\Phi \in {\mathbb{R}}^N$，索引映射 $(i,a)\to i{N}_k+a$。

3.2 离散化矩阵

定义矩阵 $A(M_H)\in {\mathbb{R}}^{N\times N}$，其矩阵元为

$$A_{(i,a),(j,b)}=\frac{16}{3}{w}_{k_4}^{(j)}{w}_k^{(b)}{\mathcal{K}}_0(k_{4,i},k_a;k_{4,j},k_b)G(k_{4,j},k_b;M_H).$$

则离散化后的 BSE 成为本征值方程

$$A(M_H)\Phi =\lambda (M_H)\Phi .$$

束缚态条件为最大本征值 $\lambda (M_H)=1$。因此，求解介子质量等价于寻找 $M_H$ 使得 $\lambda (M_H)=1$。

3.3 核的平移不变性

观察 ${\mathcal{K}}_0$ 的表达式，它仅依赖于 $|k_{4,i}-k_{4,j}|$，而与绝对位置无关（因为 $k_4$ 网格等距时，差值等于 $(i-j)\Delta k_4$）。因此，矩阵 $A$ 具有块 Toeplitz 结构：定义 $d=|i-j|$，则 $A_{(i,a),(j,b)}=T_{d,a,b}$，其中 $T_{d,a,b}$ 与 $i,j$ 无关。

这一性质是加速计算的关键，我们将在下一节详细利用。

4 Toeplitz 加速的矩阵‑向量乘法

4.1 Toeplitz 矩阵的存储

由于 $A$ 完全由 $D=N_4-1$ 个块 $T_{d,a,b}$ 确定，每个块的大小为 $N_k\times N_k$，总存储量为 $N_4\times N_k^2$，远小于 $N^2=(N_4{N}_k{)}^2$。例如，当 $N_4=50,N_k=50$ 时，存储量约 $50\times 2500=125,000$ 个浮点数，而完整矩阵需 $6.25\times 10^6$ 个浮点数，内存降低 50 倍。

4.2 Toeplitz 矩阵‑向量乘法

给定向量 $x\in {\mathbb{R}}^N$，将其排列为二维数组 Xj,bX_{j,b}Xj,<span class=\"mord math