A viscous quantum hydrodynamics model based on dynamic density functional theory.

Dynamic density functional theory (DDFT) is emerging as a useful theoretical technique for modeling the dynamics of correlated systems. We extend DDFT to quantum systems for application to dense plasmas through a quantum hydrodynamics (QHD) approach. The DDFT-based QHD approach includes correlations in the the equation of state self-consistently, satisfies sum rules and includes irreversibility arising from collisions. While QHD can be used generally to model non-equilibrium, heterogeneous plasmas, we employ the DDFT-QHD framework to generate a model for the electronic dynamic structure factor, which offers an avenue for measuring hydrodynamic properties, such as transport coefficients via x-ray Thomson scattering.

Introduction

Access to high-power laser sources, such as the Linac Coherent Light Source (LCLS)[1], National Ignition Facility (NIF)[2] and Omega Laser[3], has opened the path to investigating essential properties of non-ideal plasmas such as ionization potential depression[4], transport coefficients[5] and ionization state[6]. Understanding the dynamical properties of non-ideal plasmas is critical for modeling and designing high energy-density science experiments, including inertial-confinement fusion[7], cluster explosions[8], laser-produced ion beams[9], hypervelocity impacts[10], in nanotechnology[11,12] and astrophysics[13].

Among all the approaches to modeling heterogeneous, non-equilibrium quantum systems, quantum hydrodynamics (QHD) is a computationally attractive approach with rich history in statistical mechanics. Shortly after the development of quantum mechanics, Bloch[14] proposed the first QHD model by simply choosing the Thomas-Fermi pressure for the electrons in an otherwise classical hydrodynamics model. In 1964, Hohenberg and Kohn[15] developed ground-state density functional theory for the inhomogeneous electron gas, which was immediately generalized to finite temperature by Mermin[16]. Combining the ideas of Bloch with DFT, Ying[17] proposed a new quantum hydrodynamic model via an adiabatic generalization of the density functionals. In Ying’s model, the pressure is represented by P[n(r,t)] with n(r,t) a time-dependent density described by the continuity equation. Ying’s QHD model includes explicitly all correlation and exchange effects included in the chosen energy functional. Using an alternate approach, Gasser and Jüngel[18] derived QHD equations using the Schrodinger equation with Wentzel-Kramers-Brillouin (WKB) wave functions. This approach yields the classical momentum equation with the Bohm potential but it does not account for correlations. Correlations effects and quantum degeneracy can be included in an ad hoc manner in this model by replacing the Bohmian potential with quantum potentials[12] or self-consistently through orbital-free density functional theory (OF-DFT)[19,20]. In yet another approach, using the moment expansion of the Wigner-Boltzmann equation, Gardner[21] proposed a QHD model for semiconductor devices that extends the classical hydrodynamic model to include quantum corrections. Similar results were obtained with the Wigner-Poisson system by Manfredi and Haas[22] for a quantum electron gas. Following Levermore[23], Degond and Ringhofer[24] used a non-commutative version of the entropy externalization principle to build a QHD model starting from the quantum Liouville equation. The moment equations are closed by a quantum Wigner distribution function that minimizes the entropy.

Despite these important advances, describing collisional processes in moderately coupled quantum plasmas remains a challenge[10,22,25]. Here, we explore an alternate approach based on a new formulation of quantum hydrodynamics (QHD). QHD approaches have the advantage of including equation-of-state and transport quantities more naturally than response-function approaches. Apart from these potential modeling advantages, QHD models of DSF therefore also provide access to experimental measurements of these quantities, thereby extending the utility of DSF. We develop a QHD framework based on the extension of the classical dynamical density functional theory (DDFT)[26-28], a variant of time-dependent density functional theory (TDDFT)[29]. DDFT provides a set of hydrodynamics equations by taking the velocity moments of Liouville equation and closes the system using density functional theory[26-28]. A fundamental assumption of this theory is that the equilibrium energy functional of the system can be used to guess the correlation energy functional when the system is out of equilibrium. While DDFT has found wide use in many-body classical systems[30-32], we extend its use in quantum systems[17] to viscous quantum systems, in general, and to DSF, specifically.

