Applying the Liouville–Lanczos method of time-dependent density-functional theory to warm dense matter Open Access

The physical and chemical processes occurring in matter at elevated temperatures and densities are of great interest owing to their relevance in fields such as astrophysics,¹ materials science,² and inertial confinement fusion.^3–5 The state of materials characterized by strongly correlated electrons at temperatures around the electron Fermi temperature is commonly referred to as warm dense matter (WDM).^6–8 

Experiments designed to explore matter under extreme conditions are routinely conducted at various research centers, including the European XFEL in Germany,^9,10 SLAC in the USA,¹¹ SACLA in Japan,¹² the OMEGA laser facility in the USA,¹³ the Z Pulsed Power Facility (Z-Machine) at Albuquerque in the USA,¹⁴ and the National Ignition Facility at Livermore in the USA.¹⁵ The diagnostics of material properties in these experiments is a highly challenging task that often requires significant postprocessing efforts to extract relevant information, and typically involves the use of computer simulations.¹⁶ An important example of such a diagnostic tool is X-ray Thomson scattering (XRTS).^17–19 When the shape of the probing X-ray beam is known (specifically, the combined source-and-instrument function²⁰), XRTS measurements can provide access to the dynamical structure factor (DSF) of the electrons, S(q, ω). This information enables one to infer key information such as the temperature,^21,22 Rayleigh weights,²³ and mass density.^15,23,24

Interpreting the XRTS spectrum necessitates a theoretical analysis, which can be performed using straightforward and computationally efficient approximate models.²⁵ However, these models often result in significant uncertainties in accuracy, due both to the dependence of conditions on the model used²⁶ and to the ability of a wide range of conditions to produce similarly plausible fits to the spectrum.²⁷ Alternatively, high-fidelity (but computationally intensive) ab initio calculations can be employed, such as Kohn–Sham density functional theory (KSDFT) and path integral quantum Monte Carlo (PIMC).¹⁶ In this work, we focus on time-dependent Kohn–Sham density functional theory (TDDFT)²⁸ to calculate the DSF of electrons from first principles.

Two flavors of TDDFT are utilized to calculate the DSF in the WDM regime: one approach is real-time propagation of Kohn–Sham orbitals (RT-TDDFT),^29,30 while the other employs linear-response theory formulated using Kohn–Sham orbitals (LR-TDDFT).^31–34 

In partially ionized WDM, the number of orbitals required in the KSDFT simulations increases rapidly with temperature. In this context, RT-TDDFT is generally considered more computationally favorable than LR-TDDFT for systems with a large number of thermally occupied orbitals.¹⁶ This observation holds true when comparing the conventional LR-TDDFT method, which directly solves Dyson’s equation using the orbital representation. This method requires a significant number of unoccupied bands in addition to the fully and partially occupied orbitals of the equilibrium (ground) state.³⁵ We refer to this approach as standard LR-TDDFT. In the standard LR-TDDFT method, solving Dyson’s equation involves time-consuming inversions and multiplications of large matrices for each frequency value, which is also highly memory-intensive. Consequently, this limits the standard LR-TDDFT method in the WDM regime to relatively small systems, typically of the order of ten atoms in the simulation cell.

An alternative approach for computing the DSF is based on density-matrix representations of LR-TDDFT using the Liouville–Lanczos (LL) approach, which is also known as time-dependent density-functional perturbation theory.^36–38 This method circumvents the expensive matrix inversions and multiplications associated with standard LR-TDDFT.³⁹ Moreover, unlike standard LR-TDDFT, the LL approach does not require additional empty (virtual) bands for calculating the DSF and can access a wide energy (frequency) range, which includes the core-loss region (produced by exciting inner-shell electrons) defined by the utilized pseudopotential. The LL approach for modeling the dynamic density response function of electrons in solids has been implemented as a component in Quantum ESPRESSO under the name turboEELS,³⁹ which uses a recursive Lanczos algorithm demonstrating efficient scalability and parallelization over k-points, band groups, plane-wave schemes, and allows for restarting from previously interrupted calculations.⁴⁰ The LL approach for computing the electronic energy loss spectrum (EELS) under ambient conditions has shown good agreement with standard LR-TDDFT results derived by solving Dyson’s equation for gold,⁴¹ bismuth,⁴² carbon,³⁹ silicon, and aluminum.⁴³ 

