Theoretical foundations of quantum hydrodynamics for plasmas

Beginning from the semiclassical Hamiltonian, the Fermi pressure and Bohm potential for the quantum hydrodynamics application (QHD) at finite temperature are consistently derived in the framework of the local density approximation with the first order density gradient correction. Previously known results are revised and improved with a clear description of the underlying approximations. A fully non-local Bohm potential, which goes beyond of all previous results and is linked to the electron polarization function in the random phase approximation, for the QHD model is presented. The dynamic QHD exchange correlation potential is introduced in the framework of local field corrections, and considered for the case of the relaxation time approximation. Finally, the range of applicability of the QHD is discussed.


I. INTRODUCTION
The investigation of dynamical properties of systems containing partially or fully degenerate electrons has gained growing interest due to their relevance for such fields like dense plasmas [1,2], warm dense matter [3,4], streaming and wake effects [5,6] and plasmonics [7][8][9][10]. Interiors of giant planets, white and brown dwarfs, stellar cores, and the outer envelope of neutron stars are examples of matter in the state of a partially degenerate dense plasma (warm dense matter) [11,12]. Experimental studies of dense plasmas include the free electron laser excited plasmas [13], and inertial confinement fusion experiments [14][15][16]. On the other hand, plasmonic materials, containing fully degenerate electrons, have recently received renewed attention as the result of the advances in nanofabrication techniques [8,9,17,18]. All the above mentioned different systems are governed by the behavior of Coulomb interacting quantum electrons.
The theoretical description of quantum plasmas must take into account quantum degeneracy effects such as non-locality, spin statistics, and correlations (nonideality) appropriately on the relevant scales. A continuum description of the dynamics of the quantum electrons in the spirit of a density functional theory (DFT) is a promising approach to the problem. The possibility of such a description follows from the Hohenberg-Kohn theorem of DFT [21,22] and the Runge-Gross theorem of time-dependent (or current) density functional theory (TDDFT) [23]. The dynamics of the electrons can be described in terms of the average electron density n(r, t) and the average electron current density j(r, t) by the continuity and momentum equations [19,23]: m ∂ ∂t j (r, t) − n (r, t) eE(r, t) = −∇ · P(r, t), where E is the electric-field strength, and P stands for a tensor that contains all many-body and quantum effects, including correlations and dissipation. The exact form of P is not known, and different approaches and approximations for P have been proposed with different levels of sophistication, e.g. [24,25]. An alternative concept are first-principle approaches based on wave function methods, e.g. [26], quantumstatistical theory [1], quantum kinetic theory [20,27] or non-equilibrium Green functions [28,29]. However, as TDDFT, these methods require substantial theoretical and computational effort and are particularly important to capture electronic correlations. Therefore, in cases where correlations and their dynamics are of minor importance, simpler approaches are being used. This particularly applies to quantum hydrodynamics (QHD) that became popular as a simplified approach for quantum plasmas [30][31][32][33][34], plasmonics [8][9][10]35], and electrons in semiconductors [36,37] and is, therefore, in the focus of this paper.
The key ingredients of QHD -as it is often used in the context of quantum plasmas -are the ideal Fermi pressure, P F , and the so-called Bohm potential, V B [38]. In QHD, the closed momentum equation is, instead of Eq. (2), Manfredi and Haas [31] derived the Fermi pressure and a quantum correction in the form of the Bohm potential using a semi-classical Hartree ansatz for the N -electron wave functions with identical amplitude for all singleelectron orbitals [38]. However, in order to reach agreement with the results of the more fundamental kinetic theory in its simplest form -the random phase approximation (RPA) -both, the Fermi pressure and the Bohm potential have to be "corrected" by constant pre-factors [35,[43][44][45]. Similarly, in plasmonics, the QHD theory is used with one or, sometimes, two fitting parameters corresponding to the prefactors of the Fermi pressure and the Bohm potential, but with an additional exchange correlation potential contribution, which is valid only in the static case.
Moreover, in the context of quantum plasmas, QHD has often been used beyond the range of applicability of the model and, occasionally, even with explicitly incorrect expressions. This has led to un-physical predictions and some controversy, for a discussion, see Refs. [38][39][40][41]48].
However, introducing the above mentioned corrections factors does not solve the problem. It turned out that these fitting parameters (pre-factors) are not constants but differ depending on the characteristic wavelength and frequency of the physical problem. In addition, these coefficients are found to vary with the plasma density and temperature. This results in complicated parametric dependencies of the Bohm potential. This means that the range of the values in which correcting parameters can be chosen and the corresponding underlying physical assumptions need to be clarified.
For this reason, in this paper QHD is reconsidered in detail, for both zero temperature and finite temperature, and it is put into the context of well established approaches such as Thomas-Fermi theory and dielectric theory (such as the random phase approximation). Starting from these approaches it becomes more clear what approximations are actually being made and what is the corresponding range of validity of the model.
In Sec. II, continuity and momentum equations of the QHD model are briefly introduced for the finite temperature case. In Sec. III, the relation of the Bohm potential in the density gradient approximation to the power expansion of the inverse finite temperature RPA polarization function is presented, and the QHD potentials, i.e., the potential related to the Fermi pressure and the Bohm potential, are considered in different limiting cases. This will allow us to reproduce known results and, in part, even to improve them. In Sec. IV, the generalized nonlocal Bohm potential, based on the exact RPA polarization function, is presented. The exchange-correlation potential for the QHD application is discussed in Sec. V. The paper is concluded by a discussion of the range of applicability of the QHD model.

