Bethe-Salpeter Equation: An Overview

Updated 23 November 2025

Bethe-Salpeter Equation is a many-body perturbation theory framework that computes neutral excitations and optical spectra by accurately modeling electron–hole interactions.
It is integrated with DFT and GW methodologies in practical workflows, enabling precise prediction of excitonic and charge-transfer states in various materials.
Advances such as finite-field evaluations and low-rank approximations have significantly reduced computational costs while extending BSE applications to complex spectral analyses.

The Bethe-Salpeter Equation (BSE) is a cornerstone of many-body perturbation theory, providing a rigorous framework for computing neutral excitation energies and optical spectra in molecules and solids. The BSE describes the propagation of two-particle (electron–hole) correlations and explicitly accounts for electron–hole attraction and exchange, enabling the quantitative prediction of excitons, charge-transfer excitations, and core-level spectra in broad classes of materials.

1. Formalism and Physical Interpretation

The BSE emerges from the equation of motion for the two-particle Green's function $L(1,2;1',2')$ , capturing the correlated response of a many-electron system to external perturbations. The BSE for $L$ in operator notation is

$L = L_0 + L_0\,K\,L$

where $L_0$ is the independent-particle propagator and $K$ is the electron–hole interaction kernel, typically derived from functional differentiation of the many-body self-energy with respect to the single-particle Green's function (Blase et al., 2020). In practice, one projects the BSE onto a product basis of valence $\psi_v$ and conduction $\psi_c$ single-particle orbitals, resulting in an effective Schrödinger-like eigenvalue problem for neutral excitations: $\sum_{v'c'} \left[(\epsilon_c - \epsilon_v)\,\delta_{vv'}\delta_{cc'} + K_{vc, v'c'}\right] A^S_{v'c'} = \Omega_S A^S_{vc}$ where $\epsilon_v$ , $\epsilon_c$ are quasiparticle energies, $A^S_{vc}$ are the excitonic amplitudes, and $\Omega_S$ are excitation energies (Nguyen et al., 2019, Blase et al., 2020).

The kernel $K$ contains a repulsive bare-exchange term and an attractive direct term mediated by the statically screened Coulomb interaction $W$ , i.e., $K = K^d + K^x$ , with

$K^d_{vc,v'c'} = \iint \psi_v(r) \psi_c^*(r) W(r, r') \psi_{v'}^*(r') \psi_{c'}(r')\,dr dr'$

$K^x_{vc,v'c'} = -\iint \psi_v(r) \psi_{c'}(r) v_c(r, r') \psi_{v'}^*(r') \psi_{c}^*(r')\,dr dr'$

where $v_c(r, r')=1/|r - r'|$ and $W(r, r') = \epsilon^{-1}(r,r') v_c(r,r')$ (Nguyen et al., 2019, Rebolini et al., 2013).

2. BSE in Practical Many-Body Workflows

Realistic electronic-structure calculations typically combine:

(i) DFT ground-state: Providing orbitals and energies for the Kohn-Sham reference.
(ii) GW correction: Yielding quasiparticle energies and the statically screened Coulomb interaction $W$ (via the RPA dielectric response) (Blase et al., 2020, Gilmore et al., 2016).
(iii) BSE solution: Constructing and diagonalizing the excitonic Hamiltonian for neutral excitations, optionally within the Tamm-Dancoff approximation (TDA) (neglecting antiresonant couplings) (Blase et al., 2020, Unzog et al., 2022).

The standard workflow is:

Compute ground-state DFT orbitals.
Compute RPA screening and perform a single-shot $G_0W_0$ calculation for quasiparticle corrections.
Assemble the BSE kernel $K_{vc,v'c'}$ , including both exchange and screened-direct terms.
Solve the BSE eigenproblem (direct diagonalization, Davidson/Lanczos iterative solvers, or Haydock recursion for spectral properties).

This framework yields accurate optical gaps, absorption spectra, and can describe both localized (Frenkel) and delocalized (Wannier, charge-transfer) excitons in molecules and solids (Blase et al., 2020, Liu et al., 2019). For core-level spectroscopy (XAS, XES, RIXS), the BSE is constructed in a product basis of core and conduction states, with the same basic kernel structure (Gilmore et al., 2016, Vinson et al., 2010, Unzog et al., 2022).

3. Algorithmic Strategies and Computational Scaling

Early BSE implementations incurred steep operator scaling: full matrix build and diagonalization cost scales as $O(N^6)$ in the basis set size $N$ . Multiple algorithmic advances have reduced this bottleneck:

Finite-field evaluation: Bypasses explicit dielectric matrix inversion, obtaining $W$ via finite-field density responses in localized orbital bases, reducing required integrals and scaling toward $O(N^3)$ for large $N$ (Nguyen et al., 2019).
Low-rank/decomposition approaches: Represent BSE matrices as diagonal plus low-rank factors, enabling iterative solvers and Galerkin projections for rapid convergence in an auxiliary reduced basis (scaling typically $O(N^3)$ for low-lying excitations) (Benner et al., 2015).
Density-matrix perturbation: Employs Liouville techniques and projective dielectric eigenpotentials (PDEP), avoiding explicit empty-state summations and allowing GPU acceleration to systems with thousands of atoms (Yu et al., 23 Sep 2024).
Real-time and stochastic BSE: Recasting the BSE as a time-dependent propagation, with stochastic orbitals replacing deterministic sums, leads to quadratic ( $O(N^2)$ ) scaling for extended systems of thousands of electrons (Rabani et al., 2015).
Energy-specific iterative solvers: Divide the spectrum into small energy windows, targeting only desired excitations using energy-filtered trial spaces and orthogonalization, achieving efficient access to high-lying core-valence spectra (Hillenbrand et al., 31 Oct 2024).

These developments enable calculations on state-of-the-art molecular, condensed, and defect-containing supercells on modern workstation and HPC platforms (Yu et al., 23 Sep 2024, Nguyen et al., 2019, Benner et al., 2015, Hillenbrand et al., 31 Oct 2024).

4. Extensions, Dynamical Effects, and Kernel Approximations

The conventional BSE employs static (adiabatic) screening in the kernel and TDA for practical tractability (Blase et al., 2020). Recent advancements include:

Dynamical BSE: Incorporating frequency-dependence of $W$ , either by first-order perturbative corrections or exact reformulation as a frequency-independent eigenproblem in singles + doubles space (Bintrim et al., 2021, Loos et al., 2020). This enables access to genuinely doubly excited states, previously inaccessible in static BSE.
Beyond the GW kernel: Alternative kernels such as the T-matrix and hybrid bubble + ladder diagrams have been studied to address ground-state instabilities, underestimation of triplet energies, and lack of double excitations (Orlando et al., 2023).
BSE+RPA hybrid schemes: BSE+ combines BSE for low-energy transitions with RPA inclusion of high-energy screening, dramatically accelerating convergence of dielectric properties and electron-energy-loss spectra for extended materials (Søndersted et al., 2023).
Parameter-free approaches: Use of Koopmans-compliant functionals for the underlying single-particle basis removes the need for GW, with direct minimization for W in a maximally-localized Wannier basis (Elliott et al., 2019).

Table: Key Features and Computational Strategies

Approach	Scaling	Notable Features
Standard GW+BSE	$O(N^6)$	Direct diagonalization; high cost
Low-rank/Reduced Basis (Benner et al., 2015)	$O(N^3)$ – $O(N^4)$	Low-rank factorization; Galerkin projection
Finite field (Nguyen et al., 2019)	$O(N^4)$ – $O(N^3)$	Finite-difference screening; localization
Real-time/Stochastic (Rabani et al., 2015)	$O(N^2)$	Time-dependent stochastic sampling
Density-matrix, PDEP (Yu et al., 23 Sep 2024)	$O(N^4)$	GPU, no explicit dielectric matrices
Energy-specific (Hillenbrand et al., 31 Oct 2024)	$O(N^4)$	Sliding-window iterative solutions
Dynamical BSE (Bintrim et al., 2021)	$O(N^5)$ (DF)	Frequency-dependent kernel; doubles

5. Applications: Spectra, Excitons, and Core-Level Physics

The BSE formalism provides quantitatively accurate predictions for:

Singlet and triplet optical excitations: Typical errors for low-lying singlet excitations are 0.1–0.3 eV for molecules, matching or exceeding the accuracy of leading TD-DFT hybrid functionals (Blase et al., 2020, Liu et al., 2019).
Exciton binding energies and spectra: Both local and charge-transfer excitons in molecular aggregates and solids; directly yields oscillator strengths and absorption spectra (Nguyen et al., 2019).
Core-level spectroscopy: XAS, XES, and (R)IXS spectra for K- and L-edges—fine structure and exciton binding in agreement with experiment. BSE approaches surpass supercell core-hole techniques for resolution of fine features (Unzog et al., 2022, Gilmore et al., 2016, Vinson et al., 2010).
Pump–probe spectroscopy: BSE applied to both valence and core excitons enables the ab initio description of transient x-ray absorption, revealing signatures of photoexcited many-body states (Farahani et al., 8 Feb 2024).
Defects and nanostructures: GPU and density-matrix accelerated BSE calculations allow handling supercells with up to several thousand atoms, capturing finite-size convergence and defect physics (Yu et al., 23 Sep 2024).
Method benchmarks: Systematic comparisons to high-level quantum chemistry for small molecules (e.g., Thiel’s set, C60) show consistent agreement, with detailed basis set and GW-self-energy convergence analyses (Liu et al., 2019, Elliott et al., 2019).

6. Limitations and Directions in BSE Development

Despite its wide applicability, BSE faces known challenges:

Missing double excitations: Standard (static) BSE cannot represent states with dominant double-excitation character without dynamical kernel extensions (Bintrim et al., 2021, Loos et al., 2020, Orlando et al., 2023).
Triplet instabilities: Underestimation and possible collapse (complex frequencies) of triplet gaps in strongly correlated or bond-breaking regimes; only partially alleviated by TDA (Rebolini et al., 2013, Blase et al., 2020, Orlando et al., 2023).
Dependence on the GW starting point: Poor treatment of systems with problematic underlying GW self-energies; influences overall spectral accuracy (Blase et al., 2020, Liu et al., 2019).
Lack of analytic gradients: Routine analytic excited-state gradients are not yet broadly available, limiting geometry optimization for excited states (Blase et al., 2020).
Static approximation for W: Adiabatic treatment excludes dynamical kernel effects, with consequences for Rydberg and double excitations.

Research directions include:

Efficient O(N^2–N⁴⁾ scaling algorithms for massive systems (Yu et al., 23 Sep 2024, Hillenbrand et al., 31 Oct 2024, Nguyen et al., 2019).
Systematic inclusion of dynamical screening and higher-order kernels for correlated materials (Bintrim et al., 2021, Loos et al., 2020, Orlando et al., 2023).
Hybrid BSE+RPA schemes for faster convergence of dielectric properties (Søndersted et al., 2023).
Parameter-reduction strategies, e.g. direct BSE on Koopmans-corrected Hamiltonians (Elliott et al., 2019).
Machine-learning quasiparticle energies to reduce GW cost (Hillenbrand et al., 31 Oct 2024).

7. Summary and Outlook

The Bethe-Salpeter Equation constitutes the state-of-the-art for the ab initio calculation of neutral excitation spectra—capturing many-body effects such as excitons, screening, and multiplet structure with systematic improvability. Modern algorithmic developments allow application to ever-larger systems, while ongoing theory work seeks to address dynamical, multireference, and excited-state nuclear dynamics. Its predictive power in molecules, nanostructures, and solids underlies contemporary advances in computational spectroscopy and materials design (Blase et al., 2020, Nguyen et al., 2019, Søndersted et al., 2023, Yu et al., 23 Sep 2024, Bintrim et al., 2021).