Despite its advantages, the LL method has been largely overlooked in the WDM community and has not been utilized to its full potential. To fully exploit the capabilities of the LL method for high-accuracy calculations of the DSF, as well as other related density response and dielectric properties in the WDM regime, the first step is to benchmark and validate the LL approach at elevated temperatures and high densities.

In this work, we benchmark the DSF computed using the LL method against exact PIMC reference results⁴⁴ and standard LR-TDDFT under typical WDM conditions. First, we consider isochorically heated electrons in aluminum (Al) with a face-centered cubic (fcc) lattice. This example is relevant for X-ray-driven heating experiments, characterized by heated electrons in a cold ionic lattice.^45–48 For isochorically heated Al, we test the LL method against the standard LR-TDDFT within the projector augmented-wave (PAW) formalism⁴⁹ by considering DSFs at different wavenumbers and by considering various types of pseudopotentials. The second test case examines partially degenerate equilibrium warm dense hydrogen at metallic and solid hydrogen densities. These types of warm dense hydrogen are relevant for ICF and hydrogen jet experiments.^50,51 Since the computing dynamic properties from PIMC simulations is rather difficult,^52,53 we benchmark the LL approach using exact PIMC data for the imaginary-time density–density correlation function (ITCF) of warm dense hydrogen.⁴⁴ We also provide an analysis of the effect of Lorentzian smearing used in the LR-TDDFT calculations of ITCF. Finally, as an example of the utility of the LL approach, we examine finite-size effects in the simulation of the DSF of the warm dense hydrogen at the considered densities.

The remainder of the paper is organized as follows. In Sec. II, we present a brief discussion of the theoretical methods and outline the simulation details. The results of the simulations and their analysis are discussed in Sec. III. Finally, in Sec. IV, we conclude the paper by summarizing our findings and providing an outlook.

The primary motivation for calculating the DSF of WDM are XRTS measurements, which constitute a key method of diagnostics for extreme states of matter.^18,54,55 In these experiments, the XRTS signal is a convolution of the combined source-and-instrument function R(ω) with the DSF of electrons. Specifically, the function R(ω) takes into account broadening from the detector setup and the characteristics of the X-ray source.²⁰ 

Without loss of generality, the electronic DSF can be broken down into two components: a quasi-elastic contribution from bound electrons and the screening cloud surrounding ions,⁵⁶ and an inelastic contribution that arises from electrons transitioning between available states. In the Chihara approach, these are distinguished into those various combinations of transitions between free and bound states.^26,57,58 The TDDFT method and KS-DFT-based molecular dynamics (MD) of ions allow one to obtain the total electronic DSF, including both the quasi-elastic and inelastic features,³² but unlike the Chihara model, there is no need to distinguish bound and free electrons, since all KS states are treated on equal footing. In practice, this is done by two separate sets of simulations. The quasi-elastic component of the DSF is calculated using the Rayleigh weight,^23,56 which is determined by the static structure factors for electron–ion and ion–ion pairs from KS-DFT MD. The inelastic part of the DSF is then calculated using the TDDFT method. Using this approach under WDM conditions is warranted owing to the significant differences in the characteristic times scales of electron and ion dynamics. In this work, we limit ourselves to the inelastic part of the DSF from LR-TDDFT.

The main computational burden of LR-TDDFT lies in solving Dyson’s equation for the density response function. The standard way is to compute $χGG′0(k,ω)$ and use matrix inversion to solve Dyson’s Eq. (2). An alternative approach is the LL method, which finds the solution by employing an iterative approach derived from time-dependent density-functional perturbation theory using the quantum Liouville equation for the reduced one-electron KS density matrix (see Refs. 37, 39, and 65).

In Dyson’s Eq. (2), calculating the response function $χG,G′0(k,ω)$ requires matrix elements that represent transitions between states with a momentum difference significantly larger than q and an energy eigenvalue difference that is considerably greater than ω.³⁵ This means that a substantial number of additional bands are needed beyond those required to compute the equilibrium (unperturbed) density matrix ρ⁰. By contrast, the Lanczos algorithm used in the LL method relies only on information in ρ⁰ and does not require any extra empty bands.

