Bethe-Salpeter Equation Theory

• Bethe-Salpeter Equation Theory is a covariant, non-perturbative framework for modeling two-body dynamics and bound states in quantum systems.
• It employs an integral equation with an interaction kernel that defines bound state poles and captures key aspects of nuclear, QED, and condensed-matter interactions.
• The theory supports various computational approaches, from Salpeter reductions to quantum annealing techniques, enhancing simulations across multiple fields.

Searching arXiv for recent and foundational Bethe-Salpeter equation papers to ground the article. The Bethe-Salpeter equation (BSE) is a relativistic, covariant, and non-perturbative framework for two-body dynamics in quantum field theory. Originally developed for bound states in nuclear physics, it now appears in several closely related forms: as an integral equation for the two-particle Green function, as a homogeneous bound-state equation whose poles determine masses or excitation energies, and as a matrix eigenvalue problem for neutral or charged excitations in many-electron systems (Blase et al., 2020, Mátyus et al., 2022). Across hadron physics, QED bound states, condensed-matter spectroscopy, dynamical mean-field theory, and amplitude-based gravity, the common structure is an interaction kernel acting on a two-particle propagator or amplitude, with physical content controlled by the choice of kernel, truncation, and spectral analysis (Adamo et al., 2022).

1. Field-theoretic formulation and channels

In the space-time formulation, the central object is the interacting two-particle Green function $G(x_1,x_2;x_1',x_2')$. Bethe and Salpeter showed that it satisfies an integral equation built from one-particle propagators and an irreducible interaction kernel $K$. For bound states, poles of this Green function define the Bethe-Salpeter wave function, and the corresponding wave equation may be written schematically as

$$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$

(Mátyus et al., 2022).

In many-body perturbation theory, especially for optical excitations, the same content is expressed through the four-point susceptibility $L(1,2;1',2')$, defined as the functional derivative of the one-body Green function with respect to an external perturbation. The BSE then takes Dyson form,

$L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$

with $L_0 = GG$ and a kernel $\Xi^{\mathrm{BSE}}$ that, in the standard static form, combines bare exchange and screened direct interaction (Blase et al., 2020). This formulation makes explicit that the BSE is not restricted to bound-state spectroscopy in the narrow hadronic sense; it is a general two-particle response equation.

Different time orderings select different physical channels. The conventional electron-hole channel yields neutral excitations. A different ordering leads to a particle-particle propagator whose poles are double ionization potentials and double electron affinities. In that case, a BSE can be derived using a pairing field and anomalous propagators, with kernels built from $GW$, $T$-matrix, or second-Born self-energies (Marie et al., 2024). This suggests that “Bethe-Salpeter equation theory” is best understood as a family of kernel-based two-particle equations rather than a single universal formula.

2. Relativistic bound-state reductions and equal-time equations

The full BSE is four-dimensional and often too difficult for direct analysis. A standard reduction is the instantaneous approximation, in which the kernel is taken to be independent of the relative time and the constituents propagate as free particles, possibly with effective constituent masses. This yields the Salpeter equation, a three-dimensional eigenvalue problem,

$\mathcal{H}\Phi(\vec p)=M\Phi(\vec p),$

where $K$0 contains the confining kernel, Dirac structure, and effective masses (Lucha, 2010). For harmonic-oscillator confinement, the Salpeter equation further reduces to radial differential equations, making analytic scrutiny possible (Lucha, 2010).

The distinction between instantaneous and non-instantaneous kernels is not merely technical. The instantaneous approximation provides tractability, but it can also introduce unphysical instabilities if the Lorentz structure of the kernel is chosen carelessly. This is one of the central lessons of Salpeter-equation analysis: simplifying the time dependence does not remove the need for a controlled spectral theory (Lucha, 2010).

In QED bound-state theory, an exact equal-time variant of the BSE has been reviewed for atoms and molecules. After separating the relative-time variable and integrating the dominant instantaneous interaction, one obtains an exact equal-time equation

$K$1

with a corresponding no-pair approximation

$K$2

Here $K$3 and $K$4 are one-particle Dirac Hamiltonians, $K$5 is the instantaneous interaction, $K$6 are positive- and negative-energy projectors, and $K$7 contains retardation, pair, crossed-photon, and radiative corrections (Mátyus et al., 2022). In this setting, the BSE becomes the origin of a relativistic wave equation suitable for computational QED.

3. Kernels, symmetries, and nonperturbative gauge-theory implementations