II. QHD EQUATIONS AT FINITE TEMPERATURES
The underlying equations of the QHD model can be derived from a field theory, starting from the semi-classical Hamiltonian which, in the absence of a magnetic field, reads [35,46] H[n(r, t), w(r, t)] = E[n(r, t)] − eV ext n(r, t)dr where w is the scalar potential determining the velocity is the sum of the kinetic and the exchange-correlation energy functionals, and V ext refers to the external electric potential. Using n(r, t) and m e w(r, t) as canonically conjugate field variables, we employ Hamilton's equations [49]: and obtain the following equations of motion that form the basis of QHD [46,47]: Introducing the potential of the generalized force [48] µ[n(r, t)] = δE[n(r, t)] δn(r, t) + eϕ(r, t), with the definition ϕ(r, t) = e n(r , t) and making use of relations v = −∇w and (v · ∇)v = 1 2 ∇(∇w) 2 , we arrive at the QHD equation in terms of average density n (r, t), velocity v (r, t), and the generalized In previous works, Eq. (12) was obtained for the case of fully degenerate electrons (zero temperature limit). Equations (11) and (12) represent QHD equations in the micro-canonical ensemble, as they were derived from the semi-classical Hamiltonian (4) (alternatively, they can be obtained from a semi-classical Lagrangian [46]).
For the extension to finite temperature, we now switch to the grand canonical ensemble. There, these equations have the same form, but n and w must be understood as quantities that are averaged over the grand ensemble [50,51]. Now we generalize the momentum equation (12) to the finite temperature case where it is advantageous to use the free energy functional, F [n], instead of E[n]. Indeed, in the grand canonical ensemble we have [51] with the grand potential where µ 0 is a constant defining the chemical potential of the system in thermodynamic equilibrium, and N = n(r)dr corresponds to the mean value of the number of particles in the grand canonical ensemble. The derivation of Eq. (13) is given in the Appendix A.
It is should be noted that, in Eq. (13), E is the average value of the energy over a grand canonical ensemble. With this, we obtain for the potential of the generalized force at finite temperature: which provides the link with the standard fluid theory, cf. Eq. (2).