Theoretically, the standard LR-TDDFT method solving Dyson’s equation and the LL approach should both yield the same solution for χ(q, ω). However, variations in implementation details—such as the choice of basis sets for wavefunctions, handling of tightly bound electrons, and the numerical solvers used—can introduce numerical differences that lead to deviations between the results. As far as we are aware, the LL method for χ(q, ω) has only been implemented in Quantum ESPRESSO,³⁹ which utilizes pseudopotential-based KS-DFT with norm-conserving and ultrasoft pseudopotentials. This implementation is focused on the applications at ambient conditions to model the energy loss spectrum of electrons in solids. To test this implementation of the LL method under WDM conditions, we compare its results for the DSF with the data computed using the standard LR-TDDFT based on the PAW formulation. Furthermore, we validate the LL approach by examining the ITCF results and comparing them with the data from PIMC simulations for warm dense hydrogen.

In practice, the Lorentzian smearing parameter η is introduced into the LR-TDDFT calculations to regularize χ(q, ω) (both in the standard LR-TDDFT and in the LL method). This is done by adding a small imaginary part to the frequency, transforming it from ω to ω + iη. As a result, discrete spectral features (lines) in the DSF are broadened.³⁷ In the real-frequency domain, this technique serves as a method to compute a continuous χ(q, ω), with a negligible impact on the shape and features of the DSF S(q, ω). However, the effect of Lorentzian smearing on the ITCF F(q, τ), computed using S(q, ω) from LR-TDDFT, has not been considered in previous studies. This is a relevant question because XRTS measurements can be analyzed in the Laplace domain, which helps to circumvent the deconvolution problem of the XRTS spectrum with the source-and-instrument function R(ω).^74,79 Therefore, we also examine the impact of Lorentzian smearing employed in LR-TDDFT on the ITCF. We note that the broadening parameter η is a feature of LR-TDDFT and is not needed in the PIMC simulations of the ITCF.

The standard LR-TDDFT calculations based on the solution of Dyson’s Eq. (2) within the PAW method have been performed using the GPAW code,^35,80–84 which is a real-space implementation of the PAW method.⁴⁹ For the LL approach to computing DSF, we used Quantum ESPRESSO.^{39,40,85–87} A static (adiabatic) approximation $KG1,G2xc(k,ω=0)$ for the XC kernel was used in all calculations.

For Al, considering the [100] crystallographic direction, we used a 20 × 20 × 20 k-point grid and an energy cutoff of 500 eV. A conventional unit cell with a lattice parameter of the fcc crystal a = 4.05 was considered.⁸⁸ The electron temperature in Al was set to T = 6 eV. The results were computed using the Lorentzian smearing parameter η = 0.5 eV, which is small enough to distinguish differences between the LL method and the standard LR-TDDFT. The local density approximation (LDA) and generalized gradient approximation (GGA) level XC functionals were considered. We used N_b = 200 KS bands. In GPAW calculations, we utilized the LDA and PBE PAW setups provided by GPAW, which were generated using the Perdew–Wang⁸⁹ and Perdew–Burke–Ernzerhof (PBE)⁹⁰ functionals, respectively. For the LL calculations, the Perdew–Zunger LDA functional⁹¹ was used with the Al.pz-vbc.UPF pseudopotential, and the GGA level PBE functional was used in combination with the norm-conserving pseudopotential Al.pbe-rrkj.UPF. Both pseudopotentials are from the Quantum ESPRESSO pseudopotential database.⁹² 

