How to choose efficiently the size of the Bethe-Salpeter Equation Hamiltonian for accurate exciton calculations on supercells

TL;DR

This work addresses the prohibitively large cost of solving the Bethe-Salpeter Equation in supercells by introducing a primitive-cell (PC) based workflow to determine the minimal set of valence–conduction transitions needed for converged exciton energies. Central to the approach are Weighted Average Single-Particle Energies (WASPEs) and an energy-window criterion that quantify which transitions contribute meaningfully to each exciton, enabling a dramatic reduction in the BSE Hamiltonian size from zone-folded estimates. Applied to rocksalt LiF, the method yields a BSE kernel with only about 12% of the matrix elements required by naive zone folding, achieving about an 8.4× speed-up with an ~0.15 eV error in the first exciton binding energy, while preserving the GW bandgap within ~0.3% of the PC value. The workflow also demonstrates applicability to defect-like excitations (V$_k$-center) and provides a practical route to study excitons in large supercells, surfaces, and polarons with reduced RAM and time demands. The approach is general and can be adopted in BerkeleyGW and other BSE-capable codes, aided by an accompanying Jupyter notebook for reproducibility and deployment on a broad class of materials.

Abstract

The Bethe-Salpeter Equation (BSE) is the workhorse method to study excitons in materials. The BSE Hamiltonian size, which depends on how many valence-to-conduction band transitions are considered, needs to be chosen to be sufficiently large to converge excitons' energies and wavefunctions but should be minimized to make calculations tractable, as BSE calculations are expensive and scale with the number of atoms as $\mathcal{O}(N_{\rm{atoms}}^6)$. In particular, in the case of supercell (SC) calculations composed of $N_{\rm{rep}}$ replicas of a primitive cell (PC), a natural choice to build this BSE Hamiltonian is to include all transitions derived from PC calculations by zone folding. However, this leads to a very large BSE Hamiltonian, as the number of matrix elements in it is $(N_k N_c N_v)^2$, where $N_k$ is the number of $k$-points and $N_{c(v)}$ is the number of conduction (valence) states used. When creating a SC, the number of $k$-points decreases by a factor $N_{\rm{rep}}$ but both the number of conduction and valence states increase by the same factor, therefore the number of matrix elements in the BSE Hamiltonian increases by a factor $N_{\rm{rep}}^2$, making exactly corresponding calculations prohibitive. Here, we provide a workflow to decide how many transitions are necessary to achieve comparable results, based on only PC results. With our method, we show that to converge the first exciton binding energy of a LiF SC composed of 64 PCs, to an energy tolerance of 0.15 eV, we only need 12\% of the valence-to-conduction matrix elements that result from zone folding with a minimal set of bands. As an example, we use the number of bands from our method to obtain the absorption spectrum of LiF with a V$_k$-like defect. The procedure in our work helps in evaluating excitonic properties in large SC calculations.

Figures (8)

• Figure 1: (a) Diagram showing zone folding of the bandstructure and exciton coefficients. We show a 1D example of bands being zone folded in a BZ of a SC four times the size of the PC. Circle radius is proportional to $|A_{cv {\mathbf{k}}}|$. For the PC calculation, only one valence and one conduction band are necessary to converge this exciton wavefunction, while for the SC, two valence and two conduction bands are necessary to represent this exciton, but not four as would be expected by a zone-folding analysis. WASPEs $\bar{E}_{A,\rm{cond(val)}}$ are given by Eq. \ref{['eq:energy_window_eq']} and give an estimate of the energy scale for the transitions that need to be included. (b) Diagram showing the replication of the PC to create the SC and how this affects the BSE Hamiltonian size. Our method provides a BSE Hamiltonian smaller than the one obtained by zone folding.
• Figure 2: (a) First exciton energy evolution, in comparison to the full PC BSE calculation $\Omega^{\rm{PC}}_1$, as function of the number of valence ($N_v$) and conduction ($N_c$) bands used to build the BSE Hamiltonian. Exciton energy is calculated using Eq. \ref{['eq:omega_part']}. (b) Correlation of $\Omega^{\rm{part}}- \Omega^{\rm{PC}}_1$ with the partial sum of the exciton coefficients, following a quadratic relation.
• Figure 3: Density of States (DOS) at DFT (a) and QP (b) levels for rocksalt LiF calculations in a 2-atom PC (black solid lines) and a 128-atom SC (green dashed lines). Direct bandgap energies at the $\Gamma$-point are reported in Table \ref{['tab:comparison_energy_peaks_PC_SC']}. For DOS calculation, a Gaussian broadening of 0.1 eV was used. (c) DFT and QP bandstructure of LiF, from PC. In all panels, we set the maximum of the valence band to be zero.
• Figure 4: Energy analysis for excitons with energy ranging 11.7 to 20.0 eV. (a-b) Small dots are the conduction (a) and valence (b) WASPEs (Eq. \ref{['eq:energy_window_eq']}) for excitons with excitation energy $\Omega$. Circles are located at the energies of the states that compose the exciton, and the color intensity is proportional to $|A_{cv{\bf k}}|^2$, showing the distribution in energy of the states that compose this exciton. (c) QP bandstructure for PC around $\Gamma$ point. Horizontal green lines are the exciton energies for the first (left) and second (right) bright excitons. Red (blue) lines are the conduction (valence) WASPE corresponding to these excitons.
• Figure 5: (a) Optical absorption for PC and SC cases calculated at BSE and at QP-Random Phase Approximation (RPA) levels. For the SC, we used 137 valence and 31 conduction states to build the BSE Hamiltonian. Absorption peaks from SC calculations are blueshifted by $\sim 0.15$ eV. (b) Convergence of first exciton binding energy using PC results and Eq. \ref{['eq:omega_part']} (see text) in blue, and using SC in red, varying the BSE size (number of matrix elements). The agreement between the two data sets is excellent. (c) Same as (b) but for the second bright exciton. In the top $x$-axis, the zone-folding (ZF) BSE size is $((3\times64)\times(1\times64))^2 \sim1.5 \times 10^8$. The reference binding energies for the PC at full size are in Table \ref{['tab:comparison_energy_peaks_PC_SC']}.