We apply the DDFT-QHD model to stationary, homogeneous and isothermal plasmas for which the dynamic structure factor (DSF) is well defined. While the DSF is of interest in its own right, it is also connected to x-ray Thomson scattering (XRTS) experiments; XRTS yields much essential information about plasmas, including density, temperature and atomic physics information (e.g., ionization state[6], ionization potential depression[4], etc.).

Results

Dense strongly coupled plasmas are characterized by large collisional contributions and degenerate electrons. These features make the DDFT-QHD approach a reliable tool for accurately describing the dynamical properties of these systems. For simplicity, here, we consider a quantum plasma comprising only electrons with density distribution n and mass m interacting through a pairwise Coulomb potential . We use atomic units (i.e., ) for the remainder of this work. The hydrodynamic equations for the electrons can be written generally as which are continuity and momentum equations written in terms of a generalized force tensor . Note that the continuity equation (1) and the left-hand side of (2) are generic, with the physical properties of the quantum electron gas entering through terms on the right-hand side of (2). In the DDFT approach[17,28,33], it is assumed that the system is close enough to equilibrium that an adiabatic closure can be chosen for ; that is, . The primary assumption of this model is that the system is near equilibrium, a condition well satisfied in highly collisional plasmas. Further, the equilibrium density is forced to be consistent with the thermodynamic ground state of the system by choosing the diagonal portion of the tensor to be of the form , where Ω is the free energy of the system. When is expressed using orbital-free density-functional theory (OFDFT), that portion of is closed. The off-diagonal portion of can be written in its long-wavelength form to yield a generalized Navier-Stokes equation of the formwhere η is the shear viscosity, and ξ is the bulk viscosity; all other symbols have their usual meanings. Provided η and ξ can be expressed in terms of , the hydrodynamic equations are closed.

In DDFT, one writes the total free-energy functional aswhere is the free energy of the noninteracting system, is the Hartree free-energy functional, and is the exchange-correlation (xc) functional. The Hartree term is exactly known and is an explicit function of space and time.

A key advantage to the DDFT approach to QHD is that all thermodynamic properties are included self-consistently through the total free energy Ω, for which a wide range of approximations are available[34-37]. In fact, this approach is very similar to the well-known generalized hydrodynamics, developed by Frenkel[38], that extends the classical Navier-Stokes equation to describe the properties of both solid and liquid bodies. Furthermore, our DDFT-QHD approach can be connected with other approaches based on Bohmian dynamics. If we set the viscous terms equal to zero in (3) and choose the gradient-corrected Thomas-Fermi (TF) functional for T[n], one recovers the well-known Bohmian QHD[20] form; again, however, the DDFT approach enforces self-consistency of its form with the other terms in the free energy. The connection between DFT and the Bohmian potential will be briefly shown below.

Density fluctuations are not readily available in density-functional theories, and our DDFT-QHD approach suffers from this limitation. However, in equilibrium, the fluctuation-dissipation theorem allows us to connect the linear response of the system to density fluctuations. We write the DSF of the electrons aswhere β is the inverse electron temperature and is the susceptibility of the free electrons. A large body of literature[39] focuses on the calculation of the system DSF , with most work based on the Chihara[40] decomposition

The quantities and are the Fourier components of the density of bound and free electrons. The first term of (6) corresponds to low-frequency electron-density fluctuations arising from ion dynamics and is proportional to the ion-ion DSF . The factor in the second term describes the contribution from core electrons[41] and is modulated by the ion self-motion . The third term is the free-electron DSF in the presence of a uniform ionic background. The quantity can be obtained from the standard Lindhard dielectric function within the random-phase approximation (RPA), or extended to include collisions as proposed by Mermin[42]. Thiele et al. generalized the Mermin form to include a dynamic collision frequency within the Born approximation[43], and Arkhipov et al. generalized the Mermin form to two-component plasmas, including sum rules[44].

The RPA results were also improved by including exchange and correlations through the local field corrections[45-52]. The ionic correlations contributions in a warm dense matter have been considered by Gregori and Gericke[53]. In this scheme, the strongly coupled effects of the ions are included through the different components of the memory function constrained by the sum rules[54]. This phenomenological approach has been applied successfully in Coulomb liquid[54-57] community for systems where the memory functions have a Gaussian or exponential form. However, for more complex systems, the form of the memory becomes mathematically intractable. Schmidt and coworkers[58] have proposed a hydrodynamic model that begins with moments of the Wigner-Poisson system with a collision term added. In such an approach you cannot describe correlations properly since the resulting pressure term is of an ideal gas. The DDFT-QHD approach we introduce here accounts for self-consistently many-body physics effects and also non-local hydrodynamic effects through the choice of the free-energy functional.

