Linear-Scaling BSE Framework

• The linear-scaling Bethe–Salpeter framework is a collection of methods that reformulate the BSE to achieve near-linear scaling by leveraging localization and sparse representations.
• It employs matrix-free linear response, low-rank dielectric representations, and truncation of orbital interactions to efficiently model optical spectra and excitonic effects.
• These techniques enable scalable simulation of excitations in complex materials, facilitating studies on systems with thousands of atoms and heterogeneous properties.

Searching arXiv for the cited papers and closely related linear-scaling Bethe–Salpeter work. Searching arXiv for "GPU-Accelerated Solution of the Bethe-Salpeter Equation for Large and Heterogeneous Systems". A linear-scaling Bethe–Salpeter framework is a class of Bethe–Salpeter equation (BSE) methodologies in which the dominant operations for neutral excitations, optical spectra, or screened electron–hole interactions are reformulated so that they scale as $O(N)$ or near-linearly with system size for fixed accuracy in appropriate regimes. Across recent work, this objective is realized by different but related strategies: matrix-free linear-response formulations that avoid explicit dielectric matrices and sums over virtual states, low-rank representations of dielectric screening, localization and truncation of occupied or electron–hole degrees of freedom, sparse real-space excitonic Hamiltonians, and iterative spectral methods that avoid dense diagonalization (Yu et al., 2024, Galvani et al., 29 Jun 2026, Nguyen et al., 2019). The term therefore does not denote a single algorithm, but a family of BSE implementations tailored to exploit locality, nearsightedness, or compact excitons in gapped and often heterogeneous materials.

1. Scope and defining variants

Three representative realizations delineate the present meaning of a linear-scaling BSE framework.

Framework	Core representation	Scaling statement
WEST-BSE	Liouville/DMPT BSE with PDEP and nearsightedness-based localization/truncation	BSE-specific stages scale near-linearly in insulators and semiconductors
Real-space hBN framework	Sparse localized pair basis with KPM	Absorption spectra with $O(N)$ cost
FF-BSE	Finite-field screening with localized orbitals	A scalable pathway to linear or near-linear BSE in practice

In the WEST-BSE implementation, neutral excitations are solved within the linearized Liouville formalism and density matrix perturbation theory (DMPT) under the Tamm–Dancoff approximation (TDA). The near-linear behavior arises specifically in the BSE-specific stages, notably the construction of screened Coulomb integrals and the application of the Liouville operator, by combining DMPT, a low-rank representation of dielectric screening via projective dielectric eigenpotentials (PDEP), and nearsightedness-based localization and truncation (Yu et al., 2024).

The real-space framework introduced for disordered boron-nitride-derived systems starts from a localized tight-binding description, performs a sublattice-resolved perturbative decoupling, maps localized electron–hole pairs onto a sparse excitonic Hamiltonian, and evaluates spectra with the Kernel Polynomial Method (KPM). Its linear scaling is tied to sparse matrix–vector products at fixed spectral resolution and to a physically motivated cutoff on electron–hole separation (Galvani et al., 29 Jun 2026).

The finite-field approach replaces dielectric-matrix construction and virtual-state summations by self-consistent finite-field calculations that directly yield screening on localized pair perturbations. In combination with recursively bisected localized orbitals, it produces a sparse direct kernel and a workflow described as a scalable pathway to linear or near-linear BSE in practice, especially when paired with linear-scaling DFT (Nguyen et al., 2019).

2. Formal BSE structures and operator formulations

The common theoretical substrate is the excitonic eigenproblem, but the operational form differs substantially among frameworks. In the resonant TDA basis $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$, the usual BSE reads

$$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$

with

$$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$

The direct kernel contains the statically screened Coulomb interaction $W=\epsilon^{-1}v_c$, while the exchange kernel involves the bare Coulomb interaction (Yu et al., 2024).

WEST-BSE rewrites this problem in orbital space. Under the Liouville/DMPT formulation and TDA, the $s$-th excitation satisfies