For hydrogen, we tested various types of pseudopotentials in combination with the PBE XC functional. In the calculations using Quantum ESPRESSO, the used pseudopotentials are norm-conserving H.pbe-vbc.UPF⁹² and ultrasoft H_PBE_v1.4.USPP.F.UPF.⁹³ In addition, the PAW method was used with the PBE setup for hydrogen from GPAW. We considered 14, 32, and 100 protons in the simulation cell. As a part of the analysis of the behavior of the ITCF computed from the DSF values, we varied the Lorentzian smearing parameter in the range from η = 0.2–3.0 eV. The DSFs were computed for warm dense hydrogen in equilibrium at densities ρ = 0.08 g/cm³ (T = 4.8 eV) and 0.34 g/cm³ (12.528 eV). These correspond to density parameters r_s = 3.23 and 2.0, with r_s being the ratio of the mean interparticle distance to the Bohr radius (Wigner–Seitz radius).⁹⁴ The considered temperatures provide a partially degenerate state with the system temperature T equal to the Fermi temperature of electrons T_F. We choose these densities and temperatures to benchmark TDDFT results against the available PIMC results for the ITCF.⁴⁴ The size of the simulation box is determined by the density r_s (or ρ) and the number of protons N, expressed as $L=rs(43πN)1/3$ (in atomic units). We set the Monkhorst–Pack k-point grid to 10 × 10 × 10 for N = 14, 8 × 8 × 8 for N = 32, and 2 × 2 × 2 for N = 100. The DSF results were averaged over 10 snapshots for N = 14 and N = 32 at r_s = 2. For N = 100 at r_s = 2, the DSF values were averaged over five snapshots. At r_s = 3.23, we computed the DSF values using N = 14 and 32 protons with averaging performed over 20 snapshots. The snapshots of the position of protons were chosen randomly from the particle trajectory obtained by the KSDFT-based MD simulations. In all calculations of the DSF, we used a bi-constant extrapolation scheme after the computed number of LL iterations N_iter to obtain 10⁴ Lanczos coefficients.³⁹ For Al, we used N_iter = 3 × 10³ and for warm dense hydrogen, N_iter = 1.2 × 10⁴.

We begin by analyzing the DSF of isochorically heated Al, which has cold ions arranged in the fcc lattice structure and hot electrons at a temperature T = 6 eV. Such a transient state is generated and examined in experiments that utilize X-ray-driven heating used for XRTS measurements.^46,47,95,96 In these experiments, the incident X rays do not interact directly with the nuclei, keeping them cold during a free-electron laser pulse that lasts less than 100 fs.

The simulations were performed for wavenumbers $q=0.155Å−1$⁠, $0.7757Å−1$⁠, and $1.396Å−1$ along the crystallographic direction [100]. In the standard LR-TDDFT method, the wavenumbers must correspond to the difference between two k-points. The LL method does not have this constraint. We note that the greater flexibility in choosing q values in the LL method provides an additional advantage for analysis of the experimental XRTS data, where a wavenumber blurring effect occurs owing to the finite size of a detector.⁹⁷ To compare the results from Quantum ESPRESSO and GPAW, we select the wavenumber values in the LL calculations to correspond with those derived from the differences between two k-points at the specified parameters.

In Fig. 1, we show the convergence of the results for isochorically heated Al with respect to the number of iterations N_iter in the LL calculations at $q=1.396Å−1$⁠. It is evident that the LL results have converged after N_iter = 2000 iterations. In Fig. 2, considering the LDA-level description of the XC functional, the LL results computed using N_iter = 3000 are compared with the exact solution of Dyson’s equation within the PAW formalism. Overall, we find good agreement between the different approaches. At $q=0.155Å−1$⁠, one can see a well-defined plasmon feature, which becomes broader with increasing wavenumber, indicating a stronger Landau damping effect.⁹⁸ The pair continuum (strong plasmon damping) for Al corresponds to $q>1.3Å−1$⁠, where the plasmon is no longer well defined.^59,99 The small spikes at ω < 10 eV for $q=0.7757Å−1$ and $1.396Å−1$ are due to lattice effects leading to energy gaps at Brillouin zone boundaries.¹⁰⁰ An interesting observation from Fig. 2 is a noticeable difference in the magnitude of the DSF maximum for $q=0.7757Å−1$ between the LL method and Dyson’s equation solution within the PAW formalism. To analyze this aspect further, we compare the results of the DSF computed using the PBE XC functional with those obtained using LDA-level functionals, as shown in Fig. 3. The calculations were performed using norm-conserving pseudopotentials and PAW setups consistent with the XC functionals used. The PBE-based results from the LL method are in close agreement with the LDA-level calculations using the same method. When examining the DSFs obtained by solving Dyson’s equation within the PAW framework, we find that the maximum of the DSFs calculated with the PBE functional are larger than those derived from the LDA-level data. This larger value leads to closer agreement with the results obtained from the LL method. We conclude that the observed differences in the vicinity of the DSF maximum are due to the differing methods (pseudopotentials) used to handle core electrons.

