🔬 materials science

Analytical Forces from the Bethe-Salpeter Equation for Large-Scale Excited-State Relaxation

This paper presents an efficient, GPU-accelerated plane-wave implementation of analytical nuclear forces within the Bethe-Salpeter equation framework, enabling scalable excited-state relaxation studies for large solid-state systems and demonstrating its ability to correct semilocal TDDFT errors in defect environments like hexagonal boron nitride.

Original authors: Yu Jin, Victor Wen-zhe Yu, Marco Govoni, Giulia Galli

Published 2026-07-24

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Yu Jin, Victor Wen-zhe Yu, Marco Govoni, Giulia Galli

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the world of materials science as a giant, bustling dance floor. In this dance, electrons are the dancers, and the atoms they orbit are the partners holding their hands. Usually, we study how these dancers move when the music is calm and steady (the ground state). But sometimes, a flash of light hits the floor, and an electron gets a sudden burst of energy, jumping to a new, excited rhythm. This is an "excited state."

When an electron gets excited, it doesn't just sit still; it pulls its atomic partners along, causing the whole structure to wobble, stretch, or twist. This dance between the excited electron and the moving atoms is called "vibronic coupling." It's the reason why some materials glow with specific colors, why certain defects in diamonds can act as quantum sensors, and why solar cells might lose energy as heat. To understand these phenomena, scientists need to calculate the "forces" pushing and pulling on the atoms in these excited states. It's like trying to predict exactly how a trampoline will bounce if you jump on it while holding a heavy backpack, but the backpack is made of invisible, jittery energy.

For a long time, calculating these forces for complex materials was like trying to solve a million-piece puzzle while blindfolded. The math was so heavy and slow that scientists could only do it for tiny molecules. They had to guess or use shortcuts for bigger things like crystals or 2D materials, often missing the subtle ways the material's environment changes the dance. This is where the story of this new paper begins: a team of researchers has built a super-fast, super-accurate way to calculate these forces for huge systems, finally letting us watch the full dance of excited electrons in real, complex materials.

The Big Leap: Calculating the "Push" of Excited Electrons

The researchers, led by Yu Jin and Giulia Galli, have developed a new computer method to calculate analytical nuclear forces for excited states using a powerful mathematical tool called the Bethe-Salpeter Equation (BSE). Think of the BSE as a high-definition camera that captures the complex relationship between an excited electron and the "hole" it leaves behind (a positive charge). While previous methods could take pictures of the energy levels, they struggled to calculate the forces that move the atoms, especially in large systems with hundreds of atoms.

The team's breakthrough is a new "recipe" that combines two clever tricks: density-matrix perturbation theory and a Lagrangian approach (specifically a "Z-vector" method). In the old way of doing things, to find out how the atoms move when an electron gets excited, scientists had to solve a massive, separate math problem for every single direction every atom could move. If you had a crystal with 500 atoms, that meant solving 1,500 separate, incredibly expensive problems. It was like trying to measure the wind by asking every single leaf on a tree to report its own direction individually.

The new method is like hiring a single, super-smart wind detective who can figure out the wind's effect on the whole tree by solving just one master equation. By avoiding the need to sum up empty energy states and by using a "Z-vector" to handle the response of the electrons, the team made the calculation speed up dramatically. They also added GPU acceleration (using powerful graphics cards usually found in gaming computers) to handle the heavy lifting. The result? They can now simulate excited-state relaxations in solid materials containing hundreds of atoms, something that was previously prohibitively expensive and slow.

Testing the Dance: Two Different Defects

To prove their new method works, the team tested it on two very different "dance floors" (defects in materials) to see how the excited electrons behaved.

1. The Diamond Diamond (NV Center)
First, they looked at a nitrogen-vacancy (NV) center in diamond. Imagine a diamond lattice where a carbon atom is missing and replaced by a nitrogen atom. This is a famous defect used in quantum technology.