A. General Expressions
In order to derive a potential related to the Fermi pressure and the Bohm potential, we neglect the exchangecorrelation term, F xc [n], and turn to the local density approximation (LDA) with non-locality taken into account by the first order gradient correction to the noninteracting free energy functional [45,53,54]: where the free energy F 0 is defined via the free energy density Now we show that it is possible to connect the QHD Bohm potential with the expansion of the inverse polarization function and, thereby, to obtain a 2 . For the static case, ω = 0, this was done by Perrot [55]. Here we generalize his results to the dynamic case. We consider small localized variations of the potential, δϕ, density, δn, velocity in terms of the scalar potential, δw, and of the chemical potential, δµ id . From Eq. (7), taking the equilibrium density distribution as uniform (n 0 = const) and assuming w 0 = 0, we obtain where δñ, and δw denote the Fourier transforms of δn and δw, respectively. The resulting variation of µ id , as defined by Eq. (19), has the form [55] Thus, the linearized momentum equation (12) forw = δw and, taking into account Eq. (21), is written as: − iωm e δw = eδφ+ Finally, making use of Eq. (20), we find This is an important result that relates the density perturbation to the perturbation of the external potential. As we have assumed weak perturbations, we can make use of the results of linear response theory and identify in Eq. (23) the inverse polarization function, Π −1 ≡ eδφ/δñ: B. RPA result for the Bohm potential and Fermi pressure Equation (24) gives us the opportunity to express the r.h.s. of Eq. (24) systematically via linear response theory. The lowest order many-body approximation for Π is the random phase approximation (RPA), which reads in equilibrium, for arbitrary temperature [56] where we introduced dimensionless frequency and wave number, u = ω/(kv F ), z = k/(2k F ), and defined χ 2 0 = (πk F a B ) −1 r s /6.03, where a B is the Bohr radius, r s ≡ a/a B , with a being the mean interparticle distance, and we also introduced the Fermi wave number and the plasma frequency, k F = (3π 2 n) 1/3 , ω 2 p = 4πne 2 /m e . Finally, we defined With these explicit expressions the RPA polarization can be studied in detail. Particularly simple expressions exist for various limiting cases with respect to frequency and wave number. In the limiting cases of large or small values of z, the inverse of the real part of the RPA polarization function has the following expansion [56,59]: In particular, we obtain the long-wavelength limit of the inverse polarization function as Comparing Eqs. (24) and (27), we see that The convergence of the expansion (27) is discussed in Appendix B for the static case and k < 2k F . We stress that the coefficientsã 0 and a 2 depend on the considered limits for k and ω which directly affects the value of the Bohm potential, as will be shown below. Equations (29) and (30) are closely related to the stiffness theorem connecting the perturbation of an equilibrium state (of the energy, in the case of the ground state) due to an applied external field with an inverse linear response function. The proof of this theorem was given in Ref. [19] for the ground state and can be easily extended to finite temperatures.
Once the coefficientsã 0 [n 0 ] and a 2 [n 0 ] are defined, we allow the equilibrium density to vary in space and time, n 0 → n(r, t), according to the standard concept of LDA. The Bohm potential in the local density approximation, for QHD applications can be obtained for any degeneracy parameter from the second term on the right hand side of Eq. (17): The coefficientsã 0 [n] and a 2 [n] are expressed in terms of the density n(r) and the dimensionless chemical potential, η = βµ, where we introduced the inverse tempera- Using the local density approximation, we obtain The following partial derivatives are needed to find the Bohm potential (31), where , and the Bohm potential has the following form: where the coefficient γ sensitively depends on the considered values of the wave number and frequency.
On the other hand, the Fermi pressure is proportional to the functional derivative of the Thomas-Fermi free energy [48], which, at θ → 0, can be written as: where the coefficientᾱ depends on the considered limit on the k-ω plane. Previously, this coefficient was chosen by adjusting the QHD result for the longitudinal plasmon dispersion to the RPA prediction [43,61] for zero temperature [note that F → E, at θ → 0]. Now, we will use Eq. (29), to analyze different limits for the coefficientsã 0 [n] and a 2 [n] = −ã 2 , for different frequency-wavenumber ranges.
C. Low frequency and long wavelength limit, The results given for the static case remain valid also for low frequencies, i.e. u 1 or ω kv F . This regime is relevant for a variety of physical processes such as the screening of a test charge and the dispersion of lowfrequency waves (e.g. ion-acoustic waves) in a plasma. In a variety of publications e.g. the topic of screening of an ion in a quantum plasma were treated incorrectly using the Bohm potential for the high-frequency case rather than for the static case. This has lead to the unphysical prediction of ion-ion attraction in an equilibrium quantum plasma, for a discussion, see Refs. [40,45]. To obtain the Fermi pressure and Bohm potential in the present limiting case, we calculate the coefficientsã 0 [n] and a 2 [n] and obtain [45]:ã with In the low-temperature limit, θ → 0, we have H 1 (η) 1 and H 2 −1, and Eq. (37) gives the asymptotic results a 2 → 2 /(72mn) and γ = 1/9 for the coefficient in front of the Bohm potential in Eq. (34) [45,55,57,58].
Using the relations (32) and (33), the finitetemperature generalization of the Bohm potential can be easily found by substituting Eq. (37) into Eq. (31) and taking into account the dependence θ[n] ∼ n −2/3 . Furthermore, for the static case (ω = 0), we have f 0 ≡ f TF , where f TF is the Thomas-Fermi free energy density: where I ν is the Fermi integral of order ν. Using Eq. (38) Perrot [55] showed that the second order partial derivative of the Thomas-Fermi density of states, − 1 2 , in the static long wave length limit (k 2k F ) exactly coincides withã 0 [n] given by Eq. (36). Finally, the func-tional derivative of the Thomas-Fermi term yields: At zero temperature, η = βE F , leading to δT TF /δn = E F , in agreement with the result of Ref. [48]. By regrouping terms we can rewrite the Bohm potential (31) in the following form: where and and the last line of Eq. (42) was obtained using relations (32) and (33). In order to analyze finite temperature effects we introduce the coefficient: similar to the zero-temperature case, and use it in V 1 to obtain By taking into account that C = 4 3n θ −3/2 , we find: Finally, we introduce the coefficient: and substitute Eqs. (44) and (45) into Eq. (40). This yields the following expression for the finite-temperature Bohm potential in the static long wavelength limit: The correction coefficientγ of the first term in the Bohm potential for the finite-temperature case (47), and the dependence of γ and γγ on θ are presented in Fig. 1. In the limit of fully degenerate electrons, θ → 0, as well as in the classical limit, θ 1, the correction coefficient approaches unity,γ → 1. This correction coefficientγ is important for a partially degenerate plasma, around θ ∼ 0.5, as can be seen from Fig. 1.
Recently, Haas and Mahmood [60] considered the finite-temperature Bohm potential by analyzing ionacoustic waves on the basis of linearized QHD equations and comparing them with the RPA result. They correctly derived the coefficient γ in Eq. (43), but missed the correction coefficientγ in Eq. (46). Note that in Eq. (47), the second term of the Bohm potential is proportional to n 1 /n 0 , whereas the first term ∼ (n 0 /n 1 ) 2 , where n 1 is a small density perturbation (i.e., n 0 /n 1 1). As a result, in the linear approximation, which was considered by Haas and Mahmood [60], the information about the first term of the Bohm potential and the coefficientγ is lost.
For the degenerate electron gas, Jones and Yang [61] showed that, in the short wavelength limit, the first order gradient correction term of the non-interacting kinetic energy has the form of the von Weizsäcker gradient correction with a 2 [n] = 2 8men , which gives γ = 1 for the Bohm potential. This result is important as the von Weizsäcker gradient correction correctly reproduces Kato's cusp condition for the electron distribution close to a test charge (core) [63]. Knowledge of the change of the coefficients determining the Fermi pressure and Bohm potential at the transition from the long wavelength limit to the short wavelength case can be useful to analyze the QHD results in plasmonics [10,65], quantum plasmas [64] and the application of orbital-free DFT to dense plasmas (warm dense matter) [2,66,72].
In the low-frequency short wavelength limit, u z, z 1, we use the following expansion of the function ∆g = g(u + z) − g(u − z) [56]: where From Eqs. (48) and (49) we obtain which can be substituted into the formula for the inverse polarization function, with the result From this, we obtain the coefficientsã 0 [n] and a 2 [n] in the short wavelength case: From Eq. (53) we see that the coefficient in front of the Bohm potential in Eq. (34) becomes γ = 1 [61], and it is interesting to note that, in the short wavelength limit, the Bohm potential is independent of temperature. Consider now the Fermi pressure term. At small temperatures, θ 1, using Eq. (52) for the functional derivative of F 0 [n] yields where the first term is the result for the ground state and the second term is the finite-temperature correction. Thereby, the transition from the long wavelength to the short wavelength limit leads to a change of the factor α, from 1 to 3/5. Note that, in the short wavelength limit, or for large values of the density gradient the von Weizsäcker gradient correction is the leading term in the non-interacting free energy functional of electrons [63]. In Fig. 2 the results of the expansions of the inverse RPA polarization function in the limits of long, k/2k F 1, and short wavelengths, k/2k F 1, are compared with the exact RPA result. It can be concluded that, in the case of the uniform electron gas, the long-wavelength limit result is applicable up to k k F . E. High frequency limit, ω k 2 /2m The present limiting case is important for a variety of physical situations, such as for the description of the optical (Lengmuir) plasmon of the electrons as well as for the plasma response to a high-frequency electromagnetic field with frequencies exceeding the plasma frequency.
To obtain the correct pre-factors of the Bohm potential and Fermi pressure in the present high-frequency limit, u z, we use the formulas [56] g and obtain z ∆g The coefficientsã 0 [n] and a 2 [n] are obtained from Eqs. (57) and (51) Equation (59) shows that the coefficient in front of the Bohm potential equals γ = 1. In the high-frequency limit the coefficient a 2 [n] does not depend on temperature, which means that at finite temperature the Bohm potential is given by Eq. (34). This is explained by the fact that, at sufficiently high frequency of the perturbing electric field, the back and forth movement of the electrons is not affected by their thermal motion.
From Eq. (58), the high-frequency analogue of the Thomas-Fermi pressure can be obtained using − 2ã 0 [n]dn. In the zero-temperature limit, taking into account the relation I 3/2 (θ → 0) = 2 5 θ −5/2 leads to the following expression for the functional derivative of F 0 [n], at high frequency: Equation (60) indicates that, in the high-frequency limit, the QHD equations contain the coefficientᾱ = 9/5. We note that, in previous works [35,43], this constant coefficient was artificially added in order to reach agreement between the results of the QHD and the RPA expression for the plasmon dispersion relation. A finite-temperature correction to Eq. (60) can be obtained by expanding I 3/2 (η) around θ = 0, The different values of the factors γ andᾱ in the k-ω plane are summarized in Fig. 3. Furthermore, the information about the obtained coefficients a 2 [n], γ and a 0 [n],ᾱ is listed in tables I and II.  Tables II and I).
Using the high-frequency result forã 0 and the Bohm potential, one can derive the well-known plasmon dispersion from the continuity and momentum equations: which, at θ 1, approaches the form where we took into account that 2ã 0 [n 0 ]n 0 /m e → − 3 5 v 2 F , at θ → 0, and that, in the high-frequency limit, a 2 [n] = 2 /8m e n.
On the other hand, in the limit θ 1, taking into account thatã 0 [n 0 ] → − 3 2 k B T n0 , leads to the dispersion where v th = k B T me is the thermal velocity of the electrons.
Note that the linearized QHD equations correctly reproduce the dynamic RPA polarization function in the long wavelength limit, Π 0 (ω), Eq. (28), without any additional terms related to the fermionic pressure.