The linear susceptibility associated with a weak external potential that induces a disturbance in the electronic density is defined asThus, the susceptibility can be determined by linearizing the quantum hydrodynamics equation and using (7). To do, we first expand the density and velocities about a uniform mean as which yields the linearized QHD equations in Fourier space: where , and the tilde sign denotes the Fourier transform. By combining (10) and (11) and using (7), we obtain an expression for the electron susceptibility:where is the longitudinal viscosity. To proceed, we need to choose specific forms for the different contributions of the free-energy functional Ω[n]. The free-energy functional is typically chosen to ensure that an accurate equilibrium density is recovered, although exact analytical forms are generally not known. However, the contributions of the excess free-energy functional to the free energy of the system, , can be expressed formally in terms of the direct correlation function as follows[59]:where , and μ ex is the excess chemical potential. Once the pair potential has been specified, the self-consistent contributions of the excess free-energy functional can be calculated using the direct-correlation function via path integral quantum Monte Carlo (PIMC) simulations[60,61], integral equations[62] or analytical fits[63]. Let us now evaluate the free-energy functional of the non-interacting electron gas] T[n].

Many approximations for T[n][64-66], have been described in the literature- Thomas-Fermi (TF), Kirzhnits gradient correction (TFK)[65], von Weizsäcker funtional (vW)[64], Perrot functional[66] to name a few. Most of these models are based on an extension of the TF functionalwhereis the Fermi-Dirac integral of order p, is the chemical potential normalized with the k T and the electron density is given byWe consider here the functional form for an electron gas based on the TF functional with the finite-temperature Kirzhnits gradient correction (TFK)[65]:Here, we introduced a coefficient γ that allows to capture a variety of results. First, the systematic gradient expansion of Kirzhnits yields the prefactor γ = 1/9. Second, the von Weizsäcker result follows by a partial integration and γ = 1. The assumption in this gradient-correction expansion is that the error made by neglecting the third- and higher-order terms is very small. For high-density plasmas, interface-mixing problems or shock structures in which temperature and density gradients can be large, this expansion ceases to be valid. In such circumstances, it may be important to include higher-order terms for the thermal terms through higher gradient corrections q in the TFK functional. The non-interacting free-energy functional T[n] can also be expressed in terms of the Lindhard function, which is exactly known, instead of using OFDFT[66,67]. Let us now show the connection between DFT and the Bohm description. The functional derivative of the kinetic energy functional (17) is given by[20] where , its derivative with respect to the density is denoted . The first term of (18) is the Fermi pressure while the second term corresponds to the generalized Bohm potential in the finite-temperature regime, revealing the underlying connection between DFT and the Bohm description. The connection between DFT and the Bohmian picture has recently been discussed in a paper by Stanton and Murillo[20] where a Kirzhnits[65] correction was used to get the “Bohm” term quite generally.

After linearizing and taking the Fourier transform of (13) and (17), the contributions of the non-interacting and interacting electron gases to the susceptibility becomewhere the TF length and the parameter are given byrespectively. Substituting (19) into (12), we findHere the frequency ω is in units of the electron plasma frequency, is the wave number, is the Wigner-Seitz radius, is the Thomas-Fermi length in the units of the Wigner-Seitz radius, the viscosity is in units of , and the coupling parameter, defined as the ratio of the potential energy to the average kinetic energy, is given by .

Substituting (21) into (5) yields the free-electron DSF

Equation (22) is the main result of this work. The second factor is the usual Bose function. The denominator of the third factor includes quantum degenerate plasma effects through the direct correlation function , thermal effects with high-order gradient terms, and viscous damping through . It is worth noting that when the degeneracy parameter is very large, electrons can be considered to be in a non-degenerate, classical state. If we then replace the exponential in the Bose function by its Taylor expansion, we recover the dynamic structure of non-degenerate electrons given by the Navier-Stokes model[68].

The form of the DSF obtained from this theory, without the dissipative effects, is connected to the approaches based on the local field corrections[49]. The direct correlation function which is the main ingredient of this approach is related to the local field correction (LFC) as:

In the random phase approximation (RPA), the direct correlation function is given by consequently the LFC vanishes, . Thus, the direct correlation function describes the strongly Coulomb correlation and exchange effects beyond the RPA. Several approximations for the LFC[45,46,49,69], has been proposed starting from the formulation of local field corrections due to Coulomb correlations and exchange effects by Hubbard[45]. Utsumi and Ichimaru[49] formula has been widely used to investigated the static properties of systems at metallic densities. Holas, Aravind and Singwi[70] have suggested an expression for the dynamical LFC in strongly coupled electron gas. Although, we can use existing analytical fits for the LFC to obtain the direct correlation function, we choose here to computed this quantity directly using Ornstein-Zernike equations with the hypernetted-chain approximation closure[71]. Furthermore, matter under extreme conditions of temperature and pressure undergoes large spatial gradients (i.e., shocks structure, interface problems, etc.). The heterogeneity can greatly altered the Thompson spectra with respect to the uniform case as recently discussed by Kozlowski and coworkers[72]. In our approach, the constitutive equations will relax to the correct DFT thermodynamic ground state, which other methods cannot guarantee. This means we have the full non-local correlations absent from most other approaches opening up the possibility of studying the cases for which the usual homogeneous and isotropic forms like are not applicable. This can be done through the introduction of the inhomogeneous direct correlation function . In the next section we will focus on the characteristic features of our main result (22).

Discussion

Computing the DSF (22) requires knowledge of the viscosity and the direct correlation function. In the literature, the DSF is often expressed in terms of the local field correlation, and the latter quantity is evaluated using analytical fits[52,73,74]. The direct correlation function can also be obtained directly through numerical simulations; that is the avenue pursued in this work, using hypernetted-chain calculations[62,75], with a quantum statistical potential (QSP)[76-78]. We use QSP approach here to merely have easy access to results for which we can illustrate the DDFT method, which is the main point of the paper; other methods can be also used to get the structure information for the DDFT model. In fact, we see a strength of QSPs in this regard: the key quantity is the electron-electron , which is not accessible from DFT approaches. Electron-electron correlation functions are available from PIMC, however, and that provides validation for our input quantities. Jones and Murillo[79] have shown the theoretical underpinnings of QSPs and reviewed their extension to fully degenerate quantum systems. Dutta and Dufty[60] have compared compared QSP-based RDFs from the modified Kelbg QSP and PIMC and show that over an extremely wide range of physical conditions the QSP predictions are nearly perfect; it is only at very low densities that we can see a modest deviation. Here, we choose the QSP from the pioneering work of Hansen and McDonald[76] since they yield results similar to the more complicated Kelbg potentials. Comparisons between Coulomb, HM and Kelbg potentials are shown in the online Supplementary Material. Next, we turn to the electronic viscosity is needed, which is determined by both electron-electron and electron-ion collisions. The electron viscosity is obtained by interpolating the zero-temperature viscosity proposed by Conti and Vignale[80] and the finite-temperature viscosity for classical plasmas by Stanton and Murillo[33]. When building this fit, we considered contributions only from electron-electron interactions. However, the electron viscosity should also take into account electron-ion[81-83], contributions, which can be more significant than the electron-electron viscosity in the regimes of interest. Please see the online Supplementary Material for a detailed description of the calculation of these two quantities.

Figure 1 shows the spectra of the DSF for different values of the wave number q, the coupling strength Γ and the density parameter r . The free electron DSF is normalized by its maximum value. In Fig. 1a, the positions of the plasmon peaks (Stokes and anti-Stokes) and its amplitude remain almost unchanged when the quantum parameter r increases from 1.0 to 4.0. The reason is, in this regime , the correlations and quantum degeneracy effects are negligible and consequently have no impact on the propagation of the plasmon. Furthermore, the width of the plasmon peak reduces when q and r increase owing to the fact that the viscosity which acts to broaden the width of the peak is very sensitive to the density parameter r . Figure 1a–c show that the position and width of the plasma peak vary strongly with Γ. These figures also display some of the standard features of the plasmons peaks; they are symmetric with respect to the zero frequency, and the difference between their amplitudes gives a measure of the electron temperature through the detailed balance relation. It is worth noting that for large values of the coupling parameter and density, the plasmon peak is at a frequency smaller than the plasma frequency .

Figure 1