As an example of WDM in an equilibrium state, we examine warm dense hydrogen. We start with a metallic density with r_s = 2 and at a degeneracy parameter θ = T/T_F = 1. This is a typical WDM regime where correlations, thermal excitations, and quantum effects cannot be neglected. In Fig. 4, we show results for one snapshot of 14 protons at $q=1.376Å−1$⁠. Figure 4(a) shows the considered snapshot along with the density distribution of electrons n, expressed in terms of the mean electron density n₀. We can observe the characteristic electron structure in WDM, where the distinction between electrons that are localized around ions and those that are quasi-free is not clearly defined.¹⁰¹ In Fig. 4(b), we illustrate the convergence of the results with respect to the number of iterations N_iter in the LL calculations at $q=1.376Å−1$ and smearing parameter η = 0.5 eV, where one can see that the results have converged after N_iter = 8000 iterations. In Fig. 4(c), we provide a comparison of the data computed using norm-conserving and ultrasoft pseudopotentials in the LL method with the results computed by solving Dyson’s equation using the PAW approach. From this comparison, we see perfect agreement between the different approaches considered at 0 eV ≤ ω ≤ 80 eV. Using a logarithmic scale, in Fig. 5, we show the same data in the range 0 eV ≤ ω ≤ 500 eV. From Fig. 5, it can be seen that the DSF calculated using the standard LR-TDDFT method shows a drop starting around 150 eV. This drop results in a strong deviation from the results obtained using the LL approach. Besides that, the standard LR-TDDFT is limited to 180 eV at the number of orbitals used (N_b = 200). These features are attributed to the dependence on the number of unoccupied states in the standard LR-TDDFT method, where the maximum value of ω is limited to the largest difference in the energy of computed states. By contrast, the LL method does not require extra empty states for computing the DSF at large frequencies (energies) ω, as demonstrated in Fig. 5.

The demonstrated limitation in the frequency ranged caused by the finite number of empty states in the standard LT-TDDFT method restricts its use in predicting the DSF to relatively small energy losses. In practice, this means that many interesting features and cases are beyond the reach of the standard LR-TDDFT. For example, all elements beyond boron have at least one excitation edge beyond 200 eV.^102,103 Similarly, examining the complete DSF at high q is inaccessible to standard LR-TDDFT simulations. An increase in q leads to a shift of the DSF maximum position ω₀ to larger values, which leads to an increase in the number of required empty orbitals in the standard LR-TDDFT. In WDM, as well as metals and some nonmetallic systems, where electrons exist in or are easily excited into the conduction band, the energy of the DSF maximum, ω₀, shows $ω02−ωp2∼q2$ scaling at small wavenumbers, transitioning to a quadratic dependence ω₀ ∼ q² at large wavenumbers,^97,104–106 where ω_p denotes the electron plasma frequency. To accurately compute the DSF, it is necessary to have a sufficient number of empty bands whose energy differences cover the range ɛ ≫ ω₀. By using the density of states for a free-electron gas as an approximation, the number of electron orbitals within a specified energy range ɛ can be estimated as N_b ∼ ɛ^3/2, leading to cubic scaling N_b ∼ q³ at large wavenumbers. This cubic scaling of the required number of bands results in significant memory demands for parallel computations in the standard LR-TDDFT method. This makes the standard method impractical for modeling the DSF at large wavenumbers using a meaningful frequency grid and range. For example, at very high q, where the single-particle limit is being probed, the DSF is characterized by a Compton feature, which can be positioned at frequencies in excess of 100 eV and have a width of more than 100 eV (see, e.g., Ref. 55), which would require a huge number of bands to model. Simulating the DSF at very high q with standard LR-TDDFT would therefore be even more challenging. Furthermore, modeling the DSF of disordered systems incurs additional computational costs due to the need to average over multiple snapshots.