The Finding: In this relatively uniform environment, the new BSE method and the older, simpler method (called TDDFT) gave very similar results. Both predicted that the atoms would shift by about the same amount (a mass-weighted displacement, $\Delta Q$ , of roughly 0.2 to 0.3 in specific units) and that the light emitted (photoluminescence) would look the same.
The Takeaway: For this type of defect in a "homogeneous" (uniform) material, the simpler TDDFT method is actually good enough. The expensive BSE method confirms it, but the extra cost might not be necessary here.

2. The 2D Boron Nitride (Carbon Dimer)
Next, they moved to a trickier stage: a carbon-dimer defect in 2D hexagonal boron nitride (hBN). This is a flat, single-layer material where the environment is "inhomogeneous"—meaning the way electrons screen each other changes drastically depending on where you are.

The Finding: Here, the old TDDFT method failed spectacularly. It predicted a completely wrong dance move, suggesting the excited state would change its character and the atoms would shift wildly (a $\Delta Q$ of 0.33). The new BSE method, however, showed that the excited state stays localized (stays put) and the atoms shift much less (a $\Delta Q$ of 0.22).
The Takeaway: In materials with uneven environments, the simple TDDFT method gets the physics wrong because it misses the "screened electron-hole interaction." The BSE method, which accounts for how the material's environment shields the electron and hole, corrected the prediction. It showed that the excited state is stable and localized, matching experimental observations much better.

Why This Matters

The paper explicitly rules out the idea that simple, older methods are sufficient for all excited-state problems. While they work fine for uniform materials like diamond, they can lead to qualitative failures (getting the wrong answer entirely) in complex, 2D, or heterogeneous systems.

The authors are confident in their results because they validated their new "Z-vector" approach against standard "finite-difference" methods (which are slow but reliable) and found excellent agreement. They also showed that their method works well whether they start with simple or more complex mathematical descriptions of the electrons.

This work doesn't just solve a math problem; it opens the door to studying self-trapped excitons in perovskites (materials used in solar cells) and other defects in oxides with a level of accuracy previously impossible. By making these calculations feasible for large systems, the team has provided a scalable framework that allows scientists to predict how materials will behave when they glow, conduct, or react under light, paving the way for better quantum sensors and more efficient energy materials. The "dance" of excited electrons can finally be watched in high definition, even in the most crowded and complex ballrooms.

Technical Summary: Analytical Forces from the Bethe-Salpeter Equation for Large-Scale Excited-State Relaxation

Problem Statement
Excited-state (ES) phenomena in materials, such as the vibronic structure of optical spectra and the formation of self-trapped excitons, are governed by the coupling between electronic and vibrational degrees of freedom. Accurate theoretical descriptions require not only excitation energies but also the gradients of ES potential energy surfaces (ES forces). While Density Functional Theory (DFT) and Time-Dependent DFT (TDDFT) are widely used for ES relaxation, they often fail to capture spatially inhomogeneous screening effects, particularly when using semilocal or conventional hybrid functionals.

Many-body perturbation theory, specifically the GW approximation combined with the Bethe-Salpeter equation (BSE), provides a high-level framework for neutral excitations. However, calculating ES forces within BSE has historically been prohibitively expensive for large systems. Conventional approaches based on density-functional perturbation theory (DFPT) require solving $3N_{atom}$ self-consistent Sternheimer equations for each nuclear degree of freedom, leading to steep computational scaling ( $O(N^6)$ to $O(N^7)$ ) and limiting applications to small molecules. Existing Z-vector formulations, which reduce this to a single response equation, have been restricted to Gaussian basis sets for molecular systems, hindering their application to condensed-phase solids.

Methodology
The authors present an efficient plane-wave implementation of analytical nuclear forces for electronic excited states described by the Bethe-Salpeter equation (BSE). The method combines density-matrix perturbation theory (DMPT) with a Lagrangian approach to avoid explicit empty-state summations and the separate response calculations required by conventional DFPT.

Key methodological components include:

Z-Vector Lagrangian Formulation: Instead of solving $3N_{atom}$ Sternheimer equations, the method solves a single Handy–Schaefer Z-vector equation to account for orbital relaxation. This significantly reduces computational complexity.
Low-Rank Dielectric Representation: The implementation utilizes the projective dielectric eigenpotential (PDEP) technique within the random phase approximation. This avoids the explicit construction and inversion of the dielectric matrix and removes the need to converge excitation energies with respect to a large number of unoccupied orbitals.
Algorithmic Accelerations: To achieve scalability, the authors employ Wannier localization for integrals involving the screened Coulomb potential, an inexact Krylov solver for the Z-vector equation, and adaptively compressed exchange for hybrid-functional starting points.
GPU Acceleration: Massively parallel GPU implementations of the PDEP construction and hybrid-functional response solver enable calculations for supercells containing hundreds of atoms.
Approximations: The evaluation of the Lagrangian derivative ( $U$ ) employs two approximations: (i) the derivative of the quasiparticle (QP) Hamiltonian is approximated by that of the Kohn-Sham (KS) Hamiltonian (equivalent to a fixed scissor operator), and (ii) the derivative of the screened Coulomb interaction ( $\delta W$ ) is set to zero. These approximations are validated for specific cases in the study.

Key Results
The method is demonstrated on two point defects with distinct dielectric environments:

Nitrogen-Vacancy Center ( $NV^-$ ) in Diamond:
- This system represents a localized defect in a quasi-homogeneous dielectric host.
- Comparisons between BSE and TDDFT (using both PBE and dielectric-dependent hybrid (DDH) functionals) show that both methods yield consistent excited-state relaxations and similar photoluminescence (PL) spectra.
- The results indicate that for localized deep-level defects in homogeneous hosts, TDDFT based on a hybrid functional can serve as a reliable and cost-effective surrogate for BSE-level structural relaxation.
Carbon-Dimer Defect ( $C_{BC}^0$ ) in 2D Hexagonal Boron Nitride (hBN):
- This system features a spatially inhomogeneous dielectric environment.
- TDDFT with semilocal (PBE) functionals qualitatively fails, predicting a relaxation pattern with a significantly larger displacement and a change in the character of the excited state (from a localized defect excitation to a valence-band-to-defect transition).
- BSE, by explicitly including the screened electron-hole interaction, stabilizes the localized defect excitation. While BSE with PBE still underestimates the displacement compared to hybrid-functional results, BSE with a DDH $\alpha$ starting point yields excitation energies and relaxation patterns in close agreement with experimental values and previous cluster-model calculations.
- The study highlights that for systems with inhomogeneous dielectric screening, the explicit treatment of screened electron-hole interactions (via BSE) is necessary to correct the qualitative failures of semilocal TDDFT.

Significance and Claims
The paper claims to establish a scalable framework for BSE-level studies of excited-state relaxation and vibronic coupling in heterogeneous condensed systems. The primary significance lies in removing the computational barriers that have historically restricted BSE-level ES structural relaxation to small molecules. By combining the Z-vector Lagrangian formulation with GPU acceleration and low-rank dielectric representations, the authors enable ES structural relaxation in supercells containing several hundred atoms for the first time.

The work provides a clear delineation of when BSE forces are essential versus when TDDFT suffices:

For quasi-homogeneous hosts with localized defects (e.g., $NV^-$ in diamond), hybrid-functional TDDFT is sufficient.
For systems with spatially inhomogeneous dielectric environments (e.g., defects in 2D materials), the explicit treatment of screened electron-hole interactions via BSE is critical to obtaining correct excited-state geometries and vibronic line shapes.

The authors position this implementation as a foundation for predictive studies of ES phenomena in complex materials, including self-trapped excitons in metal-halide perovskites and defects in oxides, where electron-hole interactions, dielectric inhomogeneity, and lattice relaxation are strongly coupled. Future work is suggested to lift the fixed-scissor approximation through geometry-dependent QP corrections and to formulate analytical non-adiabatic couplings within the same framework.

The Big Leap: Calculating the "Push" of Excited Electrons

Testing the Dance: Two Different Defects

Why This Matters

More like this