In continuum QCD, the BSE is typically solved together with the quark Dyson-Schwinger gap equation. The dressed quark propagator, gluon propagator, quark-gluon vertex, and quark-antiquark scattering kernel are interdependent. A systematic approach uses the longitudinal and transverse Ward-Green-Takahashi identities to constrain the vertex structure, confirms the Ball-Chiu vertex as the correct longitudinal part, and generates anomalous chromomagnetic moment terms whose strength is proportional to the magnitude of dynamical chiral symmetry breaking (Qin, 2016). The color-singlet vector and axial-vector Ward-Green-Takahashi identities then relate the dressed quark-gluon vertex and the Bethe-Salpeter kernel, permitting symmetry-preserving truncations beyond rainbow-ladder without violating gauge symmetry or chiral symmetry (Qin, 2016).

A complementary implementation appears in Coulomb gauge at leading order. There the input interaction is a pure linear rising potential supplemented by a contact term arising from conservation of total color charge,

$K$8

Although this interaction is infrared singular, the pseudoscalar and vector BSEs can be written in terms of manifestly finite functions. The resulting equations display both dynamical chiral symmetry breaking and the leading-order heavy-quark limit (Watson et al., 2012).

Gauge-field bound states can be treated directly as well. A continuum Landau-gauge Yang-Mills formulation derives a coupled BSE for gluon-gluon and ghost-antighost components from the two-particle irreducible effective action. In this framework, ghost contributions are required by the gauge-fixed theory rather than added ad hoc. Using propagators solved on the complex momentum plane, the scalar glueball $K$9 is found at $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$0 GeV, in agreement with lattice gauge theory, whereas the pseudoscalar $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$1 emerges at $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$2 GeV, much higher than lattice values, indicating strong sensitivity to the three-gluon vertex approximation (Sanchis-Alepuz et al., 2015).

4. Spectral theory, confinement, and unstable systems

A recurring issue in instantaneous BSEs is that confining potentials do not automatically imply stable bound states. For the Salpeter equation, physically meaningful states require a real discrete spectrum bounded from below, together with self-adjointness of the relevant operator. Spectral analysis based on operator inequalities and the minimum-maximum principle shows that time-component Lorentz-vector kernels yield stable bound states for harmonic-oscillator confinement, whereas Lorentz-scalar and Lorentz-pseudoscalar kernels are problematic. For mixed kernels, stability is guaranteed if the vector component dominates the scalar, regardless of relative sign (Lucha, 2010, Lucha, 2010).

These results are important because they isolate a common misconception: “confinement” at the level of the potential does not by itself guarantee a stable hadronic spectrum in the reduced equation. The Lorentz structure of the kernel is part of the definition of a confining interaction in the Salpeter framework, not a secondary modeling choice (Lucha, 2010).

The standard homogeneous BSE also has a well-defined limitation: it describes stationary bound states, not resonances. A development aimed at unstable systems extends the formalism by considering time evolution under the total Hamiltonian and by generalizing Mandelstam’s approach to matrix elements between bound states at arbitrary final-state energy. In this resonance framework, the forward amplitude enters a dispersion relation,

$$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$3

so that the pole position acquires both a mass shift and an imaginary part associated with the total width (Chen et al., 2023). The introduction of “extended Feynman diagrams” is intended to make the analytic continuation over open and closed channels explicit (Chen et al., 2023).

A different spectral development uses Källén-Lehmann representations to solve the BSE directly in the timelike domain. In three-dimensional broken-phase scalar $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$4 theory, spectral Dyson-Schwinger equations provide dressed propagators and vertices, and the homogeneous BSE is solved as an eigenvalue problem $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$5 with on-shell condition $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$6. In the infinite-coupling limit, the lowest bound-state mass approaches $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$7, in very good agreement with the lattice value $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$8, while tree-level kernels produce tachyonic solutions (Eichmann et al., 2023).

5. Electronic excitations, susceptibilities, and molecular properties

In condensed-matter physics and quantum chemistry, the BSE is the state-of-the-art approach for neutral excitations from first principles. Combined with the $$\Psi_T(x_1,x_2) = -i \int dy_1\,dy_2\,dy_1'\,dy_2'\, G_1(x_1,y_1)G_2(x_2,y_2)K(y_1,y_2;y_1',y_2')\Psi_T(y_1',y_2')$$9 approximation for quasiparticle energies, it is routinely written as a non-Hermitian matrix problem of the form

$L(1,2;1',2')$0

where $L(1,2;1',2')$1 and $L(1,2;1',2')$2 contain quasiparticle energy differences, bare Coulomb terms, and screened Coulomb terms (Blase et al., 2020). For molecular singlet excitations, the combination $L(1,2;1',2')$3+BSE yields a typical error of $L(1,2;1',2')$4–$L(1,2;1',2')$5 eV, with a computational cost similar to time-dependent density-functional theory using hybrid functionals, but without the need for empirical exchange-correlation functional choice. At the same time, the standard static kernel does not describe double and higher excitations, triplet energies tend to be underestimated, and analytical gradients are not yet robustly available (Blase et al., 2020). Recent work extends the $L(1,2;1',2')$6-BSE toolbox to molecular properties beyond excitation energies, including coupled-perturbed response, excited-state dipole moments, gradients, and NMR spin-spin couplings, while retaining static screening as the practical choice for most current applications (Holzer et al., 12 Feb 2025).