An example of averaging over snapshots is illustrated in Fig. 6(a) for $q=1.528Å−1$ using η = 0.2 eV. Here, we present the DSF values computed from ten different snapshots, along with the corresponding averaged results. The results for warm dense hydrogen are compared with data from the uniform electron gas (UEG) model, which was computed using an effective static approximation for the local field correction.^107,108 From Fig. 6(a), we observe that the UEG model significantly overestimates the DSF magnitude at its peak and exhibits a faster decay at frequencies ω beyond the DSF maximum. In Fig. 6(b), we compare the shifted ITCF $F̃(q,τ)$ [as defined in Eq. (9)] calculated using the DSF values obtained from LL-based LR-TDDFT simulations—both averaged and from individual snapshots—against the $F̃(q,τ)$ from the PIMC simulations of warm dense hydrogen⁴⁴ and the data from the PIMC simulation of the UEG. This comparison reveals that the averaged LR-TDDFT data for $F̃(q,τ)$ align excellently with the PIMC results for warm dense hydrogen. By contrast, the UEG model exhibits substantial differences from the results for warm dense hydrogen, particularly around τ = 0.5β.

The shape of the LR-TDDFT results depends on the chosen Lorentzian smearing parameter η. A larger η leads to a stronger broadening of the DSF results. To demonstrate that the observed agreement between the LR-TDDFT results (calculated using the LL method) and the PIMC data for $F̃(q,τ)$ is not an artifact of the selected η value, we investigate the effects of varying η on S(q, ω) and $F̃(q,τ)$⁠. The results for S(q, ω) and $F̃(q,τ)$⁠, averaged over snapshots, for the range 0.1 eV ≤ η ≤ 0.9 eV are presented in Fig. 7. It is evident that decreasing η leads to noisier results for the DSF, alongside a greater magnitude of the maximum in the DSF. However, for $F̃(q,τ)$⁠, a decrease in η results in closer agreement between the LR-TDDFT results obtained using the LL method and the PIMC data. It is noteworthy that the Laplace transform acts as a noise filter,¹⁰⁹ resulting in smooth ITCF results even at smaller η values.

As previously mentioned, the LL approach to LR-TDDFT is particularly valuable for DSF calculations at large wavenumbers. To test the application of the LL method for large q values, we conducted simulations for $q=4.584Å−1$⁠. The results for this wavenumber, at r_s = 2 and θ = 1, are presented in Fig. 8. In Fig. 8(a), we show results for ten individual snapshots at η = 2 eV, along with the averaged value over these snapshots and data from the UEG model. It is evident that the differences between the DSFs of individual snapshots and the averaged value are much less pronounced compared with the case of $q=1.528Å−1$⁠. Furthermore, at $q=4.584Å−1$⁠, the deviation of the UEG model results from the LR-TDDFT data is significantly smaller than that observed at $q=1.528Å−1$⁠. Both of these findings are readily explained in terms of the characteristic length scale λ = 2π/q. Specifically, the limit of large q corresponds to the structureless single-particle regime, where correlations between electrons and nuclei effectively vanish.

In Fig. 8(b), we present the results for the snapshot-averaged DSF at $q=4.584Å−1$ for 0.1 eV ≤ η ≤ 3.0 eV. Notably, we observe significantly stronger quasi-periodic noise in the DSF for η < 1 eV compared with that at $q=1.528Å−1$ (see Fig. 7). A smooth curve for the DSF is produced using η ≳ 2 eV. We note that the magnitude of η should be gauged taking into account the characteristic frequency (energy) scale ω of the DSF. For example, the full width at half maximum (FWHM), denoted by Δω₀, at $q=4.584Å−1$ is ∼3.5 times larger than that at $q=1.528Å−1$⁠. The relative value of η/Δω required for damping the noise in the DSF is similar for both wavenumbers, with an optimal value around η₀ = 10⁻²Δω₀. This degree of smearing results in a minor effect on the overall shape of the DSF, with the magnitude of the DSF maximum being underestimated by about 5%. At η ≲ η₀, the FWHM Δω₀ and the position of the DSF maximum, ω₀, remain nearly unchanged for both $q=4.584Å−1$ and $q=1.528Å−1$⁠. An increase in η leads to a gradually deterioration in the DSF maximum prediction, owing to the DSF asymmetry relative to ω₀.