Effects of the coupling parameter on the spectra. We show the variation in the DSF for different values of the coupling parameter Γ and the normalized wavelengths . The coupling parameter Γ ranges from 0.2 to 0.8, and r ranges from 1.0 to 4.0. The dynamic structure factor is normalized by its maximum value. The two plasmon peaks are symmetric and the ratio of their amplitudes gives a measure of the electrons temperature.

For a given value of the wave number q, the peak of corresponds to the dispersion relation of the plasmon . According to (22), the dispersion relation is approximately given by

In the RPA limit, with a pure Coulomb potential, the direct correlation function is given by . By substituting this result into (24) and setting the viscosity equal to zero, we recover the Bohm-Gross dispersion relation[84] for an isothermal plasma. The dispersion relation of the plasmon is shown in Fig. 2 for and for Γ = 1.0 and Γ = 0.7. The red triangles lines show the DDFT-QHD result (24), and the data points indicated with blue line corresponds to (24), with the viscosity set to zero, η  = 0. From our basic result (24) we can explore several limits that yield other models. For example, because the direct correlation function is related to the local-field correction through the relation , we can neglect the three higher-order terms in (24) to obtain the local-field correction (LFC) result. We show this limit in Fig. 2 as a green line. Next, we can retain the second and third terms to example quantum corrections to the LFC result, and that model is shown as blue triangles. The full result, including viscosity is shown as the red line; note that the viscous correction is large, suggesting that the power series in q of (24) is probably not converged at the largest values of q in the plot.

Figure 2

Plasmon dispersion relation. We show the frequency as a function of the wave number q for and (a) Γ = 1.0 and (b) Γ = 0.7. DDFT-QHD refers to (24), which accounts for strong correlations and viscosity. The label “DDFT-QHD: ” corresponds to (24) with the viscosity set to zero. The green curve shows the local field correction dispersion relation.

The dispersion relation (24) implies that the width of the peak of the DSF is broadened by the viscosity. Therefore, by measuring the width of the peak of a spectrum, information about the electron viscosity can be inferred. We can obtain a good estimate of the width in the following way. We know the location of the peak from the dispersion relation (24). Near that peak, we know what ω is, and this information can be put into to obtain the width of the DSF. This scaling gives the width in terms of both η and , which means that η cannot be determined without knowledge of the direct correlation function. However, by fitting the entire spectrum of the DSF, all quantities can be obtained: the density, the temperature, and the viscosity. The DSF spectra suggest that measurements at a few q values is best.

Concluding Remarks

A general framework for electron dynamics is provided with DDFT-QHD. Our model is the specific approximation of TDDFT in which the generalized-force functional is replaced by the equilibrium functional. We established the connection between DDFT-QHD and DSF through the fluctuation-dissipation theorem, allowing for improved QHD models to be compared with experimental data. The predicted DSF spectrum exhibits strong correlations and collisions that are built self-consistently into the model; this result differs from those obtained with more common Lindhard approaches[42-44], in which collisions enter through a dynamic collision frequency[43,44], or though local field corrections[85]. Our result suggests that the electronic viscosity can be determined experimentally by measuring the electron DSF.

Our approach is a full hydrodynamics model that can be used to simulate non-equilibrium, heterogeneous dense plasmas[72]. For example, we could investigate shock physics, fluid instabilities, and large-scale experiments. Most other methods are based on stationary, homogeneous/isotropic approximations; this is explicit in functions such as . Thus, while our DDFT formulation of QHD is significantly beyond these simple linear response functions, we show here that we are able to make contact with the XRTS community and connect the scattering spectrum to transport coefficients in a direct way with a hydrodynamic approach.

Finally, our model still lacks an equation for the energy fluctuations. Past experience suggests[54-57,86], that energy fluctuations might cause a zero frequency mode. This latter was experimentally measured in a liquid lithium by Sinn and coworkers[87] confirming molecular dynamics simulations results performed by Canales, González, and Padró[86]. Our DDFT-QHD approach would miss any mode originating from thermal fluctuations[54,58,88], because it is based on an isothermal assumption. We think this approach can incorporate an energy equation, but this is work in progress. In future work, it would be useful to include other transport quantities, such as the viscoelastic relaxation time and the thermal conductivity. Extension of this model to the full XRTS form factor with an electron-ion generalization of the DDFT-QHD equations[89,90], is left for the future.