The static approximation can be relaxed perturbatively. A renormalized first-order dynamical correction computed with full frequency-dependent random-phase-approximation screening, and explicitly beyond the plasmon-pole approximation, systematically red-shifts singlet and triplet excitation energies. For valence $L(1,2;1',2')$7 and $L(1,2;1',2')$8 transitions, the correction is $L(1,2;1',2')$9–$L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$0 eV; for charge-transfer and Rydberg excitations it is only a few hundredths of an eV. On benchmark sets, the mean absolute error is reduced from $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$1 to $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$2 eV for singlets and from $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$3 to $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$4 eV for triplets (Loos et al., 2020).

Several developments target the cost of large-scale calculations. A finite-field BSE formulation computes screened Coulomb integrals from self-consistent density differences, avoiding explicit dielectric matrices, virtual states, and even the random-phase approximation if desired; with localized orbitals obtained by recursive bisection, the screened-exchange workload scales as $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$5 and can be applied to water and ice samples with hundreds of atoms (Nguyen et al., 2019). A stochastic approach evaluates the action of the effective Coulombic interaction $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$6 by time-dependent Hartree propagation of only $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$7 stochastic orbitals, reaching at most cubic scaling and enabling spectra and exciton densities for a $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$8-electron carbon-nanohoop–fullerene system using less than $L = L_0 + L_0\,\Xi^{\mathrm{BSE}}\,L ,$9 core hours (Bradbury et al., 2022).

The BSE also appears in two-particle response beyond optical absorption. Within dynamical mean-field theory, the dual Bethe-Salpeter equation separates local and non-local contributions to the lattice susceptibility and rewrites the response in terms of the fully reducible impurity vertex $L_0 = GG$0,

$L_0 = GG$1

Because the non-local bubble decays as $L_0 = GG$2, the dual formulation achieves cubic convergence with respect to the Matsubara cutoff rather than the linear convergence of the conventional approach, and uses a vertex free from irreducible-vertex divergences (Loon et al., 2023). In a different channel, the particle-particle BSE derived with anomalous propagators gives access to double ionization potentials and double electron affinities, with $L_0 = GG$3, $L_0 = GG$4-matrix, and second-Born kernels assessed for singlet and triplet valence double ionization and for double core-hole states (Marie et al., 2024).

6. Matrix structure, algorithms, and expanding frontiers

Once discretized, many BSEs become structured eigenvalue problems. For crystalline systems, the Hamiltonian has the form

$L_0 = GG$5

with a positive-definite BSE Hamiltonian $L_0 = GG$6. This guarantees real eigenvalues in $L_0 = GG$7 pairs and permits reformulation in terms of the product $L_0 = GG$8, where $L_0 = GG$9 and $\Xi^{\mathrm{BSE}}$0. On that basis, Cholesky and Cholesky-plus-SVD solvers were developed: the CHOL method requires less computational effort while retaining the same degree of accuracy as the square-root approach, and the CHOL+SVD method improves the expected accuracy, especially for small eigenvalues relevant to absorption edges and excitons (Benner et al., 2020).

On the Matsubara axis, two-particle objects can be compressed in the intermediate representation. Because the singular values decay faster than exponentially, Green functions and vertices admit sparse expansions with only $\Xi^{\mathrm{BSE}}$1 basis functions. The resulting BSE solver has $\Xi^{\mathrm{BSE}}$2 time and $\Xi^{\mathrm{BSE}}$3 memory complexity, with exponential convergence demonstrated for the Hubbard atom, a multi-orbital weak-coupling model, and a realistic impurity problem (Wallerberger et al., 2020).

Quantum-computing-oriented reformulations have also appeared. The homogeneous BSE in ladder approximation can be discretized into a generalized eigenvalue problem with one symmetric and one non-symmetric matrix, transformed into a quadratic unconstrained binary optimization problem, and then solved with a quantum annealer. Numerical studies up to matrix dimension $\Xi^{\mathrm{BSE}}$4 show good agreement with classical solvers and indicate nontrivial scalability features (Fornetti et al., 2024).

At the formal edge of the subject, the BSE language has been extended to classical gravitational bound states. In that construction, the space of classical amplitudes is defined by quotienting out symmetrization over internal graviton exchanges, the classical kernel exponentiates in impact-parameter space, and the resulting analytic structure makes explicit the continuation between scattering and bound observables, including a relativistic analogue of Sommerfeld enhancement (Adamo et al., 2022). This suggests that the BSE has become less a single model-specific equation than a general organizing principle for two-body dynamics, kernel resummation, and spectral extraction across relativistic and many-body theory.