X-Ray Thomson scattering without the Chihara decomposition

X-Ray Thomson Scattering (XRTS) is an important experimental technique used to measure the temperature, ionization state, structure, and density of warm dense matter (WDM). The fundamental property probed in these experiments is the electronic dynamic structure factor (DSF). In most models, this is decomposed into three terms [Chihara, J. Phys. F: Metal Phys. {\bf 17}, 295 (1987)] representing the response of tightly bound, loosely bound, and free electrons. Accompanying this decomposition is the classification of electrons as either bound or free, which is useful for gapped and cold systems but becomes increasingly questionable as temperatures and pressures increase into the WDM regime. In this work we provide unambiguous first principles calculations of the dynamic structure factor of warm dense beryllium, independent of the Chihara form, by treating bound and free states under a single formalism. The computational approach is real-time finite-temperature time-dependent density functional theory (TDDFT) being applied here for the first time to WDM. We compare results from TDDFT to Chihara-based calculations for experimentally relevant conditions in shock-compressed beryllium.

Warm dense matter (WDM) arises in many contexts ranging from planetary science [1][2][3] to the implosion stage of inertial confinement fusion [4][5][6]. While there are no sharp pressure, temperature and density boundaries for the WDM regime, it is generally viewed as an intermediate state between a condensed phase and an ideal plasma where Fermi degeneracy is present, and the Coulomb coupling and thermal energy are comparable in magnitude [7].
Experimental characterization of warm dense matter is challenging due to the difficulty of producing uniform samples at extreme conditions and developing diagnostic techniques that can provide accurate and independent measurements of these conditions for transient samples opaque to optical photons. X-ray Thomson Scattering (XRTS) [7,8], one such diagnostic technique, exploits the scattering of hard coherent x-rays to directly probe the system's dynamic structure factor (DSF). Through the fluctuation-dissipation theorem, the DSF is related to the system's density-density response, and consequently XRTS provides direct insight into electron dynamics.
XRTS experiments have been performed on a variety of materials including beryllium [9,10], lithium [11], carbon [12], CH shells [13], and aluminum [14,15]. With recent improvements in source brightness [15] producing increasingly high resolution and high signal-to-noise data, the full DSF is expected to become routinely available. In anticipation of these advances, it is critical and timely to examine the theoretical constructs underpinning the interpretation of these experiments.
The most common model of XRTS experiments relies on an additive form of the DSF due to Chihara [16,17]: S(q, ω) = |f I (q)+ρ(q)| 2 S ii (q, ω)+Z f S ee (q, ω)+S bf (q, ω) (1) The DSF varies with momentum and energy transfers (q and ω) and is partitioned into 3 features that can be interpreted in terms of x-ray scattering processes. These include scattering from electrons bound to and adiabatically following ions (|f I (q) + ρ(q)| 2 S ii (q, ω)), from Z f free electrons per ion (Z f S ee (q, ω)), and from bound electrons that are photo-ionized (S bf (q, ω)). While successfully applied to many systems, this model relies on numerous approximations and assumptions. Most critically, the electrons are separated into bound and free populations, a distinction that is often ambiguous in the WDM regime. Each term in Eqn. 1 is subject to different models potentially leading to under-constrained fits to experimental data [18]. The ionic feature typically relies upon a decomposition into a product of an ion-ion structure factor, S ii (q, ω), and an average atomic form factor, f I (q) + ρ(q), with the first term describing the unscreened bound electrons and the latter the screening cloud, which must be treated carefully [19]. However, recent work has focused on moving past this decomposition of the elastic peak [20].
In this work, we transcend the Chihara decomposition by explicitly simulating the real-time dynamics of warm dense matter using a finite temperature form of timedependent density functional theory (TDDFT) [21,22] and the Projector Augmented-Wave (PAW) formalism [23,24]. In PAW, the all-electron Kohn-Sham orbitals (and their associated density) can be accessed via an explicit linear transformation on the smoother pseudo orbitals. By performing calculations with none of the core states frozen, we avoid making any assumptions between bound and free electrons, and treat all electrons similarly.
We work with a real-time implementation of TDDFT in the Vienna Ab-Initio Simulation Package (VASP) [24][25][26] (see Supplemental Material [27]) that provides a number of attractive features. Physically, higher-order response phenomena and Ehrenfest molecular dynamics are accessible in this framework. Computationally, the orthogonalization bottleneck that limits standard DFT approaches is removed, as it is only explicitly required arXiv:1512.05795v2 [cond-mat.mtrl-sci] 1 Apr 2016 for the calculation of the initial state of the Kohn-Sham orbitals. This leads to excellent strong scaling [28,29], and we have observed near-perfect scaling up to 65,536 cores in our implementation.
We next outline the details of our TDDFT calculations, noting that we use Hartree atomic units (m e = e 2 = = 1/4π 0 = 1) unless otherwise indicated. Time-dependent quantities evaluated at t = 0 are indicated by the addition of a subscript 0 and the absence of a temporal argument. Fourier transformed quantities are indicated by a diacritic tilde. All calculations are spin unpolarized.
The equation of motion in real-time TDDFT is the time-dependent Kohn-Sham (TDKS) equation: . The orbitals, indexed by band and Bloch wave number, are such that their weighted sum produces the time-dependent density, ρ(r, t). The external potential, v ext , includes contributions due to the Coulomb field of the bare nuclei as well as a model of the x-ray probe, v probe (r, t), which is quiescent until t = 0. v H and v xc are the Hartree and exchange-correlation potentials, with accurate and efficiently computable approximations to the latter being a central theoretical concern of TDDFT. The initial conditions for the {φ n,k (r, t)} from Eqn. 2 are the self-consistent solution to a Kohn-Sham Mermin DFT calculation [30] at electron temperature T e in a supercell of volume Ω sc . The initial density is then: where f n,k is a composite weight consisting of the measure of the specific Bloch orbital weighting and the Mermin weight encoding temperature dependence according to a Fermi-Dirac distribution. It is important to note the implicit dependence of these initial orbitals on the equilibrium electron temperature, T e . Evolving these orbitals under Eqn. 2, the time-dependent density becomes: The weights are not time-evolved, as we do not expect them to change in the linear response regime [31]. For stronger perturbations, additional formalism might be required for the weights. That the Mermin formalism is sensible within TDDFT in the linear response regime is supported by recent foundational work [32]. The action of a probe potential, v probe (r, t), turned on at t = 0 leads to a change in the time-dependent density. The linear density-density response function, χ ρρ (r, r , t), encodes the relationship between these two quantities: δρ(r, t) = ∞ 0 dτ Ωsc dr χ ρρ (r, r , τ )δv ext (r , t − τ ) (5) where we use the notation δf (r, t) = f (r, t) − f 0 (r) for δρ(r, t) and δv ext (r, t). If we fix ionic positions at time t = 0, which evolve slowly on the relevant attosecond timescales, then δv ext (r, t) = v probe (r, t) and we can construct a probe potential that can be used to extract the DSF, similar to Sakko, et. al. [33]. The real-time density response to such a probe potential is shown in Fig. 1 FIG. 1. The density response of warm dense beryllium (density 5.5 g/cm 3 and Te = 13 eV, movie in Supplemental Material [27]) due to (a) a perturbing potential with the illustrated envelope observed at times coinciding with (b) the peak of the perturbation, (c) the peak of the density response, and (d) the plasmons continuing to ring around the system after the perturbation has ceased. Red (blue) isosurfaces bound volumes of charge accumulation (depletion) and yellow isosurfaces indicate the nodal surface. The amplitude of the perturbation is within the linear response regime, a discussion of which is in the Supplemental Material [27].
In principle, any sufficiently weak analytic probe potential will allow us to extract the response function. For convenience, we choose v probe (r, t) = v 0 e iq·r f (t), where f (t) is a Gaussian envelope and v 0 is related to the probe intensity. The Fourier transformed response function, χ ρρ (q, −q, ω) = δρ(q, ω)/v 0f (ω), is then related to the DSF through the fluctuation-dissipation theorem: We are careful to note that Fourier transforms are normalized such that δρ(q, ω) has units of inverse frequency. As a proof of principle, we report calculations of the DSF for 3x-compressed beryllium consistent with the conditions reported in [10]. We consider the same momentum transfers as in [34] to study a range of excitations spanning the collective regime to the beginning of noncollective regime. For convenience of presentation, each q value is mapped onto an XRTS scattering angle, θ, relative to the 2Å probe wavelength in [10] (see Supplemental Material [27] for more information). Our results come from averaging the response of electronic densities generated from several static uncorrelated ionic configurations sampled from thermally equilibrated DFT-MD calculations. These calculations were performed on 32 and 64 atom supercells with a four electron beryllium PAW potential within the local density approximation (LDA), with electrons and ions thermostatted at T = 13 eV, and k-point sampling at ( 1 4 , 1 4 , 1 4 ), analogous to the Baldereschi mean-value point for cubic supercells. For these conditions a plane wave cutoff of 1400 eV was required to converge the pressure to within 1%, and 576 (32 atom)/1152 (64 atom) Kohn-Sham orbitals were needed to represent the thermal occupation (f n,k (T e ) ≥ 10 −5 ).
Each sample configuration is used to seed a TDDFT calculation of the DSF, utilizing the same cutoff and number of orbitals. The initial Mermin electronic state is recomputed using a denser k-point sampling on a 3×4×4 (32 atom) or 2×2×2 (64 atom) Monkhorst-Pack grid. To assess the effect of the frozen core approximation (FCA), we consider electronic initial conditions and dynamics generated using both two and four electron beryllium PAW potentials. Results of our calculations are illustrated in Figs. 2 and 3. Details of the averaging procedure used to generate results, and information concerning the satisfaction of sum rules, can be found in the Supplemental Material [27]. ∆t=4 as ∆t=2 as ∆t=1 as ∆t=0.5 as FIG. 2. The DSF of warm dense beryllium (density 5.5 g/cm 3 and Te = 13 eV) at the scattering angles (θ) considered in [34] with (dashed lines) and without (marked lines) the frozen core approximation. All TDDFT calculations utilize the adiabatic LDA. Fig. 2 directly compares the DSF computed with and without the FCA. While the PAW method still includes the proper all-electron density in aggregate, only the orbitals tied to the two outermost valence electrons are included in the time-evolved response within the FCA. This effectively removes the dynamics of the core states from the density response and the high energy shoulder above 80 eV in the DSF is removed. For the temperatures and densities being considered, this is roughly equivalent to partitioning the inner two and outer two electrons into bound and free groups in the Chihara picture, such that the two-electron response corresponds to S ee (q, ω). However, there are important distinctions to keep in mind. First, in the four-electron calculation, all electrons are being treated identically, whereas it is typical to treat S ee (q, ω) and S bf (q, ω) using different levels of theory and without self-consistency in the Chihara framework. Second, even in the two-electron calculation, the response of the outer two electrons is still aware of the two frozen core states tied to each atom through their screening of the nuclear potential. Finally, these calculations are based upon explicit simulations of the realtime electron dynamics of a bulk supercell of warm dense beryllium rather than a phenomenological model of the response based upon a jellium plus average-atom picture.
Based upon our observation that the two-electron response roughly corresponds to S ee (q, ω), we can also extract a quantity akin to S bf (q, ω) by differencing the fourelectron and two-electron DSFs. The effective S ee (q, ω) and S bf (q, ω) computed within TDDFT are illustrated in Fig. 3. Here we compare our TDDFT calculations to calculations done using state-of-the-art models for S ee (q, ω) and S bf (q, ω). The former is treated with an RPAlevel model dielectric function with lifetime effects taken from four-electron DFT-MD calculations of the optical conductivity; the Mermin approximation-ab initio collision frequencies (MA-AICF) method in [34]. The latter is calculated with the formalism developed in [35] and a quantum mechanical average-atom ion-sphere model with Slater exchange. As we are interested in studying energies relevant to the electronic response, we ignore S ii (q, ω), though it is necessarily present in the Chiharaindependent TDDFT calculation.
Examining the dispersion of the primary plasmon peak in Fig. 3a, we see that TDDFT predicts a slight (∼5 eV) blue shift relative to the MA-AICF calculation θ = 20 • and 40 • , whereas it predicts a stronger (∼10 eV) red shift relative to MA-AICF at θ = 60 • and 90 • . We attribute this shift to exchange/correlation and band structure effects, not present in the MA-AICF dielectric function. Previously, comparisons of inelastic x-ray scattering spectra to RPA and LDA in cold free electron metals (Na and Al) indicate a similar trend in which both LDA and experiment are red shifted relative to RPA [36]. However, these calculations also indicate that the addition of lifetime effects to the LDA are necessary to totally reconcile theory and experiment. While non-adiabatic exchange IG. 3. Comparison of the DSF of warm dense beryllium computed using TDDFT (marked lines) and individual terms in a Chihara decomposed theory (dashed lines). (a) Illustrates the two-electron TDDFT response compared to the MA-AICF method [34] for Z f See(q, ω). (b) Illustrates the difference of the four-and two-electron TDDFT responses compared to an average atom treatment [35] of S bf (q, ω). correlation kernels are available for energy domain linear response TDDFT, no such time-domain exchange correlation potentials are currently available for real-time TDDFT. As warm dense matter requires a large number of thermally occupied states such that energy domain TDDFT may become computationally prohibitive, warm dense matter may provide compelling motivation for the development and testing of these functionals.
Considering the bound-free shoulder in Fig. 3b, we see that TDDFT is generally in good agreement with the average-atom model treatment of S bf (q, ω) with some minor differences. Such average atom models agree well with real-space Green's function methods for cold solid beryllium [37]. We observe a trend opposite the free-free feature, in which there is a red shift of the TDDFT result at small angles, and a blue shift at large angles. We expect that LDA will do an increasingly poor job of describing the Compton scattering limit at large θ due to its well-established self-interaction error. Applying a func-tional with some fraction of Fock exchange should remedy this behavior and will be the subject of future investigations. Further, the TDDFT bound-free feature computed by differencing the two-and four-electron DSFs has a small negative dip below the 80 eV onset of the core feature. It is difficult to determine whether this is due to core-polarization suppressing the response of the valence electrons in the four-electron calculation, or potentially an artifact of the different pseudization procedures used to generate the two PAW potentials used for this comparison. This motivates further investigations of the PAW formalism applied to both real-time TDDFT and specifically, warm dense matter in which thermal effects will start to blur the line between core and valence electrons.
We presented a method for the direct calculation of the DSF for warm dense matter, independent of the Chihara decomposition, by applying real-time TDDFT to configurations drawn from thermal Mermin DFT-MD calculations. Comparison of our results with state-ofthe-art models applied within the Chihara picture illustrates some subtle differences between the two, though it generally supports the use of the Chihara formalism as an inexpensive alternative to the very detailed and computationally intensive TDDFT calculations. We anticipate that TDDFT may provide a powerful discriminating tool for arbitrating disagreements between these more phenomenological theories and experiment, especially as experimental data becomes more highly resolved. Our framework enables future explorations of systems in which the partition between bound and free electrons is more ambiguous. It also provides a platform for studying the impact of recent foundational developments in DFT at non-zero temperature [32,[38][39][40][41][42][43]. Supplemental Materials: X-ray Thomson scattering without the Chihara decomposition

IMPLEMENTATION DETAILS
We have implemented real-time TDDFT using a plane wave basis and the Projector Augmented Wave (PAW) formalism [23,24] within version 5.3.5 of VASP [25,26]. As in other implementations using ultrasoft pseudo-potentials [44] or PAW [45] we chose a Crank-Nicolson (CN) time integration scheme to propagate the time-dependent Kohn-Sham (TDKS) equations numerically rather than chosing an integrator to achieve a high-order convergence in the local error. The unitarity of the discrete propagator was deemed an essential factor to guarantee the satisfaction of charge conservation and other sum rules associated with the dynamic structure factor (DSF) itself. At each timestep, the CN-discretized TDKS equations are solved using the Generalized Minimal Residual method (GMRES) [46]. Given the perturbation on identity form of the discrete propagator we expect and observe rapid convergence in the number of iterations.
As many details of our implementation have not been published elsewhere, we begin by demonstrating that it is robust. We first test stability of the time integration. Specifically, we illustrate that our implementation is not susceptible to any significant or uncontrolled errors, given the intrinsic nonlinearity of the TDKS equations and the use of an iterative solver on the linearized problem. To do so, we time propagate a Mermin thermal equilibrium set of orbitals in the absence of a probe potential or ionic motion. Here, the time-dependent potential should be constant in time, and any variations in the instantaneous total free energy or supercell charge are due to the accumulation of numerical errors.
Results are reported for 32 beryllium atoms under the warm dense conditions in the paper and using the adiabatic zero-temperature local density approximation throughout. The initial condition was generated from a Mermin DFT calculation done on a single atomic configuration drawn from a thermalized DFT-MD run with a four-electron PAW potential, 576 thermally occupied orbitals, and a plane wave cutoff of 1400 eV. The time propagation utilized the same plane wave cutoff and number of orbitals, and a time step (∆t) of 1 attoseconds (as) carried out for 8000 steps (8 fs). We varied the relative tolerance of the iterative solver, and kept the absolute tolerance fixed to be 2 orders of magnitude smaller. As the right hand side of the CN-discretized TDKS system is of order unit norm, the absolute tolerance is effectively irrelevant. The resultant free energy per atom and total charge density of the supercell are reported in Figure S1.   Here, we see that the free energy drift per atom is −2.6 × 10 1 µeV/(atom·fs) and that the drift in the total charge of the unit cell is 8.7 × 10 −6 electrons/fs in the most permissive case (8 digits). In the least permissive case (14 digits), these quantities are 1 × 10 −3 µeV/(atom·fs) and 1.3 × 10 −11 electrons/fs, with the 12 digit case producing drifts of a similar order of magnitude. Going from 8 digits to 14 digits, the average CPU time per step varies linearly from 0.85 s/time step to 1.64 s/time step, i.e., convergence is exponential in the number of GMRES iterations. As the results are practically indistinguishable from 14 digits and the CPU time per time step is 20% shorter, a relative tolerance of 12 digits was used in all subsequent calculations. Pertinent to the satisfaction of sum rules on δρ(r, t), in all calculations in this work, we have verified that charge conservation is guaranteed to a relative accuracy of greater than 8 digits.
We next demonstrate convergence of the integrated density response, the observable of interest, in ∆t. To do so, we consider the same 32 beryllium atom initial condition described above, this time applying a time-varying scalar perturbation of the form described in the main body of the paper. Here with t d = 10 as, t w = 2 as, v 0 = 0.001 eV·fs, and q = 1.091xÅ −1 . ∆t is varied from from 4 as to 0.25 as, and the associated convergence of the real-time density response is illustrated in Figure S2. From these results, it is evident that the density response exhibits first order convergence. Further, the convergence supports the decision that ∆t = 1 as (attosecond) provides a reasonable balance between accuracy and time required per calculation for generating production results.
Next, we consider the determination of parameters for v probe (r, t). We chose f (t) to take the form of a Gaussian envelope to ensure that the exciting potential and its response are approximately band-and time-limited. To this end, the pulse width (t w ), determines the bandwidth of the response and was chosen to ensure that modes of the density response with energies on the order of 100s of eV are excited with appreciable amplitude. The delay, t d , was chosen to ensure that the excitation is approximately quiescent at t = 0 as. The remaining parameter, v 0 , determines the effective intensity of the probe potential, and must be chosen to be large enough that the response is not dominated by numerical noise, yet small enough to remain in the linear response regime. Varying v 0 over 4 orders of magnitude from 1 eV·fs to 0.001 eV·fs, we do not observe numerical noise to be a problem at any value. The post-processed DSF computed using these probe amplitudes for a single 32 atom configuration at |q| = 1.091Å −1 (θ = 20 • ) are illustrated in Figure S3. The results for v 0 ranging from 0.001 eV·fs to 0.1 eV·fs are indistinguishable, while the distortion of the results at 1 eV·fs indicate the onset of physics beyond linear response. That we can easily access this regime is one of the benefits of real-time TDDFT, though we do not explore this further in this work.

GENERATION OF INITIAL CONDITIONS
Each TDDFT calculation requires a set of Kohn-Sham orbitals and occupancies as initial conditions. These are generated using DFT-MD as implemented in the standard version of the VASP software package [24][25][26]. To assess the impact of the supercell shape and the underlying ionic positions on our DSF calculations, for each value of q we run 2 separately thermalized DFT-MD trajectories on different supercells and draw 5 independent sets of initial conditions from each. Each set of initial conditions is used to seed independent TDDFT calculations of the DSF at a fixed q.
Separate DFT-MD configurations are needed for each value of q due to the requirements of the probe potential. For the DSF calculations, v probe (r, t) must be commensurate with our supercell, and consequently any realizable value of q must be in the reciprocal lattice of our supercell. To precisely specify the value of q we work with tetragonal supercells in which the perturbing q is directed along the c-axis. Then q can be set by scaling the c-axis and the a-axes can be adjusted to ensure that the desired mass density is realized. Given the liquid-like ordering in the extreme conditions under investigation, we do not anticipate that this biases our results. Table I gives the dimensions of the 8 supercells used in this study. It is worth noting that all DFT-MD trajectories are generated with the full four-electron PAW potential. However, both the two-electron and four-electron TDDFT calculations are initialized from this same set of ionic configurations, with the Kohn-Sham orbitals and occupancies being recomputed for each set of fixed ionic positions in the two-electron case. This is done to assess the impact of the frozen core approximation on the DSF in isolation.
All DSF results reported are averaged over 10 configurations sampled at each value of q. The sample size is necessarily small due to the significant computational resources required for each TDDFT calculation, with each production calculation being done on 1,152 cores. However, because the electronic density of warm dense beryllium is relatively uniform, we do not see much variability in the DSF from configuration to configuration. In an effort to quantify this variability, we apply jackknife resampling to estimate the variance and indicate single standard deviation intervals. For the scattering angles considered in the manuscript, we show the estimated error intervals for the DSF computed with four-electron PAW potentials in Figure S4. ] IG. S4. Jackknife error estimates of the DSF for warm dense beryllium at 4 different scattering angles. The colored shaded regions bound ±4σ about the mean, which is a solid black line. The ±1σ intervals are barely perceptible, and we simply present the mean values in the manuscript.

