Thermodynamic and transport coefficients from the dynamic structure factor of Yukawa liquids

The ion-ion dynamic structure factor (DSF) of warm dense matter and dense plasmas contains information on collective ionic modes and various thermodynamic and transport coefficients, which are important for modeling the interiors of giant planets or the dense plasmas occurring in inertial confinement fusion. Here, it is demonstrated, using the Yukawa liquid as a reduced model, that the complete hydrodynamic information encoded in the DSF can be extracted with an accuracy comparable to that of dedicated methods. This is achieved by applying a generalized hydrodynamic and a viscoelastic model and extrapolating the results at finite wave numbers into the hydrodynamic limit. Very good agreement with previous data is obtained for the sound speed and the viscosity. The thermal diffusivities deduced from different methods exhibit somewhat larger deviations.


I. INTRODUCTION
Warm dense matter (WDM) is an unusual state of matterbetween a solid and a fully ionized plasma-with solid densities, megabar pressures, and temperatures in the eV range [1]. At keV temperatures, one enters the regime of a hot dense plasma [2]. Matter at such extreme conditions is characterized by partial ionization, moderately coupled and partially degenerate electrons, and strongly coupled classical ions. The thermodynamic and transport properties of WDM and dense plasmas are of high importance for a variety of fields, including planetary science [3][4][5] and astrophysics [6,7], and for the advancement of inertial confinement fusion concepts [8][9][10].
Various methods have been devised for the creation and diagnostic of these exotic states [11]. In particular, x-ray scattering provides detailed insights into their properties [12][13][14][15][16][17]. The analysis of collective modes of the electrons [18], e.g., gives access to the electron density and temperature [19,20], and the electrical conductivity [21]. A significant step forward can be expected from the LCLS at Stanford [22] or the European XFEL [23]. In particular, it should become possible to diagnose collective modes of the ions [24], which carry information on fundamental material properties such as the viscosity but are very difficult to resolve. Within the Chihara decomposition [25] of the total electron dynamic structure factor, which determines the scattering signal [13], they manifest themselves in the ion-ion dynamic structure factor (DSF).
In recent years, significant efforts have been made to compute the ion-ion DSF from ab initio simulations [26], combining a classical description of the ions with a quantum treatment of the electrons [27][28][29][30][31]. The computational cost can be reduced by employing effective ion-ion pair potentials [32][33][34]. Several theoretical models have been proposed to describe the DSF in dense plasmas and WDM. Mithen et al. [35] have used the hydrodynamic result to model the DSF of a strongly coupled Yukawa liquid and later extended their investigations via the memory function formalism [36,37]. Vorberger et al. [38] performed a comparison of various models [35,[39][40][41] for the DSF of shocked Silicon, computed from simulations with a modified Yukawa potential. Recently, Choi et al. [42] performed a similar comparison for Coulomb and Yukawa liquids. Arkhipov et al. [43] used the moments method to investigate the DSF of Coulomb and Yukawa one-component plasmas.
Rather than attempting to reproduce the DSF with as few free parameters as possible, a different goal is pursued in this work. In the hydrodynamic limit, the DSF contains information on the fundamental thermodynamic and transport properties of a liquid [44][45][46]. Here, it is demonstrated that this can be exploited to extract various hydrodynamic properties of strongly coupled Yukawa liquids solely from the DSF. The Yukawa model is chosen for two reasons. First, it has been used as a simplified model for a reduced description of WDM and dense plasmas [2,34,47]. For a direct comparison with experimental or ab initio simulation data, a short-ranged repulsive potential can be added to the Yukawa potential to capture effects related to bound electrons and to better match the ion structural properties [20,38,48]. Second, and more importantly, a large amount of accurate data are available for the thermodynamic [49][50][51] and transport properties [52][53][54][55][56], which will be used to demonstrate the feasibility of the method. The methodology could provide a practical means for the determination of hydrodynamic coefficients (e.g., viscosity) from experimental and simulation data for the ion-ion DSF of WDM and dense plasmas, which are of high importance for simulations of inertial confinement fusion implosions or the modeling of planetary interiors. Since the Yukawa liquid is a very general model system, the results could also be relevant to ultracold neutral plasmas [57], dusty (complex) plasmas [58], or colloidal dispersions [59]. In the latter two cases, however, frictional damping is important.
The analysis of Mithen et al. [35] has shown that the Rayleigh-Brillouin triplet [44,46] accurately describes the DSF of Yukawa liquids at small wave numbers k. It consists of a central Lorentzian line, representing a diffusive heat mode, and two shifted Lorentzians, representing propagating sound modes. The hydrodynamic parameters determine their widths and intensities. When fitted to the DSF at a single wave number [35,37], the model gave hydrodynamics coefficients that were in reasonable agreement with the available data at the time. In this work, the DSF is analyzed using two extensions of the hydrodynamic model. In the generalized hydrodynamic (GH) model [60], the frequency dependence of the hydrodynamic DSF is left unchanged but the peak intensities and widths can have an arbitrary k dependence. The viscoelastic (VE) model extends hydrodynamics to finite frequencies by introducing the concept of viscoelasticity. Here, the DSF contains an additional central Lorentzian peak, which vanishes in the hydrodynamic limit [60]. By (i) fitting the model parameters to simulation data and (ii) extrapolating the results to k → 0 (similar to recent work by Guarini et al. [61] and Silvestri et al. [62]), the complete thermodynamic and transport properties contained in the DSF are obtained, with an accuracy comparable to dedicated other methods [52][53][54][55][56]. The VE model allows one to go beyond a purely hydrodynamic description and yields results for the viscoelastic properties.
This paper is organized as follows. A brief description of the simulation method is given in Sec. II. The two models used for the DSF and their relation to the hydrodynamic limit are introduced in Sec. III. Results for the DSF and the hydrodynamic parameters are presented and discussed in Sec. IV. Section V provides a brief summary and outlook.

II. SIMULATIONS
Molecular dynamics (MD) simulations are carried out for particles of mass m interacting via the Yukawa potential, v(r) = q 2 exp(−r/λ s )/r, where q is the particle charge and λ s the screening length. The screening parameter, κ = a/λ s , where a = (3/4π n) 1/3 is the Wigner-Seitz radius and n the particle density, is chosen as κ = 2. The reason is that a large body of data for the transport coefficients are available [52][53][54][55][56]. The values for the coupling parameter, = q 2 /(a k B T ), where T is the temperature, are in the range = 10, . . . , 200, covering a wide range of conditions -from a moderately coupled system at / m ≈ 0.02 up to a strongly coupled liquid at / m ≈ 0.5, where m ≈ 440 is the melting point [49]. Two sets of simulations were carried out using N = 3800 and N = 5000 particles, which gives access to wave numbers with ka 0.23. 1 The dynamic structure factor, which is defined in terms of the intermediate scattering is computed from the Fourier transform of n k (t ) = N i=1 e −ik·r i (t ) [63].

III. THEORY
As discussed above, the hydrodynamic result for the DSF can be written as the sum of a central Lorentzian (will be referred to as "heat mode") and two shifted Lorentzians ("sound modes") [44,46]. Compared to some of the previous applications in the field of WDM [32,35,38,64], the asymmetry of the latter [45] is taken into account, which leads to a finite second frequency moment [60]. Following Bafile et al. [60], the generalized hydrodynamics (GH) model is written as where S(k) is the static structure factor, I h [I s = (1 − I h )/2] denotes the intensity of the heat (sound) peak and z h (z s ) its ] describes the asymmetry of the sound peaks. In the hydrodynamic limit, the parameters are given by, to lowest order in k, Here, D T is the thermal diffusivity, s the sound attenuation coefficient, γ = c p /c v the adiabatic index (with specific heats at constant pressure and volume, c p and c v , respectively), ν l = (4 η s /3 + η b )/(nm) the longitudinal kinematic viscosity (with the shear viscosity η s and the bulk viscosity η b ), and c s the adiabatic sound speed [45,60].
The VE model is typically formulated in the framework of the memory function formalism [44][45][46], where the instantaneous decay of the viscous term in the hydrodynamic model is replaced by an exponential decay [60]. It has been shown [60] that it can be expressed, at small k, in the same form as the GH model but with an additional central Lorentzian (intensity The intensity I 2 can be written in terms of z h , z s , z 2 , ω s , and I h as In the limit k → 0, I 2 should vanish as I 2 ∝ k 4 while the width z 2 should turn into a constant, z 2 → 1/τ , to be consistent with the hydrodynamic limit. Here, τ is the k → 0 limit of the viscoelastic relaxation time. Theory shows that for consistency with hydrodynamics, τ is given by τ = ν l /(c 2 L − c 2 s ), where c L is the k = 0 limit of the infinite-frequency (longitudinal) sound speed [60]. A liquid shows elastic response for frequencies ω τ −1 while viscous behavior is found at low frequencies, ω τ −1 . The widths and intensities of the heat and sound peak have the same limits as above. Compared to the GH model, the VE model has a finite fourth moment [60]. Both models can be considered few-term approximations of an exact expansion of the DSF in a series of Lorentzian functions [65]. The parameters of the memory function can be reconstructed from the above parameters, if necessary. In particular, Bafile et al. [60] provide explicit formulas for ω 2 k (k) = ω 2 /S(k) and ω 2 L (k) = ω 4 / ω 2 (only VE), where ω k are the frequency moments of the DSF.
A least-squares fit is used to determine the widths z h , z s , z 2 , the frequency of the sound mode, ω s , and the intensity I h . In addition, the static structure factor is used as a fit parameter. All other parameters follow from the model. Since the DSF decays very rapidly at high frequencies [42,66], the fit region for Eq. (2) is restricted to the central peak and the main body of the sound peak [ω ω s (k) + z s (k)]. The VE model is applicable at higher frequencies, and a larger fit range has been chosen, ω < 0.4 ω p , where ω p = 4π q 2 n/m is the plasma frequency.

A. Dynamic structure factor
A comparison between the simulation data and the GH fit is shown in Fig. 1. 2 While the model provides a very good fit function for the low-frequency part, the DSF decays faster at high frequencies. Also shown are the individual contributions of the heat and sound peaks from Eq. (2). At ka = 0.23, the central peak is almost fully reproduced by the heat mode while at ka = 0.68 the sound mode makes an important contribution at ω = 0. A similar comparison is shown in Fig. 2 for the VE model. One observes that the VE model accurately describes the decay of the DSF immediately following the sound peak. Figure 2(d) shows the intensities of the various peaks. The intensity of the second central peak, I 2 , is consistent with a k 4 decay in the k → 0 limit, as anticipated from theoretical considerations [60] and discussed above.

B. Fit parameters and extrapolation
The fit parameters are shown in Fig. 3. Their k → 0 limits yield values for (a) the sound speed, (c) the sound attenuation coefficient, (d) the inverse of the viscoelastic relaxation time, and (e) the thermal diffusivity. Figure 3(b) shows I h /(2 I s ), which reduces to the Landau-Placzek ratio, I h /(2 I s ) → γ − 1, in the hydrodynamic limit [46]. The hydrodynamic scaling, ω s → c s k, z s → s k 2 , and z h → D T k 2 , has been removed. For the isotropic liquid, a simple polynomial ansatz with even powers ofk = ka, f (k) = c 0 + c 2k 2 + c 4k 4 , is used for the estimation of the hydrodynamic coefficient, c 0 . In some cases, a simple average over the first few values has been employed instead. In particular for intermediate coupling strengths, 30 70, the hydrodynamic scaling is excellent, and the remaining dependence on k is typically weak, see 2 Note that the number of data points has been reduced for a better presentation. = 50. The thermal diffusivity is the most difficult quantity to determine (e.g., = 10 and 200), which is possibly related to the smallness of the central peak, in particular at strong coupling. Interestingly, the sound dispersion, ω s (k), appears to change from negative to positive for 10 and 200. However, larger simulations would be required to confirm this behavior.

C. Thermodynamic and transport data
The thermodynamic coefficients obtained from the extrapolation are shown in Fig. 4. For comparison, recent data for the sound speed by Silvestri et al. [62] and data obtained from fits for the internal energy of a single-component Yukawa fluid (without background terms) from Khrapak and Thomas [50,51] are included. Note that the sound speed and the adiabatic index are related by c 2 s = γ /(n m χ T ), where χ T is the isothermal compressibility [44,46]. Excellent agreement is found between all data sets for the sound speed (deviations 1%). The adiabatic index from the present analysis is slightly smaller than the analytical result for 20. For larger , the agreement is very good.  [50,51] and the MD results of Silvestri et al. [62] are shown for comparison. Figure 5 depicts the results for the transport coefficients. Recent data for the thermal diffusivity D T are available from the work of Ott et al. [56]. In addition, D T can be computed from the thermal conductivity λ and c p via D T = λ/(n c p ) [45,60]. While results for λ are available from Refs. [52,55], c p is calculated from the fits of Refs. [50,51]. The results for D T obtained in the present work are in very good agreement with the corresponding values of Ott et al. (OBHD) [56], albeit slightly lower (≈ 5%-10%). Good agreement is also observed when the nonequilibrium results for λ of Donkó and Hartmann [52] are used (DH-KT). The Green-Kubo data of Ott et al. (OBD-KT) [55] yield somewhat larger values (roughly 10%−20%, 25% larger than GH at = 70), as noted previously [56]. The longitudinal viscosity, ν l = 2 s − (γ − 1)D T , is shown in Fig. 5(b). The mutual agreement between data in the literature [53,54] and the present results is very good, with particularly small deviations in the range ≈ 30, . . . , 100. The GH and VE model yield almost identical results and even allow one to determine the famous viscosity minimum. The interpolation formula of Khrapak [67] provides an accurate analytical expression for the viscosity. Figure 3(d) shows that the inverse relaxation time, τ −1 = lim k→0 z 2 (k), decreases with the coupling strength, i.e., the transition from viscous to elastic response takes place at increasingly lower frequencies. Using the previous data for ν L , c s , and τ in the relation τ = ν l /(c 2 L − c 2 s ) (see Sec. III), one may compute a value for c L , which can be compared with the theoretical value. The latter is related to the high-frequency shear (G ∞ ) and bulk (K ∞ ) elastic moduli and the k → 0 limit of ω 2 L (k) = ω 4 / ω 2 [45,60],  [55] have been used to calculate D T via D T = λ/(n c p ), see the text for details. Also shown are the results for D T from Ref. [56] (OBHD). The shear viscosities of Donkó and Hartmann (DH) [53] and Daligault, Rasmussen and Baalrud (DRB) [54], and the interpolation formula of Khrapak [67] have been used to compute the longitudinal viscosity from ν l ≈ 4 η s /(3 n m). The bulk viscosity is negligible [68] and has been neglected.

D. Relaxation time and infinite-frequency sound speeds
The elastic moduli have been calculated as specific integrals over the pair distribution function involving derivatives of the Yukawa potential [45,69], Here, p is the pressure and g(r) the pair distribution function. The results are shown in Fig. 6. Figure 6(a) shows that the calculation of c L from the extrapolation of the various fit parameters is consistent with the k → 0 limit of ω L (k)/k, using the provided formula in Ref. [60]. The same applies to the data taken directly from the moments of the DSF and the calculation of c L via the pair distribution function. However, there is a shift between these data, mainly caused by an overestimation of the fourth moment in the VE model. 3 This could be related to the very fast decay of the DSF in the high-frequency limit (see also Ref. [42]), which becomes increasingly important for the high-order moments. The resulting deviations for c L are on the order of a few percent and decrease with the coupling strength, see Fig. 6(b). Even though the deviations from the theoretical value are rather 3 Similarly, the GH model overestimates the second moment.  small, they can significantly affect τ as c L and c s are of the same order [see Fig. 6(b)] and τ ∝ (c 2 L − c 2 s ) −1 , see the discussion below.
A related quantity is the infinite-frequency sound speed c ∞ = √ K ∞ /(mn) [70,71], which involves the (highfrequency) bulk modulus only. Note that in the present case, the kinetic terms ∼nk B T are included in the calculation of the elastic moduli. As shown in Fig. 6(b), c ∞ is in excellent agreement with the adiabatic sound speed c s , in line with previous results for Yukawa [70] and soft inverse power-law interactions [71]. Using the approximations c s ≈ c ∞ and ν l ≈ 4 η s /(3nm) (neglecting the bulk viscosity), and making use of the relation c 2 L = c 2 ∞ + 4 G ∞ /(3 nm), one obtains, for the theoretical value of the relaxation time, Here, τ M is the Maxwell relaxation time [45,72], which describes shear relaxation in a liquid. As shown in Fig. 7, the theoretical value for τ practically coincides with τ M . As could be anticipated from the discussion above, the VE value for τ is somewhat lower, in particular at intermediate coupling strengths. On the other hand, it is in rather good agreement with the inverse crossover frequency of the real and complex parts of the generalized shear viscosity determined by Goree et al. [72], which was used as an empirical measure of τ M . A more detailed investigation of these aspects would be required but is beyond the scope of this work.

E. Discussion
In summary, the analysis of the DSF at small frequencies and wave numbers can yield comprehensive insights into the hydrodynamic properties of Yukawa liquids. Compared to other approaches such as nonequilibrium or Green-Kubo methods, which typically only yield one particular transport coefficient, the method used here allows one to estimate several coefficients at the same time. On the other hand, it relies (i) on a specific model for the DSF and (ii) an extrapolation scheme, which must be sufficiently accurate to avoid systematic errors. Comparisons with a variety of available data sets show that the deviations from the results of other methods are relatively small, in particular for the sound speed ( 1%) and the viscosity (∼5%-10%). The largest deviations among different methods are found for the thermal diffusivity, which requires further investigation. For reference, the adiabatic sound speed and the transport coefficients are summarized in Table I.

V. CONCLUSION
It has been demonstrated that fundamental thermodynamic and transport properties of Yukawa liquids can be obtained with high accuracy from an analysis of the DSF. The VE model even allows one to quantify viscoelasticity in the strongly coupled liquid regime, giving access to properties beyond a purely hydrodynamic description. It remains to test the applicability of the method for a wider range of screening and coupling parameters. The analysis of Mithen et al. [35] suggests that it could be easier to apply the method at stronger screening, where the applicability range of the hydrodynamic model becomes larger. On the other hand, in the very weakly coupled limit, it has recently been shown [62] that the acoustic peak disappears below a critical , rendering the applicability of the hydrodynamic approach doubtful in this regime. Note that other extensions of hydrodynamics to finite ω and k (see, e.g., Refs. [36,60]) could be more appropriate than the GH or VE model under certain conditions.
For ab initio (density functional theory based) MD simulations in the WDM regime, it may become difficult to reach sufficiently small wave numbers since the number of particles that can be simulated is much smaller than for simulations with a pair potential. Thus, one may resort, e.g., to effective potentials deduced from average atom models [32,33]. The wave numbers required for an appropriate extrapolation, however, will depend on the applicability of the GH and VE model under WDM conditions. Further, treating electrons dynamically [31] has an effect on the DSF. Provided the very narrow acoustic and diffusive peaks in the ion-ion DSF at small wave numbers can be resolved in future experiments, which requires meV resolution, it could become possible to accurately determine transport quantities such as the viscosity from experimental data [24].

ACKNOWLEDGMENTS
The author would like to thank Zoltan Donkó, Torben Ott, and Luciano Silvestri for providing data from their works, Michael Bonitz for helpful discussions, and the referees for constructive feedback and suggestions. 033287-6