The Bethe–Salpeter equation

The Bethe–Salpeter equation (BSE) describes excitons, the bound electron–hole pairs of semiconductors and insulators. Xatu solves it in a basis of localized orbitals, with a static screened interaction.

Excitons as electron–hole pairs

An exciton of centre-of-mass momentum \(\bm{Q}\) is written as a superposition of electron–hole pairs, with an electron in conduction band \(c\) at \(\bm{k}+\bm{Q}\) and a hole in valence band \(v\) at \(\bm{k}\):

\[|X\rangle = \sum_{v,c,\bm{k}} A^{\bm{Q}}_{vc}(\bm{k})\; c^\dagger_{c,\bm{k}+\bm{Q}}\, c_{v,\bm{k}}\, |GS\rangle .\]

Projecting the interacting Hamiltonian onto these pairs gives an eigenvalue problem for the coefficients \(A^{\bm{Q}}_{vc}(\bm{k})\):

\[\sum_{v',c',\bm{k}'} H_{vc,v'c'}(\bm{k},\bm{k}',\bm{Q})\, A^{\bm{Q}}_{v'c'}(\bm{k}') = E_X\, A^{\bm{Q}}_{vc}(\bm{k}) .\]

Working in the band basis of a mean-field Hamiltonian (Hartree–Fock, DFT or tight-binding) absorbs the self-energy corrections into the quasiparticle energies. The equation Xatu solves is then

\[\left( \varepsilon_{c,\bm{k}+\bm{Q}} - \varepsilon_{v,\bm{k}} \right) A^{\bm{Q}}_{vc}(\bm{k}) + \sum_{v',c',\bm{k}'} K_{vc,v'c'}(\bm{k}, \bm{k}', \bm{Q})\, A^{\bm{Q}}_{v'c'}(\bm{k}') = E_X\, A^{\bm{Q}}_{vc}(\bm{k}),\]

where

  • \(\varepsilon_{n,\bm{k}}\) are the quasiparticle band energies (plus the optional # scissor);

  • \(A^{\bm{Q}}_{vc}(\bm{k})\) are the exciton coefficients, written to .states — exciton eigenstates;

  • \(K = -(D - X)\) is the interaction kernel, with the direct term \(D\) (screened electron–hole attraction) and the optional exchange term \(X\) (# exchange).

This is the Tamm–Dancoff approximation: coupling between resonant and anti-resonant pairs is neglected.

Self-energy correction

With # selfenergy true, the band energies entering the BSE are corrected by a self-energy computed from the same model interaction (# selfenergy.potential):

\[\varepsilon_{n,\bm{k}} \;\to\; \varepsilon_{n,\bm{k}} + \Sigma_n(\bm{k}) .\]

\(\Sigma_n(\bm{k})\) is a Hartree–Fock-type term: a direct (Hartree) contribution minus an exchange contribution, summed over the valence bands of the window and the k-mesh. It follows Eqs. (2.10)–(2.12) of arXiv:2510.25009. It is available for the real-space method, and can be written to .selfenergy — band self-energy with -i.

The size of the problem, the BSE dimension, is \(N_v \times N_c \times N_{\bm{k}}\): bands from # bands/# bandlist, and k-points \(N_{\bm{k}} =\) ncells\(^d\).

Interaction matrix elements

The kernel is built from the electron–hole interaction potential chosen with # potential (see Interaction and screening). The matrix elements can be computed in two ways:

Real space (default)

The potential is evaluated between orbitals of the lattice, treated as point-like charges at their centres. The divergence at \(r=0\) is removed with # regularization.

Reciprocal space

The potential is summed over reciprocal-lattice vectors \(\bm{G}\) with \(|\bm{G}|\) below # gcutoff, which enables this method (see Exciton file).

Screened RPA potential

With # potential rpa, the direct term uses the numerical RPA screening. In the strictly 2D approach, the screened Coulomb potential is a matrix at each \(\bm{q}\) of the Brillouin zone:

\[W_{\bm{G}\bm{G}'}(\bm{q}) = \sqrt{v_c(\bm{q}+\bm{G})}\; \epsilon^{-1}_{\bm{G}\bm{G}'}(\bm{q})\; \sqrt{v_c(\bm{q}+\bm{G}')} ,\]

where \(\epsilon^{-1}_{\bm{G}\bm{G}'}(\bm{q})\) is the inverse RPA dielectric matrix of Interaction and screening.

If a thickness \(d_\perp\) is given in the screening file, the quasi-2D (Q2D) mode is used instead:

\[\bar{W}_{\bm{G}\bm{G}'}(\bm{q}) = \sqrt{\bar{v}_c(\bm{q}+\bm{G})}\; \bar{\epsilon}^{-1}_{\bm{G}\bm{G}'}(\bm{q})\; \sqrt{\bar{v}_c(\bm{q}+\bm{G}')} ,\]

with the Coulomb potential averaged over the thickness of the layer,

\[\bar{v}_c(\bm{q}) = \int_{-d_{\perp}/2}^{d_{\perp}/2} \int_{-d_{\perp}/2}^{d_{\perp}/2} v_c(\bm{q}, z-z')\, \mathrm{d}z\, \mathrm{d}z' .\]

Here \(v_c(\bm{q}, z-z')\) is the Coulomb potential in the mixed \((\bm{q}, z)\) representation: the in-plane Fourier transform of \(V(\bm{r}-\bm{r}')\), keeping \(z\) and \(z'\).

Solving the BSE

The BSE matrix is Hermitian. It is built and diagonalized with Armadillo; the solver is chosen with -m:

-m

Method

diag

Full diagonalization (default). All states; memory and time grow fast with the BSE dimension.

davidson

Iterative Davidson solver for the lowest -n states.

sparse

Lanczos-based sparse diagonalization (ARPACK) for the lowest -n states.

From the solution, Xatu derives the energies, wavefunctions in real and reciprocal space, spin, and optical matrix elements (see Output files).