In Fig. 8(c), we show the shifted ITCF $F̃(q,τ)$ data computed using the LL method, the PIMC results for warm dense hydrogen, and the ITCF data for the UEG model from PIMC simulations at $q=4.584Å−1$ (with r_s = 2 and θ = 1). It is evident that the agreement between the LR-TDDFT results and the PIMC data for $F̃(q,τ)$ of warm dense hydrogen is very good when η < 0.5 eV. However, the UEG model shows significant deviations from both of these results. As the smearing parameter increases, the LR-TDDFT results computed using the LL method diverge from the PIMC data for hydrogen. This indicates that converged results for the ITCF can be obtained by decreasing the value of η, despite the noise that arises in the DSF.

We validated the results computed using the LL method by considering a system with N = 14 protons, for which the PIMC data are available. Let us now test the effect of the system size on the DSF. In Fig. 9, we compare the DSF results at r_s = 2 and T = 12.58 eV computed using 14, 32, and 100 protons at $q=1.528Å−1$ and $4.584Å−1$⁠. The overall shape of the DSF, the position of the maximum, and the FWHM are in good agreement for the simulations with different numbers of particles.

Next, we test the LL approach to LR-TDDFT by examining the more challenging case of warm dense hydrogen at solid density, with r_s = 3.23 and T = 4.8 eV. This regime has stronger coupling between particles and enhanced localization of electrons around protons compared with the case with r_s = 2.¹¹⁰ 

In Figs. 10 and 11, we present the LL-based LR-TDDFT results for S(q, ω) and $F̃(q,τ)$ at $q=0.94617Å−1$ (with η = 0.3 eV) and $q=2.838Å−1$ (with η = 0.5 eV), respectively. The calculations utilized 20 snapshots, and we include the corresponding snapshot-averaged values. The LR-TDDFT results are compared with the PIMC results for warm dense hydrogen from Ref. 44. Additionally, we provide data computed using the UEG model.

It is evident from Figs. 10 and 11 that the LR-TDDFT data averaged over snapshots for $F̃(q,τ)$ align excellently with the PIMC results at both wavenumbers considered. By contrast, the UEG model exhibits substantial discrepancies when compared with the LR-TDDFT and PIMC results for hydrogen. We note that S(q, ω) and $F̃(q,τ)$ in this case exhibit similar dependences on η, as discussed in Sec. III B. These findings for solid density hydrogen further validate the LL approach to LR-TDDFT under warm dense matter conditions, characterized by strong localization of electrons around ions.^101,111

To investigate the effect of system size on the DSF of hydrogen at r_s = 3.23, we compare the results for 14 particles with those obtained using 32 particles at $q=0.94617Å−1$ and $2.838Å−1$ in Fig. 12. The results indicate that increasing the wavenumber reduces the sensitivity of the results to system size as we eventually reach the single-particle limit. The DSF data for 14 and 32 particles show close agreement at $q=2.838Å−1$⁠. At $q=0.94617Å−1$⁠, the DSF results for 14 and 32 particles exhibit only minor deviations from each other to the left of the DSF maximum, with the data for 14 particles being somewhat noisier.

We have validated and tested the LL method, as implemented in Quantum ESPRESSO (compatible with both CPUs and GPUs), for WDM applications. This step is crucial for expanding the resources available for modeling matter under extreme conditions, given the importance of XRTS measurements as a diagnostic tool in WDM experiments.

To test the LL method, we first considered isochorically heated Al, which features hot electrons in an fcc lattice structure of cold ions. Such transient states are generated and examined in experiments that utilize X-ray-driven heating^46,47,95,96 (see also Ref. 48 for detailed discussions). For the fcc Al, we found good overall agreement between the results from the LL method and the solution of Dyson’s equation in the standard LR-TDDFT within the PAW formalism. However, we noted differences in the DSF values near the DSF maximum between the two methods. This discrepancy arises from the different approaches in approximating the interaction of core electrons with valence electrons. Specifically, the pseudopotential or PAW setup employed defines the interaction of valence electrons with the crystal lattice, affecting the band structure and consequently the DSF as well as EELS.^48,112–115 Our analysis of the DSF for fcc Al demonstrated that if relatively subtle features of the DSF are of interest, particular care must be taken by cross-validating results using different approaches to TDDFT calculations. This has also been pointed out in a recent work by Hentschel et al.¹¹⁶ considering pseudopotentials with 3 and 11 valence electrons.