IV. GENERALIZED NON-LOCAL BOHM POTENTIAL IN LINEAR RESPONSE
In the previous section, the Bohm potential of quantum hydrodynamics was derived on the basis of the RPA polarization function, using the expansion of the latter in powers of the wavenumber. This has allowed us to improve the QHD model, depending on the wavenumber and frequency range in three relevant cases, by involving the dynamic LDA and the gradient correction, on the basis of the RPA On the other hand, QHD can be used to independently obtain the response function of (uncorrelated) electrons, Π id QHD (k, ω), for arbitrary frequencies and wavenumbers. In particular, at large wavenumbers the previous power expansion is not applicable. It is, therefore, instructive to inquire whether this expansion can be entirely avoided. The idea is again the enforce agreement with the polarization function of linear response theory, but now in the whole frequency-wavenumber range. The simplest solution is, again, to use the full (non-local) RPA polarization and to require This will result in a generalized non-local Bohm potential. We now derive a relation between the free energy of the system and the RPA polarization. To this end, we consider an equilibrium density profile n 0 (r) that is current-free, w 0 (r) = 0 and which follows, as before, from δF [n 0 (r)] δn 0 (r) where µ 0 is a constant. The equilibrium configuration is now exposed to a weak external perturbation giving rise to n(r, t) = n 0 (r) + n 1 (r, t), w(r, t) = w 1 (r, t), ϕ(r, t) = ϕ 0 (r) + ϕ 1 (r, t). The first order perturbations follow by taking the continuity equation in first order in the perturbations, and, similarly, the momentum equation [46] m e ∂w 1 (r, t) ∂t where w 1 determines the velocity perturbation via v 1 = −∇w 1 , and the last term on the right hand side of Eq. (68) represents the total potential of the fermionic pressure which includes the Fermi pressure, Bohm potential and, in general, the exchange-correlation terms. Now we again assume the equilibrium density to be uniform and consider small fluctuating quantities n 1 , w 1 . From Eqs. (67) and (68) we obtain, after Fourier transform to frequency and wavenuber (k, ω) space, which we denote by F −iωñ 1 + k 2 n 0w1 = 0, From Eqs. (69) and (70), immediately follows the QHD result for the inverse polarization function Π −1 QHD (k, ω) = eφ 1 /ñ 1 : Neglecting the exchange-correlation contribution to the free energy and taking Π QHD = Π RPA , it is straightforward to deduce from Eq. (71) that: where we used the definition (28) of the long-wavelength limit, Π 0 (ω), of the RPA polarization. As it was mentioned previously, Eq. (72) incorporates both the Fermi and Bohm potentials. In QHD, Eq. (72) can be used without additional separation of different contributions. Now we verify that the generalized result, Eq. (72), correctly reproduces the results for the Bohm potential for the special cases that were obtained in the previous section. We do not use the (ideal) free energy density in the form of a gradient expansion, Eq. (17), but, instead, start from the more general non-local form [21,45]: (73) where the kernel K is symmetric with respect to permutation of r and r . Using Eqs. (18) and (29), Eq. (73) can be written as: where K ([n]; k) is the Fourier transform of the kernel K. Substituting Eq. (74) into Eq. (72) we find for K, after the inverse Fourier transformation, and, for the generalized non-local Bohm potential, we have where we have used the abbreviation n 1 (r) = n 1 (r, t).
Utilizing the expansion (27) we arrive at and obtain, after inverse Fourier transformation Considering only the first term on the right-hand-side of Eq. (78), one can find the Bohm potential from (76) in the form of Eq. (31). Furthermore, higher order terms give rise to higher order gradient corrections to the noninteracting free energy density functional [45,[67][68][69][70]. For more details on the convergence of the gradient expansion in the ground state, we refer the reader to Ref. [71] and the references therein. The closure relation Eq. (72) for the QHD equations (11) and (12) provides a unified general picture for the understanding of the complex parametric dependencies of the pre-factors of both Fermi pressure and Bohm potential on frequency, wavenumber, density and temperature. Moreover, it indicates ways how to systematically go beyond both the local density approximation and the model of an ideal Fermi gas. The inclusion of exchange and correlation effects will be performed in Sec. V whereas further issues of non-locality will be discussed in Sec. VI.
On the other hand, the local version of QHD, i.e. LDA with first order density gradient corrections, has now also been clarified. Indeed, there is no inconsistency in using one set of pre-factors in front of the Bohm potential and Fermi pressure for computing the equilibrium (static) density profile and another one for the study of the time-dependent perturbation around the equilibrium equilibrium state. But within this approach, one can use more sophisticated free energy density functionals that were recently developed in orbital-free DFT [72,73,81] for the calculation of the equilibrium density distribution, taking into account correlation effects. As the next step one can use the general expression (72) or the approximation discussed in Sec. III.E for the consideration of the time-dependent density perturbation. The only question remaining is how to include in the most consistent way correlation effects into the QHD description of the density perturbation of arbitrary frequency. This question is discussed in the following section.

V. EXCHANGE-CORRELATION POTENTIAL FOR QHD APPLICATION
Now we make progress in another direction: we include, in addition to the ideal free energy, also the exchange-correlation free energy functional F xc . This will allow us to generalize the polarization function from the ideal to the interacting case, Π id QHD → Π QHD . We note that for T = 0 exchange correlation contributions were included in QHD phenomenologically in Ref. [42].

A. QHD and local field corrections
It is known from linear response theory and DFT that the exchange-correlation free energy can be obtained from the local field correction [73,74]. Now we use the same approach to derive the exchange-correlation potential for the QHD application.
Taking into account the result from Eq. (72) for the second order functional derivative of F id [n], we find from Eq. (71): Interestingly, this expression has the same form as the density response function of a correlated electron gas, e.g. [75] or the polarization function, within the formalism of local field corrections [74]: withũ(k) = 4πe 2 /k 2 being the Fourier transform of the Coulomb potential. From the requirement that the correlated QHD polarization is in agreement with the latter result, i.e. Π QHD (k, ω) ≡ Π LFC (k, ω), we now have the possibility to derive the exchange-correlation free energy in terms of local field corrections, for application in the QHD equation (68): Note that a similar result was obtained in the framework of TDDFT [19,23]. The analytic properties of the dynamic local field correction, such as the asymptotic expansion, were studied in detail by Kugler [76]. The interpolation formula for the dynamic local field correction was considered in Refs. [77] and [78]. The static local field correction G(k) can be obtained using the finite temperature STLS approximation [79][80][81] or by ab initio quantum Monte Carlo simulations [4]. In the latter case, G(k) can be represented by a fit formula for small and large wave numbers [82,83], , where the coefficients A and B can be obtained using analytical fits for the exchangecorrelation free energy per electron, f xc , and the pairdistribution function of the uniform electron gas, g(r), at r → 0, based on quantum Monte Carlo data [73]. For the ground state, θ → 0, an accurate analytical formula of G(k) which has been fitted to quantum Monte Carlo data [84] was presented by M. Corradini et al. [85]. For the finite temperature case, an accurate parametrization of the exchange-correlation free energy has been provided recently by Groth et al. [86], and first ab initio calculations of the local field correction were presented in Ref. [88].
Further, knowing the static local field correction, the dynamic result for G(k, ω) can be calculated, for instance, on the basis of the method of moments [90]. Alternatively, the dynamic STLS approximation can be used to calculate G(k, ω) for both ground state [91] and at finite temperature [92,93].

B. Collision effects in relaxation-time approximation
In order to illustrate the effect of correlations in the most simple way, let us consider the polarization function in the relaxation time approximation (Mermin polarization function) [54,94]: with the definitioñ and ν being the electron collision frequency. Comparison of Eqs. (80) and (82), leads to: .
In the long-wavelength limit, k → 0, we use the expansion (27), and derive from Eq. (84) whereã 0 [n 0 ] andã 2 [n 0 ], are the (frequency-dependent) expansion coefficients of the RPA polarization in the long-wavelength limit that were given in tables I and II, andã 0 0 [n 0 ] andã 0 2 [n 0 ] are the respective zero-frequency limits.
This example corresponds to the simplest form of the dynamic exchange-correlation potential. However, this approximation appears to be very useful for the description of dense plasmas and warm dense matter if one extends the model to a dynamic collision frequency, ν = ν(ω) [96][97][98][99]. In this case, ν(ω) can be computed taking into account quantum and non-ideality effects by other techniques such as quantum kinetic theory or Green functions [96], or semiclassical molecular dynamics simulations [100]. Now we can write down the general form of the QHD momentum equation: where the notation n 0 → n 0 (r) emphasizes that the initially constant equilibrium density is allowed to vary in space, but only after performing the inverse Fourier transformation to real space.

VI. SUMMARY
In this paper, first of all, previously known results on QHD have been revisited and extended. To this end we started with a general expression for the free energy as a functional of the density and used the local density approximation together with gradient corrections. Explicit results were obtained in different limiting cases and compared to the RPA polarization function. The factorsᾱ and γ in the equation for the Fermi pressure and the Bohm potential have been derived consistently connecting the linear density response function in the RPA to the Thomas-Fermi theory complemented by the first order density gradient correction This goes substantially beyond many previous works where these pre-factors were included empirically by modifying the equation of state of the ideal electron gas [30,35] and the Bohm potential in order to reach agreement between the QHD and the RPA results for the plasmon dispersion.
Secondly, a generalized non-local Bohm potential was derived in linear response and linked to the RPA polarization function via the second-order functional derivative of the non-interacting free energy density, Eq. (72). This has allowed us to avoid the gradient expansion entirely and, thus, constitutes a crucial step in the further improvement of QHD. In fact, this approach has allowed us subsequently to systematically include exchange-correlation contributions (terms beyond RPA) of the free energy.
Using the obtained relation, one possible non-local form of the Bohm potential was proposed in Eq. (76), making use of the ansatz Eq. (73). It is worth noting that, on the basis of Eq. (72), one may find different forms of the non-local Bohm potential by utilizing different approximations for the non-interacting free energy density functional. In the static case, as it is known from orbital-free density functional theory, there exist several different apporximations, based on the relation between the second-order functional derivative of the free energy and the inverse RPA polarization func-tion, such as a density-dependent (or independent) kernel [101], and the so-called two-point and single-point functionals [102,103]. However, any choice of an ansatz must be checked for consistency and numerical stability [104] to avoid un-physical results.
As a third result, the exchange-correlation potential for the QHD application in linear response has been analyzed and expressed in terms of the dynamic local field correction. This has allowed us to use the result of previous studies of the dynamic local field correction in the QHD theory. In the present paper, the Bohm and exchangecorrelation potentials were first considered for the case of the uniform electron gas and, after determining the potentials, for the case of a spatially varying equilibrium density. This means that, regarding an equilibrium density variation in space, the result was obtained in the adiabatic approximation. From TDDFT, it is known however, that, because of the long-ranged behavior of the exchange-correlation potential of an inhomogeneous electron system, a frequency-dependent adiabatic local density approximation does not exist [19,23]. In other words, the exchange-correlation kernel (the second-order functional derivative of F xc ) is not a short-ranged function of |r − r |. This feature is known as the ultra-nonlocality problem of time-depended DFT. However, the exchange-correlation potential in LDA can be used if the characteristic scale on which the equilibrium density distribution changes is much larger than that of the timedependent potential [19,105] as the exchange-correlation kernel of the homogeneous system is a short-ranged function of |r − r |.
Information about the applicability of QHD for the description of electrons with a non-uniform equilibrium density distribution can be deduced by considering the continuity equation (67) ∂n 1 (r, t) ∂t = ∇(n 0 ∇w 1 ) = n 0 ∇ 2 w 1 1 + ∇n 0 n 0 When one turns to the case of a uniform electron gas, the information about the term in brackets on the right hand side is lost. This means that the obtained result is valid only if ∇n0 n0 ∇w1 ∇ 2 w1 1. Reformulating this condition in terms of the velocity, we obtain: The condition (88) means that the length scale of the equilibrium density variation must be much larger than the length scale of the velocity variation. Therefore, the QHD model under consideration is designated for description of quantum electrons with a weakly inhomogeneous equilibrium density distribution, but possibly, with strong correlations.
Another important point is that it is straightforward to incorporate a magnetic field into the QHD model via the minimal coupling approach, i.e. by redefining the velocity [or the scalar field w, in Eq. (8)] as −∇w(r, t) = v(r, t) → v(r, t) − 1/c A(r, t), e.g. [106], where A is the vector potential.
Finally, we note that the presented recipe for the consistent derivation of the quantum potential can be used for formulation of a QHD model for electrons confined in lower dimensions.

ACKNOWLEDGMENTS
We are grateful to T. Dornheim and S. Groth for helpful discussion on local field corrections. Zh.A. Moldabekov gratefully acknowledges funding from the German Academic Exchange Service (DAAD) under program number 57214224. This work has been supported by the Ministry of Education and Science of Kazakhstan under Grant No. 0263/PSF (2017).

APPENDIX A
For the partial derivative of a Hamiltonian with respect to a parameter a, the proof of the relation ∂H ∂a = ∂Ω ∂a is given in Ref. [51], where . . . denotes averaging over the grand canonical ensemble, and H is a Hamiltonian of the system at rest and in the absence of an external field, i.e w = 0 and V ext = 0. Here we extend this relation to the case of the functional derivative δH δn = δE δn . We consider the density distribution as the sum of the mean density n 0 = const and the modulation n 1 (r) around n 0 , with n 1 (r)dr = 0. The mean density n 0 is understood as the smoothed density distribution with respect to microscopic density fluctuations. The semi- n 1 (r ) + . . . , (90) where δn = δn 1 . Now we take into account that, for a functional of the form F [f ] = exp f (x)g(x)dx , the functional derivative reads and use Eqs. Here the convergence of the expansion (27) for the static case is discussed. This is of interest for the description of the screening of an ion by electrons as well as in the context of the construction of the non-interacting free energy density functional for orbital-free DFT applications [45].
In Fig. 4, curves of the inverse static polarization function in the RPA in units of its value at k = 0, Π −1 RPA (0, 0) = 2ã 0 , are shown, for different values of the degeneracy parameter. At a fixed wavenumber, k/2k F > 1, the dimensionless inverse static RPA polarization function monotonically decreases with θ, as is demonstrated in Fig. 4. In contrast, in the long wavelength limit, k/2k F < 1, the dependence of Π −1 RPA (k, 0) on θ is non-monotonic. In this limit, with increase in the degeneracy parameter up to ∼ 0.75 the value of Π −1 RPA (k, 0)/Π −1 RPA (0, 0) increases. In contrast, at θ 0.75, the increase in the degeneracy parameter leads to the decrease of the dimensionless Π −1 RPA (k, 0). Further, we discuss the long-wavelength limit. In Fig. 5, the convergence of the expansion, Eq. (27), for the static case, Π −1 RPA (k, 0) = 2 ã l k l with the coefficientsã l given in Ref. [45], is demonstrated (where all odd terms are equal to zero). As is seen from Fig. (5), already the first non-zero correction, with l = 2, gives a good description of Π −1 RPA (k, 0), at k/k F < 1. The agreement of the expansion with the exact curve becomes better with increase in θ. From this one can conclude that, at θ ∼ 1 and k 2k F , taking into account the first order correction, 2ã 2 k 2 , gives already a good description of the static polarization function, at least in the case of the homogeneous electron gas. For θ 1, the accurate interpolation of the inverse static polarization function in the RPA at k < 2k F is provided by taking into account the second non-zero correction (l = 4): Π −1 RPA (k, 0) = 2(ã 0 +ã 2 k 2 +ã 4 k 4 ).