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. 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 scienceGuillot (1999); Fortney and Nettelmann (2009); Kraus et al. (2010) to the implosion stage of inertial confinement fusionMatzen et al. (2005); Drake (2009); Hu et al. (2014). 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 magnitudeGlenzer and Redmer (2009).
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)Gregori et al. (2003); Glenzer and Redmer (2009), 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 berylliumGlenzer et al. (2007); Lee et al. (2009), lithiumSaiz et al. (2008), carbonBrown et al. (2014), CH shellsFletcher et al. (2014), and aluminumMa et al. (2013); Fletcher et al. (2015). With recent improvements in source brightnessFletcher et al. (2015) 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 DSF varies with momentum and energy transfers ( 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 (), from free electrons per ion (), and from bound electrons that are photo-ionized (). 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 dataSouza et al. (2014). The ionic feature typically relies upon a decomposition into a product of an ion-ion structure factor, , and an average atomic form factor, , with the first term describing the unscreened bound electrons and the latter the screening cloud, which must be treated carefullyPlagemann et al. (2015). However, recent work has focused on moving past this decomposition of the elastic peakVorberger and Gericke (2015).
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 time-dependent density functional theory (TDDFT)Ullrich (2011); Marques et al. (2012) and the Projector Augmented-Wave (PAW) formalism Blöchl (1994); Kresse and Joubert (1999). 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)Kresse and Furthmüller (1996a, b); Kresse and Joubert (1999) (see Supplemental Material 111See Supplemental Material, which includes Refs. Qian et al. (2006); Ojanperä et al. (2012); Saad (2003); Mahan (2000); Weissker et al. (2010); Brown (1970); Perdew and Wang (1992); Gori-Giorgi et al. (2000); Ceperley and Alder (1980); Ortiz et al. (1999)) 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 for the calculation of the initial state of the Kohn-Sham orbitals. This leads to excellent strong scalingAndrade et al. (2012); Schleife et al. (2014), 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 () unless otherwise indicated. Time-dependent quantities evaluated at are indicated by the addition of a subscript 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:
in which . The orbitals, indexed by band and Bloch wave number, are such that their weighted sum produces the time-dependent density, . The external potential, , includes contributions due to the Coulomb field of the bare nuclei as well as a model of the x-ray probe, , which is quiescent until . and 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 from Eqn. 2 are the self-consistent solution to a Kohn-Sham Mermin DFT calculationMermin (1965) at electron temperature in a supercell of volume . The initial density is then:
where 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, . 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 regimeGiuliani and Vignale (2005). 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 workPribram-Jones et al. (2015).
The action of a probe potential, , turned on at leads to a change in the time-dependent density. The linear density-density response function, , encodes the relationship between these two quantities:
where we use the notation for and . If we fix ionic positions at time , which evolve slowly on the relevant attosecond timescales, then and we can construct a probe potential that can be used to extract the DSF, similar to Sakko, et. al.Sakko et al. (2010). The real-time density response to such a probe potential is shown in Fig. 1.
In principle, any sufficiently weak analytic probe potential will allow us to extract the response function. For convenience, we choose , where is a Gaussian envelope and is related to the probe intensity. The Fourier transformed response function, , is then related to the DSF through the fluctuation-dissipation theorem:
We are careful to note that Fourier transforms are normalized such that 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 Lee et al. (2009). We consider the same momentum transfers as in Plagemann et al. (2012) to study a range of excitations spanning the collective regime to the beginning of non-collective regime. For convenience of presentation, each value is mapped onto an XRTS scattering angle, , relative to the 2 Å probe wavelength in Lee et al. (2009) (see Supplemental Material Note1 () 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 eV, and -point sampling at , analogous to the Baldereschi mean-value point for cubic supercells. For these conditions a plane wave cutoff of eV was required to converge the pressure to within , and 576 (32 atom)/1152 (64 atom) Kohn-Sham orbitals were needed to represent the thermal occupation ().
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 -point sampling on a (32 atom) or (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 Note1 ().
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 . 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 and 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 real-time 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 , we can also extract a quantity akin to by differencing the four-electron and two-electron DSFs. The effective and computed within TDDFT are illustrated in Fig. 3. Here we compare our TDDFT calculations to calculations done using state-of-the-art models for and . The former is treated with an RPA-level 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 Plagemann et al. (2012). The latter is calculated with the formalism developed in Johnson et al. (2012) 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 , though it is necessarily present in the Chihara-independent 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 and , whereas it predicts a stronger (10 eV) red shift relative to MA-AICF at and . 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 RPACazzaniga et al. (2011). 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 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 with some minor differences. Such average atom models agree well with real-space Green’s function methods for cold solid beryllium Mattern et al. (2012). 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 functional 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-of-the-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 temperatureLi and Tong (1985); Li and Li (1985); Brown et al. (2013); Karasiev et al. (2014); Pribram-Jones et al. (2014a, b, 2015).
Calculations were completed on Chama and Skybridge at Sandia National Laboratories and Sequoia at Lawrence Livermore National Laboratory. We are grateful for discussions with Stephen Bond, Kieron Burke, Hardy Gross, Ryan Hatcher, Neepa Maitra, Thomas Mattsson, Normand Modine, Jonathan Moussa, Kai-Uwe Plagemann, Aurora Pribram-Jones, Kenneth Rudinger, Travis Sjostrom, and Sam Trickey. A.D.B, L.S., M.P.D., and R.J.M. were supported by Sandia’s Laboratory Directed Research and Development (LDRD) Project 165731. S.B.H. was supported by the U.S. Department of Energy, Office of Science Early Career Research Program, Office of Fusion Energy Sciences. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.
- Guillot (1999) T. Guillot, Science 286, 72 (1999).
- Fortney and Nettelmann (2009) J. J. Fortney and N. Nettelmann, Space Science Reviews 152, 423 (2009).
- Kraus et al. (2010) R. Kraus, S. Stewart, A. Seifter, and A. Obst, Earth and Planetary Science Letters 289, 162 (2010).
- Matzen et al. (2005) M. K. Matzen, M. A. Sweeney, R. G. Adams, J. R. Asay, J. E. Bailey, G. R. Bennett, D. E. Bliss, D. D. Bloomquist, T. A. Brunner, R. B. Campbell, G. A. Chandler, C. A. Coverdale, M. E. Cuneo, J.-P. Davis, C. Deeney, M. P. Desjarlais, G. L. Donovan, C. J. Garasi, T. A. Haill, C. A. Hall, D. L. Hanson, M. J. Hurst, B. Jones, M. D. Knudson, R. J. Leeper, R. W. Lemke, M. G. Mazarakis, D. H. McDaniel, T. A. Mehlhorn, T. J. Nash, C. L. Olson, J. L. Porter, P. K. Rambo, S. E. Rosenthal, G. A. Rochau, L. E. Ruggles, C. L. Ruiz, T. W. L. Sanford, J. F. Seamen, D. B. Sinars, S. A. Slutz, I. C. Smith, K. W. Struve, W. A. Stygar, R. A. Vesey, E. A. Weinbrecht, D. F. Wenger, and E. P. Yu, Physics of Plasmas 12, 55503 (2005).
- Drake (2009) R. P. Drake, Physics of Plasmas 16, 55501 (2009).
- Hu et al. (2014) S. X. Hu, L. A. Collins, T. R. Boehly, J. D. Kress, V. N. Goncharov, and S. Skupsky, Physical Review E 89, 043105 (2014).
- Glenzer and Redmer (2009) S. Glenzer and R. Redmer, Reviews of Modern Physics 81, 1625 (2009).
- Gregori et al. (2003) G. Gregori, S. H. Glenzer, W. Rozmus, R. W. Lee, and O. L. Landen, Physical Review E 67, 026412 (2003).
- Glenzer et al. (2007) S. .H. Glenzer, O. .L. Landen, P. Neumayer, R. W. Lee, K. Widmann, S. W. Pollaine, R. J. Wallace, G. Gregori, A. Höll, T. Bornath, R. Thiele, V. Schwarz, W.-D. Kraeft, and R. Redmer, Physical Review Letters 98, 065002 (2007).
- Lee et al. (2009) H. J. Lee, P. Neumayer, J. Castor, T. Döppner, R. W. Falcone, C. Fortmann, B. A. Hammel, A. L. Kritcher, O. L. Landen, R. W. Lee, D. D. Meyerhofer, D. H. Munro, R. Redmer, S. P. Regan, S. Weber, and S. H. Glenzer, Physical Review Letters 102, 115001 (2009).
- Saiz et al. (2008) E. G. Saiz, G. Gregori, D. O. Gericke, J. Vorberger, B. Barbrel, R. Clarke, R. R. Freeman, S. Glenzer, F. Khattak, M. Koenig, et al., Nature Physics 4, 940 (2008).
- Brown et al. (2014) C. Brown, D. O. Gericke, M. Cammarata, B. Cho, T. Döppner, K. Engelhorn, E. Förster, C. Fortmann, D. Fritz, E. Galtier, et al., Scientific reports 4 (2014).
- Fletcher et al. (2014) L. Fletcher, A. Kritcher, A. Pak, T. Ma, T. Döppner, C. Fortmann, L. Divol, O. Jones, O. Landen, H. Scott, et al., Physical review letters 112, 145004 (2014).
- Ma et al. (2013) T. Ma, T. Döppner, R. Falcone, L. Fletcher, C. Fortmann, D. O. Gericke, O. Landen, H. Lee, A. Pak, J. Vorberger, et al., Physical review letters 110, 065001 (2013).
- Fletcher et al. (2015) L. Fletcher, H. Lee, T. Döppner, E. Galtier, B. Nagler, P. Heimann, C. Fortmann, S. LePape, T. Ma, M. Millot, et al., Nature Photonics 9, 274 (2015).
- Chihara (1987) J. Chihara, Journal of Physics F: Metal Physics 17, 295 (1987).
- Chihara (2000) J. Chihara, Journal of Physics: Condensed Matter 12, 231 (2000).
- Souza et al. (2014) A. N. Souza, D. J. Perkins, C. E. Starrett, D. Saumon, and S. B. Hansen, Physical Review E 89, 023108 (2014).
- Plagemann et al. (2015) K.-U. Plagemann, H. R. Rüter, T. Bornath, M. Shihab, M. P. Desjarlais, C. Fortmann, S. H. Glenzer, and R. Redmer, Physical Review E 92, 013103 (2015).
- Vorberger and Gericke (2015) J. Vorberger and D.O. Gericke, Physical Review E 91, 033112 (2015).
- Ullrich (2011) C. A. Ullrich, Time-dependent density-functional theory: concepts and applications (Oxford University Press, 2011).
- Marques et al. (2012) M. A. Marques, N. T. Maitra, F. Nogueira, E. Gross, and A. Rubio, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 837 (2012).
- Blöchl (1994) P. E. Blöchl, Physical Review B 50, 17953 (1994).
- Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
- Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996a).
- Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996b).
- (27) See Supplemental Material, which includes Refs. Qian et al. (2006); Ojanperä et al. (2012); Saad (2003); Mahan (2000); Weissker et al. (2010); Brown (1970); Perdew and Wang (1992); Gori-Giorgi et al. (2000); Ceperley and Alder (1980); Ortiz et al. (1999).
- Andrade et al. (2012) X. Andrade, J. Alberdi-Rodriguez, D. A. Strubbe, M. J. Oliveira, F. Nogueira, A. Castro, J. Muguerza, A. Arruabarrena, S. G. Louie, A. Aspuru-Guzik, et al., Journal of Physics: Condensed Matter 24, 233202 (2012).
- Schleife et al. (2014) A. Schleife, E. W. Draeger, V. M. Anisimov, A. A. Correa, and Y. Kanai, Computing in Science & Engineering 16, 54 (2014).
- Mermin (1965) N. Mermin, Physical Review 137, A1441 (1965).
- Giuliani and Vignale (2005) G. Giuliani and G. Vignale, Quantum theory of the electron liquid (Cambridge University Press, 2005).
- Pribram-Jones et al. (2015) A. Pribram-Jones, P. E. Grabowski, and K. Burke, arXiv preprint arXiv:1509.03071 (2015).
- Sakko et al. (2010) A. Sakko, A. Rubio, M. Hakala, and K. Hämäläinen, The Journal of chemical physics 133, 174111 (2010).
- Plagemann et al. (2012) K.-U. Plagemann, P. Sperling, R. Thiele, M. P. Desjarlais, C. Fortmann, T. Döppner, H. J. Lee, S. H. Glenzer, and R. Redmer, New Journal of Physics 14, 55020 (2012).
- Johnson et al. (2012) W. R. Johnson, J. Nilsen, K. T. Cheng, Physical Review E 86, 036410 (2012).
- Cazzaniga et al. (2011) M. Cazzaniga, H.-C. Weissker, S. Huotari, T. Pylkk, P. Salvestrini, G. Onida, and L. Reining, 075109, 1 (2011).
- Mattern et al. (2012) B. A. Mattern, G. T. Seidler, J. J. Kas, J. I. Pacold, and J. J. Rehr, Physical Review B 85, 115135 (2012).
- Li and Tong (1985) T.-C. Li and P.-Q. Tong, Physical Review A 31, 1950 (1985).
- Li and Li (1985) T. Li and Y. Li, Physical Review A 31, 3970 (1985).
- Brown et al. (2013) E. W. Brown, J. L. DuBois, M. Holzmann, and D. M. Ceperley, Physical Review B 88, 081102 (2013).
- Karasiev et al. (2014) V. V. Karasiev, T. Sjostrom, J. Dufty, and S. B. Trickey, Physical Review Letters 112, 076403 (2014).
- Pribram-Jones et al. (2014a) A. Pribram-Jones, Z.-H. Yang, J. R. Trail, K. Burke, R. J. Needs, and C. A. Ullrich, The Journal of chemical physics 140, 18A541 (2014a).
- Pribram-Jones et al. (2014b) A. Pribram-Jones, S. Pittalis, E. K. U. Gross, and K. Burke, in Frontiers and Challenges in Warm Dense Matter, Lecture Notes in Computational Science and Engineering, Vol. 96, edited by F. Graziani, M. P. Desjarlais, R. Redmer, and S. B. Trickey (Springer International Publishing, 2014) pp. 25–60.
- Qian et al. (2006) X. Qian, J. Li, X. Lin, and S. Yip, Physical Review B 73, 035408 (2006).
- Ojanperä et al. (2012) A. Ojanperä, V. Havu, L. Lehtovaara, and M. Puska, The Journal of chemical physics 136, 144103 (2012).
- Saad (2003) Y. Saad, Iterative methods for sparse linear systems (Siam, 2003).
- Mahan (2000) G. D. Mahan, Many-Particle Physics (Springer Science & Business Media, 2000).
- Weissker et al. (2010) H.-C. Weissker, J. Serrano, S. Huotari, E. Luppi, M. Cazzaniga, F. Bruneval, F. Sottile, G. Monaco, V. Olevano, and L. Reining, Physical Review B 81, 085104 (2010).
- Brown (1970) R. T. Brown, Physical Review A 2, 614 (1970).
- Perdew and Wang (1992) J. P. Perdew and Y. Wang, Physical Review B 46, 12947 (1992).
- Gori-Giorgi et al. (2000) P. Gori-Giorgi, F. Sacchetti, and G. B. Bachelet, Physical Review B 61, 7353 (2000).
- Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, Physical Review Letters 45, 566 (1980).
- Ortiz et al. (1999) G. Ortiz, M. Harris, and P. Ballone, Physical Review Letters 82, 5317 (1999).
Supplemental Materials: X-ray Thomson scattering without the Chihara decomposition
Appendix A Implementation Details
We have implemented real-time TDDFT using a plane wave basis and the Projector Augmented Wave (PAW) formalism Blöchl (1994); Kresse and Joubert (1999) within version 5.3.5 of VASP Kresse and Furthmüller (1996a, b). As in other implementations using ultrasoft pseudo-potentials Qian et al. (2006) or PAW Ojanperä et al. (2012) 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) Saad (2003). 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 () 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 eV/(atomfs) and that the drift in the total charge of the unit cell is electrons/fs in the most permissive case (8 digits). In the least permissive case (14 digits), these quantities are eV/(atomfs) and 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 shorter, a relative tolerance of 12 digits was used in all subsequent calculations. Pertinent to the satisfaction of sum rules on , 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 . 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 where with as, as, eVfs, and Å. is varied from from as to 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 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 . We chose 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 (), determines the bandwidth of the response and was chosen to ensure that modes of the density response with energies on the order of s of eV are excited with appreciable amplitude. The delay, , was chosen to ensure that the excitation is approximately quiescent at as. The remaining parameter, , 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 over 4 orders of magnitude from eVfs to eVfs, 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 Å are illustrated in Figure S3. The results for ranging from eVfs to eVfs are indistinguishable, while the distortion of the results at eVfs 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.
Appendix B 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 Kresse and Furthmüller (1996b, a); Kresse and Joubert (1999). To assess the impact of the supercell shape and the underlying ionic positions on our DSF calculations, for each value of 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 .
Separate DFT-MD configurations are needed for each value of due to the requirements of the probe potential. For the DSF calculations, must be commensurate with our supercell, and consequently any realizable value of must be in the reciprocal lattice of our supercell. To precisely specify the value of we work with tetragonal supercells in which the perturbing is directed along the c-axis. Then 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 1 gives the dimensions of the 8 supercells used in this study.
|a (Å)||c (Å)||# of atoms||# of bands||(Å)|
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 . 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.
Appendix C Satisfaction of Sum Rules
The units on the density response and DSF are such that the following form of the f-sum rule Mahan (2000) is satisfied:
Here, 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 Sakko et al. (2010) 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 Weissker et al. (2010). 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 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 (), (), (), and (). 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 (), (), (), and (). 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 , , and (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.
For the bound-free component we compare to results obtained using a configuration-interaction expansion by Brown for doubly- and triply-ionized beryllium Brown (1970). For and , the TDDFT gives results that are closer to doubly-ionized beryllium, consistent with physical intuition for the thermodynamic conditions under consideration. For and , 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 a.u., consistent with the experimentally determined free-electron density Lee et al. (2009). In evaluating Eqn. S2, we remove the elastic peak in a region of width centered at 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 Ortiz et al. (1999) by Gori-Giorgi, et. al. Gori-Giorgi et al. (2000). 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 Vorberger and Gericke (2015). 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 ( Å 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.