Second, we considered warm dense hydrogen as an example of WDM in equilibrium. We tested the results from the LL method by comparing them with LR-TDDFT data computed within the PAW formalism. Additionally, by examining the ITCF derived from the DSF results, we confirmed the high quality of the LL method by benchmarking the data against the PIMC results for hydrogen. This opens the door for further detailed studies of warm dense hydrogen properties such as the DSF, density response functions, and the dynamic dielectric function using the LL method across a wide range of parameters. This direction of research is important given the relevance of warm dense hydrogen for ICF and astrophysical applications.

Moreover, our analysis of the effects of pseudopotentials, the number of iterations in the LL method, and the Lorentzian smearing parameter provides a valuable foundation for future applications of the LL method in the WDM regime.

The primary computational limitation of standard LR-TDDFT arises from the memory requirements to store matrices on simulation cores.¹¹⁷ This can reach several terabytes when around 10³ bands are used. When a relatively small number of bands is required, such as for wavenumbers $q≲qqF$ and small-sized systems (typically containing about N ∼ 10 atoms), the standard LR-TDDFT approach can be up to two orders of magnitude more efficient than the LL method. However, in cases involving a large number of bands or atoms, as discussed in this work, employing standard LR-TDDFT becomes computationally unfeasible or highly impractical. By contrast, the LL method has lower memory requirements and does not require empty states, allowing efficient parallelization over thousands of cores. This makes it capable of handling larger wavenumbers, broader frequency ranges, and larger systems. A detailed comparison of computational costs between the LL method and LR-TDDFT, using bulk diamond as an example, can be found in Ref. 39.

We would like to emphasize the utility of the LL method for calculating the DSF at large wavenumbers. In the standard approach to solving Dyson’s equation within LR-TDDFT, increasing the wavenumber requires a significant increase in the number of empty bands and memory resources for calculations. Consequently, the standard LR-TDDFT method faces considerable challenges when calculating the DSF at large wavenumbers. By contrast, the LL method does not require empty orbitals, enabling the DSF calculations at large wavenumbers. The XRTS spectra at large q values are particularly important for WDM experiments, since backward scattering geometry measurements are essential for diagnostics at key facilities, including the NIF,^15,118 which operates the most powerful laser platform for WDM studies in the world.

In summary, the LL approach to LR-TDDFT is valuable for ab initio simulation of WDM properties, with applications in XRTS analysis and the calculation of dielectric properties. For instance, the LL method can be employed to plan and design new experiments. Furthermore, it can be used to test other TDDFT implementations, such as stochastic RT-TDDFT¹¹⁹ and orbital-free TDDFT,^119–121 for calculating the linear response properties of WDM across various parameters. Finally, we note that LR-TDDFT provides access to a wide range of density response functions (including the ideal Kohn–Sham response and the random phase approximation with and without local field effects), which is helpful for methodological advances in fundamental WDM theory.

This work was partially supported by the Center for Advanced Systems Understanding (CASUS), financed by Germany’s Federal Ministry of Education and Research (BMBF) and the Saxon State Government out of the State Budget approved by the Saxon State Parliament. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2022 research and innovation programme (Grant Agreement No. 101076233, “PREXTREME”). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work has received funding from the European Union’s Just Transition Fund (JTF) within the project Röntgenlaser-Optimierung der Laserfusion (ROLF), Contract No. 5086999001, co-financed by the Saxon State Government out of the State Budget approved by the Saxon State Parliament. Computations were performed on a Bull Cluster at the Center for Information Services and High-Performance Computing (ZIH) at Technische Universität Dresden and at the Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen (HLRN) under Grant No. mvp00024.

The authors have no conflicts to disclose.

Zhandos A. Moldabekov: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Sebastian Schwalbe: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Thomas Gawne: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Thomas R. Preston: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Jan Vorberger: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Tobias Dornheim: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Validation (equal); Writing – original draft (equal).

The data that support the findings of this study are available within the article.