$(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$

where $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$ defines the linear change of the density matrix

$$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$

This implicit application of the kernels avoids explicit electron–hole basis construction and operates through projectors in occupied and conduction subspaces (Yu et al., 2024).

The finite-field formulation also adopts a Liouville viewpoint, but frames the response problem as a linear system for the density matrix under a monochromatic perturbation,

$O(N)$0

with $O(N)$1. In batch form, this yields non-homogeneous linear systems for occupied-state response orbitals $O(N)$2, with the screened integrals $O(N)$3 entering through the direct kernel. The same work also states the full BSE in block form with $O(N)$4 and $O(N)$5 matrices, so the framework is not intrinsically restricted to TDA, even though the Hermitian TDA form is typically used for computational convenience (Nguyen et al., 2019).

The real-space hBN formulation starts from the TDA BSE but translates it into a localized pair basis

$O(N)$6

where $O(N)$7 and $O(N)$8 denote hole and electron positions on distinct sublattices. The excitonic Hamiltonian is written as

$O(N)$9

with $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$0 neglected in the present work. The direct interaction is approximated as diagonal in the localized pair basis,

$$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$1

and direct optical absorption probes the $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$2 block (Galvani et al., 29 Jun 2026).

3. Mechanisms that produce linear or near-linear scaling

The essential reduction in complexity is achieved by eliminating dense global objects and retaining only physically relevant local couplings. In WEST-BSE, DMPT removes explicit sums over empty and occupied states in the response construction, while PDEP replaces the full dielectric matrix by a low-rank spectral representation

$$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$3

with rank $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$4 chosen to converge screening and in practice equal to “a few times” the number of electrons, significantly smaller than the full plane-wave dimension $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$5. This reduces the cost of applying $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$6 in $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$7 to $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$8 per vector instead of $$\{|vc\rangle = |\psi_v\rangle \otimes |\psi_c\rangle\}$$9 (Yu et al., 2024).

The decisive near-linear effect in WEST-BSE comes from nearsightedness. For gapped systems, the one-particle density matrix decays exponentially,

$$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$0

Occupied orbitals are localized by the JADE algorithm, and screened Coulomb integrals are truncated using the normalized overlap

$$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$1

If $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$2, then $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$3. At fixed density and fixed $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$4, each localized orbital overlaps only with a constant number $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$5 of neighbors, so the number of $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$6 integrals scales as $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$7 rather than $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$8. This is identified explicitly as the primary source of near-linear scaling in the BSE-specific stages (Yu et al., 2024).

The finite-field framework obtains sparsity differently but with the same objective. Occupied Bloch states are transformed to a localized basis by recursive bisection. Only pairs of localized orbitals with significant spatial overlap, or within a cutoff radius, are retained. The retained screened-integral matrix is sparse, the number of non-zero pairs scales linearly with system size in insulating systems with fixed accuracy, and both storage and kernel application become $$\sum_{v'c'} H^{BSE}_{vc,v'c'} A_{v'c'} = \Omega A_{vc},$$9 per iteration. The screening stage itself remains limited by the cost of the underlying self-consistent finite-field solves unless linear-scaling DFT is also used (Nguyen et al., 2019).

In the real-space hBN framework, the electron–hole pair basis is truncated by a separation cutoff $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$0. For excitons with Bohr radius $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$1, the approximation is controlled by convergence in $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$2. The excitonic Hamiltonian dimension then scales as

$$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$3

and the number of nonzeros scales linearly with that dimension because the effective electron and hole Hamiltonians are short-ranged and $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$4 is local in the retained basis. KPM then evaluates absorption and density of states through sparse matrix–vector products, so for fixed spectral resolution the computational cost scales effectively as $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$5 (Galvani et al., 29 Jun 2026).

4. Algorithmic realizations and computational workflows

The WEST-BSE workflow begins from plane-wave DFT eigenvalues and eigenvectors, with KS orbitals imported from Quantum ESPRESSO and ground-state calculations performed with SG15 ONCV pseudopotentials, PBE, a 60 Ry cutoff, and $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$6-only supercells. Quasiparticle corrections are introduced either through full-frequency $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$7 or, for $$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$8–$$H^{BSE}_{vc,v'c'} = (\epsilon_c-\epsilon_v)\delta_{vv'}\delta_{cc'} + K^{dir}_{vc,v'c'} - K^{exc}_{vc,v'c'}.$$9 bonded systems, by a scissor operator. The dielectric response is constructed in the PDEP basis by iterative DMPT application of $W=\epsilon^{-1}v_c$0 through Sternheimer-type equations in the projected conduction space, thereby avoiding both virtual orbitals and dielectric-matrix inversion. The production implementation uses TDA for both spin-conserving and spin-flip BSE, Davidson for low-lying vertical excitation energies, and Liouville–Lanczos recursion for spectra (Yu et al., 2024).

Its GPU implementation is organized hierarchically through images, pools, and band groups. Excitonic vectors are block-distributed across band groups; transformations between KS and localized bases require MPI_Allgather across band groups, and non-blocking variants are overlapped with the evaluation of $W=\epsilon^{-1}v_c$1 and $W=\epsilon^{-1}v_c$2, with completion enforced before $W=\epsilon^{-1}v_c$3. The software stack uses cuFFT and cuBLAS for FFT and BLAS kernels, OpenACC for bespoke loops, and GPU-aware MPI for direct device communications. Calculations are reported on up to 4096 GPUs, with bottlenecks at very high GPU counts attributed to MPI communications and FFT collective patterns (Yu et al., 2024).

The finite-field workflow starts from a semilocal or hybrid DFT reference, optionally corrected at the eigenvalue level by $W=\epsilon^{-1}v_c$4. For each non-negligible localized pair, one first forms the bare pair potential

$W=\epsilon^{-1}v_c$5

then solves two independent self-consistent problems with $W=\epsilon^{-1}v_c$6, obtains the perturbed densities $W=\epsilon^{-1}v_c$7, forms the central finite-difference response

$W=\epsilon^{-1}v_c$8

and finally reconstructs the screened integral

$W=\epsilon^{-1}v_c$9

This yields matrix-free access to screening without explicit dielectric matrices or unoccupied-state sums (Nguyen et al., 2019).

The hBN framework begins from a localized $s$0-orbital tight-binding Hamiltonian and performs a second-order perturbative decoupling into effective electron and hole Hamiltonians on the B and N sublattices,

$s$1

In disordered systems, diagonal Anderson disorder contributes a pair-space term $s$2. After restricting to the $s$3 block and truncating pair states with $s$4, absorption is computed from a dipole-prepared vector by Chebyshev recursion,

$s$5

followed by KPM reconstruction with a damping kernel. The implementation uses KWANT for KPM and ARPACK for selected eigenstates (Galvani et al., 29 Jun 2026).

5. Validation, demonstrated scales, and physical results

These frameworks have been validated in markedly different physical regimes. WEST-BSE was demonstrated for spin defects in wide-band-gap materials and for a large heterogeneous defect complex. For $s$6 in diamond, GW-BSE@PBE yielded vertical excitation energies of $s$7 eV for $s$8, $s$9 eV for $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$0, and $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$1 eV for $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$2 in a $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$3 supercell, with agreement with previous GW-BSE within $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$4 eV. For $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$5 in diamond, singlets showed weak supercell-size dependence, while triplets exhibited stronger finite-size effects due to defect–bulk mixing. For $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$6 in 3C-SiC, the dilute-limit vertical excitation energy was approximately $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$7 eV and the lowest-state radiative lifetime approximately $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$8 ns, close to experiment at $(D + K^{1e} - K^{1d}) A_s = \omega_s A_s,$9 ns. The same implementation was then applied to a diamond lattice with 1727 atoms, where an $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$0 center at a dislocation core showed a $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$1 splitting greater than $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$2 eV and vertical excitation energies of $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$3 and $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$4 eV (Yu et al., 2024).

The reported performance metrics underscore that the scaling claim is stage-specific. For diamond $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$5 with 999 atoms and 3998 electrons, $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$6 and $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$7 scale to 512 nodes (2048 GPUs), while $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$8 and $A_s=\{|a_{s,v}\rangle : v=1,\dots,N_{occ}\}$9 are two orders faster and scale to 64 nodes (256 GPUs). For bulk Si with 1024 atoms and 4096 electrons, $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$0 scales efficiently to 1024 nodes (4096 GPUs), with a Lanczos chain length $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$1. JADE localization for Si-512 with 2048 occupied states requires approximately 3 minutes on a single A100, and the $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$2 and $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$3 stages complete under approximately 400 s on 64 nodes (256 GPUs) for $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$4 with 999 atoms (Yu et al., 2024).

The real-space hBN framework was applied to Anderson-disordered monolayer hexagonal boron nitride with up to $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$5 orbitals, with disorder statistics on approximately $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$6-atom cells using about 200 configurations per $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$7. It revealed disorder-induced asymmetric broadening of bright excitons, activation of dark spectral weight under broken translational symmetry, and a crossover in the redshift of the main absorption peak from $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$8 at small disorder to $$\Delta \rho_s = \sum_v |a_{s,v}\rangle\langle \psi_v|.$$9 at larger disorder, beyond approximately $O(N)$00 eV. The excess FWHM beyond the kernel broadening $O(N)$01 meV was fit by $O(N)$02 with $O(N)$03 up to $O(N)$04 eV. Real-space analysis showed Anderson localization of the exciton center of mass while the relative coordinate remained compact; the pristine binding energy was approximately $O(N)$05 eV, calculations used $O(N)$06 Å, and the compact exciton radius was reported as $O(N)$07 Å (Galvani et al., 29 Jun 2026).

The finite-field framework was benchmarked on molecules and condensed phases. For the Thiel set of 28 molecules, lowest singlet excitation energies agreed closely with high-level quantum chemistry best estimates, and using PBE0 as the starting point improved accuracy relative to $O(N)$08. For $O(N)$09, 120 doubly occupied valence states imply 14,400 pairs in principle, yet only 4,284 screened integrals were required with negligible change in the absorption spectrum. The method was also applied to liquid water with 64–256 molecules and proton-disordered hexagonal ice with 96 molecules. Reported exciton binding energies were $O(N)$10 eV and $O(N)$11 eV for water, and $O(N)$12 eV and $O(N)$13 eV for ice, with size effects from 64 to 256 molecules described as modest, around $O(N)$14–$O(N)$15 eV (Nguyen et al., 2019).

6. Applicability, limitations, and common misconceptions

A persistent misconception is that “linear-scaling BSE” implies that every stage of an end-to-end ab initio workflow is strictly linear. The evidence is more specific. In WEST-BSE, the BSE-specific stages $O(N)$16 and $O(N)$17 exhibit near-linear scaling for insulators under fixed $O(N)$18 and density, but the overall GW-BSE workflow remains dominated at large $O(N)$19 by dielectric construction through PDEP and, when used, full-frequency GW, both of which are superlinear. In the finite-field framework, sparse kernel application is near $O(N)$20, but with conventional plane-wave SCF the formal total work for screening remains $O(N)$21–$O(N)$22, becoming genuinely near-linear only when the finite-field SCF is combined with linear-scaling DFT. In the hBN real-space framework, the $O(N)$23 claim refers to a sparse model with fixed $O(N)$24 and fixed spectral resolution, rather than to conventional dense ab initio GW+BSE (Yu et al., 2024, Nguyen et al., 2019, Galvani et al., 29 Jun 2026).

All three variants rely on locality assumptions that can weaken outside gapped, compact-exciton regimes. WEST-BSE states explicitly that near-linear scaling hinges on nearsightedness and degrades in metals, very small-gap systems, or highly heterogeneous interfaces where the localization length $O(N)$25 grows and the number of overlapping neighbors $O(N)$26 increases. The finite-field approach likewise depends on exponential localization of occupied states and screened interactions; in metals, small-gap systems, or systems with very large exciton radii, more pairs survive the cutoff and the near-linear behavior deteriorates. The hBN framework requires sufficiently large sublattice asymmetry $O(N)$27, moderate inter-sublattice couplings, compact excitons with $O(N)$28, and short-ranged effective hopping; strong inter-sublattice hybridization, weak binding, or very long-range interactions can degrade accuracy (Yu et al., 2024, Nguyen et al., 2019, Galvani et al., 29 Jun 2026).

The approximations also differ materially. WEST-BSE currently uses TDA in production, with full BSE beyond TDA identified as future work, and spin–orbit coupling listed among planned extensions. The hBN framework neglects exchange, uses static Rytova–Keldysh screening with fixed $O(N)$29 Å, and omits dynamical screening and phonon coupling. The finite-field approach uses static screening and typically the Hermitian TDA form, although it can include exchange–correlation kernel effects beyond RPA through the self-consistent finite-field response when hybrid functionals are used (Yu et al., 2024, Galvani et al., 29 Jun 2026, Nguyen et al., 2019).

Taken together, these works indicate that a linear-scaling Bethe–Salpeter framework is best understood as a set of physically structured reductions of the BSE problem. In one direction, low-rank dielectric eigenpotentials, localization, and GPU-distributed Liouville solvers make first-principles GW-BSE tractable for supercells of 1000–1700 atoms. In another, sparse real-space excitonic Hamiltonians and KPM extend excitonic spectroscopy to disordered systems with up to $O(N)$30 orbitals. In a third, finite-field screening and localized orbitals remove dielectric matrices and empty-state sums and provide a route to sparse, scalable BSE kernels. A plausible implication is that future linear-scaling BSE developments will continue to be domain-specific, with the form of locality—density-matrix nearsightedness, pair-space truncatability, or sparse real-space screening—determining both the attainable scaling and the class of materials for which it is reliable.