SATISFACTION OF SUM RULES
The units on the density response and DSF are such that the following form of the f-sum rule [47] is satisfied: Here, N e is 2 for the two-electron PAW, and 4 for the four-electron PAW. Similar real-time calculations in the gas phase have reported satisfaction to within 5% [33] and we report similar results. Forms of the DSF based upon model dielectric functions may or may not be consistent with various sum rules by construction. Our DSF is derived strictly from a time-evolved electronic density for which charge conservation is numerically guaranteed to high precision, and makes no assumptions about the form of the equivalent dielectric response. To this end, the primary source of error in our satisfaction of sum rules is due to numerical errors in the post-processing, e.g., the numerical evaluation of Fourier integrals with integrands which become increasingly oscillatory at higher energies, and thus more prone to small errors in the response. When checking sum rules produced using linear response TDDFT, it is common practice to fit the high energy tail of the DSF to force it to go to zero in a way that is consistent with the integrand in the above sum rule [48]. This is also critical for real-time simulations where small phase errors between the real and imaginary parts of the time domain density response are amplified at high energies relative to low energies when post-processing to generate the energy domain response. Rather than force the tail to fit some form, we simply impose a high energy cutoff on the integrand once the DSF has decayed to < 1% of its peak value. We simply seek to verify that the majority of the weight of our response is consistent with the sum rule, and could resort to fitting tails if necessary to improve agreement. Applying this technique, we verify that we satisfy the f-sum rule for our four-electron data with a relative error of −7% (θ = 20 • ), −3% (θ = 40 • ), −2% (θ = 60 • ), and −4% (θ = 90 • ). In all cases, we underestimate the value of the f-sum rule, giving us confidence that a fit of the slowly-decaying high energy tail might be used to improve this agreement. Applying this same technique to the two-electron data we find relative errors of 8% (θ = 20 • ), 5% (θ = 40 • ), 0.4% (θ = 60 • ), and −2% (θ = 90 • ). Here, we do not uniformly underestimate the sum rule as was the case with the four-electron data. The two-electron data decays much more abruptly at high energies, and does not need to represent the response of slowly-decaying core states, so we do not view this as a deficiency (i.e., fitting a tail would not have as strong of an impact here). However, this seems to point to the two-electron response as overestimating the free-free peak, consistent with the small negative dips in the bound-free data in Fig. 3b for θ = 20 • , 40 • , and 60 • (the calculations for which the sum rule was over-estimated). Whether or not this can be improved with a different two-electron PAW may be an interesting topic of further study.
Perhaps a more interesting sum rule to study is the one that defines the static structure factor through the integral of the DSF: Computing this integral for the effective free-free and bound-free data from TDDFT, we can extract structure factors that we can compare against physically intuitive models. Results are presented in Fig. S5. The static structure factor computed using Eqn. S2 applied to the effective free-free and bound-free responses extracted from TDDFT. The bound-free structure factor is compared to that of doubly-and triply-ionized beryllium [49].. The free-free structure factor is compared to an analytic fit of QMC data for jellium with rs = 1.3 a.u [51].
For the bound-free component we compare to results obtained using a configuration-interaction expansion by Brown for doubly-and triply-ionized beryllium [49]. For θ = 60 • and 90 • , the TDDFT gives results that are closer to doublyionized beryllium, consistent with physical intuition for the thermodynamic conditions under consideration. For θ = 20 • and 40 • , the TDDFT gives results that are closer to triply-ionized beryllium, though the difference between the doubly-and triply-ionized structure factors are smaller at these angles. We note that by excluding the negative difference data below 80 eV in our integration (evident in Fig. 3b), we can improve agreement with the doubly ionized curve at all angles. This indicates that the underestimation of the TDDFT structure factor may be due to differences between the free electron response for the two-electron and four-electron PAWs, which may or may not be physical. This highlights the importance of being able to self-consistently compute the four-electron response without the Chihara decomposition using our methodology.
For the free-free component (our frozen core result) we compare to results for the static structure factor of jellium with r s = 1.3 a.u., consistent with the experimentally determined free-electron density [10]. In evaluating Eqn. S2, we remove the elastic peak in a region of width ∼ k B T e centered at ω = 0 and apply a simple linear interpolant between the DSF data on either side of the excluded region. We compare the resultant free-free static structure factor to an analytic fit to QMC data [53] by Gori-Giorgi, et. al. [51]. The TDDFT free-free structure factor exhibits the same trend as the QMC fit, but is not expected to agree perfectly. We only expect qualitative agreement because the external potential experienced by the free electrons in warm dense beryllium is not a uniform neutralizing background as is the case in jellium.
This qualitative agreement between results from TDDFT and QMC stands in contrast to Vorberger and Gericke's recent work in which DFT-MD is used to compute a free electron density, from which the free-free static structure factor of warm dense beryllium is extracted [20]. In their work, the static structure factor computed from DFT differs greatly from QMC, namely it does not go from one to zero for momentum transfers less than 2k F (k F = 2.88Å −1 for our conditions), but instead oscillates slightly about unity. The authors postulate that this is due to the mean field nature of the Kohn-Sham equations, and that this may be beyond the capability of DFT. Our results indicate that this basic physics is in fact within the grasp of TDDFT. To this end, it is worth noting that TDDFT has provided us with a convenient means of computing the static structure factor as an exact functional of the time-dependent density. In this case, we mean exact in the sense that if we are given a representation of the exact interacting time-dependent density, we can map it onto the exact DSF and the exact structure factor through Eqn. S2.