Quantum hydrodynamics (QHD) theory for finite temperature plasmas is consistently derived in the framework of the local density approximation of the free energy with 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 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–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–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 (non-ideality) 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}

where $E$ is the electric-field strength, and $\Sigma $ stands for a tensor that contains all many-body and quantum effects, including correlations and dissipation. The exact form of $\Sigma $ is not known, and different approaches and approximations for $\Sigma $ have been proposed with different levels of sophistication, e.g., Refs. 24 and 25.

An alternative concept are first-principles approaches based on wave function methods, e.g., Ref. 26, quantum-statistical 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–34} plasmonics,^{8–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 single-electron 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–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–41 and 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 non-local 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}

where *w* is the scalar potential determining the velocity field by $ v = \u2212 \u2207 w , \u2009 E [ n ] = E id [ n ] + E xc [ n ] $ 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}

with the definition

and making use of relations $ v = \u2212 \u2207 w $ and $ ( v \xb7 \u2207 ) v = 1 2 \u2207 ( \u2207 w ) 2 $, we arrive at the QHD equation in terms of average density $ n ( r , t ) $, velocity $ v ( r , t ) $, and the generalized force $ \u2212 \u2207 \mu [ r , t ] $^{46–48}

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 = \u222b n ( r ) d r $ corresponds to the mean value of the number of particles in the grand canonical ensemble. The derivation of Eq. (13) is given in 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

where $ F [ n ] = F id [ n ] + F xc [ n ] $ is decomposed into the ideal (non-interacting) part $ F id [ n ] $ and the exchange correlation part $ F xc [ n ] $. In the static long wavelength limit (see below), the generalized force, $ \u2212 \u2207 \mu [ n ( r , t ) , T ] $, is related to the pressure tensor $ P \xaf \xaf $^{52}

which provides the link with the standard fluid theory, cf. Eq. (2).

## III. DERIVATION OF THE BOHM POTENTIAL AND FERMI PRESSURE

### A. General expressions

In order to derive a potential related to the Fermi pressure and the Bohm potential, we neglect the exchange-correlation 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 $ \delta n \u0303 $ and $ \delta w \u0303 $ denote the Fourier transforms of *δn* and *δw*, respectively. The resulting variation of $ \mu id $, as defined by Eq. (19), has the form^{55}

Thus, the linearized momentum equation (12) for $ w \u0303 = \delta w \u0303 $ and, taking into account Eq. (21), is written as

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, $ \Pi \u2212 1 \u2261 e \delta \phi \u0303 / \delta n \u0303 $

### 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 = \omega / ( k v F ) , \u2009 z = k / ( 2 k F ) $, and defined $ \chi 0 2 = ( \pi k F a B ) \u2212 1 \u2243 r s / 6.03 $, where *a _{B}* is the Bohr radius, $ r s \u2261 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 \pi 2 n ) 1 / 3 , \u2009 \omega p 2 = 4 \pi n e 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

The convergence of the expansion (27) is discussed in Appendix B for the static case and *k* < 2*k _{F}*.

We stress that the coefficients $ a \u0303 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 for the ground state was given in Ref. 19. For completeness of the consideration, the proof of the *stiffness* theorem at a finite temperature for the case of the density perturbation by an external electric field and the discussion of the consistency of the results of this paper with this theorem are given in Appendix C.

Once the coefficients $ a \u0303 0 [ n 0 ] $ and *a*_{2}[*n*_{0}] are defined, we allow the equilibrium density to vary in space and time, $ n 0 \u2192 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 $ a \u0303 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 temperature, *β* = 1/(*k _{B}T*).

Using the local density approximation, we obtain $ n ( r ) = 2 m 3 / 2 \pi 2 \beta 3 / 2 \u210f 3 I 1 / 2 [ \eta ( r ) ] $. The following partial derivatives are needed to find the Bohm potential (31):

where $ C = 2 \pi 2 \beta 3 / 2 \u210f 3 2 m 3 / 2 = 4 3 n \theta \u2212 3 / 2 $. As is shown below, at $ \theta = k B T / E F \u223c T \xd7 n \u2212 2 / 3 \u2192 0 $, we have $ a 2 [ n ] = \gamma \u210f 2 / ( 8 m n ) $, 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 $ \theta \u2192 0 $, can be written as

where the coefficient $ \alpha \xaf $ 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,62} for zero temperature (note that $ F \u2192 E $, at $ \theta \u2192 0 $).

Now, we will use Eq. (29), to analyze different limits for the coefficients $ a \u0303 0 [ n ] $ and $ a 2 [ n ] = \u2212 a \u0303 2 $, for different frequency-wavenumber ranges.

### C. Low frequency and long wavelength limit, $ \omega \u226a \u210f k 2 / 2 m , k \u226a 2 k F $

The results given for the static case remain valid also for low frequencies, i.e., $ u \u226a 1 $ or $ \omega \u226a k v F $. This regime is relevant for a variety of physical processes such as the screening of a test charge and the dispersion of low-frequency 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 was treated incorrectly using the Bohm potential for the high-frequency case rather than for the static case. This has led to the unphysical prediction of ion-ion attraction in an equilibrium quantum plasma; for a discussion, see Refs. 40 and 45. To obtain the Fermi pressure and Bohm potential in the present limiting case, we calculate the coefficients $ a \u0303 0 [ n ] $ and $ a 2 [ n ] $ and obtain^{45}

with $ H 1 ( \eta ) = \theta [ n ] 2 \u2009 I \u2212 1 / 2 ( \eta ) , \u2009 \u2003 H 2 ( \eta ) = 1 2 \theta [ n ] \u2009 I \u2212 3 / 2 ( \eta ) $. In the low-temperature limit, $ \theta \u2192 0 $, we have $ H 1 ( \eta ) \u2243 1 $ and $ H 2 \u2243 \u2212 1 $, and Eq. (37) gives the asymptotic results $ a 2 \u2192 \u210f 2 / ( 72 m n ) $ and $ \gamma = 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 finite-temperature generalization of the Bohm potential can be easily found by substituting Eq. (37) into Eq. (31) and taking into account the dependence $ \theta [ n ] \u223c n \u2212 2 / 3 $. Furthermore, for the static case (*ω* = 0), we have $ f 0 \u2261 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, $ \u2212 1 2 \u2202 2 f T F [ n ] \u2202 n 2 $, in the static long wave length limit ( $ k \u226a 2 k F $) exactly coincides with $ a \u0303 0 [ n ] $ given by Eq. (36). Finally, the functional derivative of the Thomas-Fermi term yields

At zero temperature, *η* = *βE _{F}*, leading to

*δ T*

_{TF}/

*δn*=

*E*, in agreement with the result of Ref. 48.

_{F}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 3 n \theta \u2212 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 $ \gamma \xaf $ of the first term in the Bohm potential for the finite-temperature case (47) and the dependence of *γ* and $ \gamma \gamma \xaf $ on *θ* are presented in Fig. 1. In the limit of fully degenerate electrons, $ \theta \u2192 0 $, as well as in the classical limit, $ \theta \u226b 1 $, the correction coefficient approaches unity, $ \gamma \xaf \u2192 1 $. This correction coefficient $ \gamma \xaf $ is important for a partially degenerate plasma, around $ \theta \u223c 0.5 $, as can be seen from Fig. 1.

Recently, Haas and Mahmood^{60} considered the finite-temperature Bohm potential by analyzing ion-acoustic 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 $ \gamma \xaf $ 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 $ \u223c ( n 0 / n 1 ) 2 $, where *n*_{1} is a small density perturbation (i.e., $ n 0 / n 1 \u226a 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 $ \gamma \xaf $ is lost.

### D. Short wavelength limit, $ k \u226b 2 k F $ at low frequencies, $ \omega \u226a \u210f k 2 / 2 m $

For the degenerate electron gas, Jones and Young^{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 ] = \u210f 2 8 m e n $, 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 \u226a z , z \u226b 1 $, we use the following expansion of the function $ \Delta g = g ( u + z ) \u2212 g ( u \u2212 z ) $:^{56}

where

which can be substituted into the formula for the inverse polarization function, with the result

From this, we obtain the coefficients $ a \u0303 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, $ \theta \u226a 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 $ \alpha \xaf $, 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 / 2 k F \u226a 1 $, and short wavelengths, $ k / 2 k F \u226b 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 \u2243 k F $.

### E. High frequency limit, $ \omega \u226b \u210f k 2 / 2 m $

The present limiting case is important for a variety of physical situations, such as for the description of the optical (Langmuir) 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 \u226b z $, we use the formulas^{56}

and obtain

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 $ \u2212 \u222b 2 a \u0303 0 [ n ] d n $. In the zero-temperature limit, taking into account the relation $ I 3 / 2 ( \theta \u2192 0 ) = 2 5 \theta \u2212 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 $ \alpha \xaf = 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 $ \alpha \xaf $ in the *k*–*ω* plane are summarized in Fig. 3. Furthermore, the information about the obtained coefficients $ a 2 [ n ] , \gamma $ and $ a \u0303 0 [ n ] $, $ \alpha \xaf $ is listed in Tables I and II.

. | a_{2}[n]
. | γ
. |
---|---|---|

$ \omega \u226a \u210f k 2 / 2 m e $ | $ \u2212 \u210f 2 H 2 ( \eta ) 72 m e n H 1 2 ( \eta ) $ | 1/9 |

$ k \u226a 2 k F $ | ||

$ \omega \u226a \u210f k 2 / 2 m e $ | $ \u210f 2 8 m e n $ | 1 |

$ k \u226b 2 k F $ | ||

$ \omega \u226b \u210f k 2 2 m e $ | $ \u210f 2 8 m e n $ | 1 |

. | a_{2}[n]
. | γ
. |
---|---|---|

$ \omega \u226a \u210f k 2 / 2 m e $ | $ \u2212 \u210f 2 H 2 ( \eta ) 72 m e n H 1 2 ( \eta ) $ | 1/9 |

$ k \u226a 2 k F $ | ||

$ \omega \u226a \u210f k 2 / 2 m e $ | $ \u210f 2 8 m e n $ | 1 |

$ k \u226b 2 k F $ | ||

$ \omega \u226b \u210f k 2 2 m e $ | $ \u210f 2 8 m e n $ | 1 |

. | $ a \u0303 0 [ n ] $ . | $ \alpha \xaf $ . |
---|---|---|

$ \omega \u226a \u210f k 2 2 m e $ | $ \u2212 \pi e 2 2 k F 2 \chi 0 2 H 1 ( \eta ) $ | 1 |

$ k \u226a 2 k F $ | ||

$ \omega \u226a \u210f k 2 2 m e $ | $ \u2212 4 \pi e 2 k F 2 \chi 0 2 \xd7 3 16 I 3 / 2 ( \eta ) \theta 5 / 2 $ | $ 3 5 $ |

$ k \u226b 2 k F $ | ||

$ \omega \u226b \u210f k 2 2 m e $ | $ \u2212 4 \pi e 2 k F 2 \chi 0 2 \xd7 9 16 I 3 / 2 ( \eta ) \theta 5 / 2 $ | $ 9 5 $ |

. | $ a \u0303 0 [ n ] $ . | $ \alpha \xaf $ . |
---|---|---|

$ \omega \u226a \u210f k 2 2 m e $ | $ \u2212 \pi e 2 2 k F 2 \chi 0 2 H 1 ( \eta ) $ | 1 |

$ k \u226a 2 k F $ | ||

$ \omega \u226a \u210f k 2 2 m e $ | $ \u2212 4 \pi e 2 k F 2 \chi 0 2 \xd7 3 16 I 3 / 2 ( \eta ) \theta 5 / 2 $ | $ 3 5 $ |

$ k \u226b 2 k F $ | ||

$ \omega \u226b \u210f k 2 2 m e $ | $ \u2212 4 \pi e 2 k F 2 \chi 0 2 \xd7 9 16 I 3 / 2 ( \eta ) \theta 5 / 2 $ | $ 9 5 $ |

Using the high-frequency result for $ a \u0303 0 $ and the Bohm potential, one can derive the well-known plasmon dispersion from the continuity and momentum equations

which, at $ \theta \u226a 1 $, approaches the form

where we took into account that $ 2 a \u0303 0 [ n 0 ] n 0 / m e \u2192 \u2212 3 5 v F 2 $, at $ \theta \u2192 0 $, and that, in the high-frequency limit, $ a 2 [ n ] = \u210f 2 / 8 m e n $.

On the other hand, in the limit $ \theta \u226b 1 $, taking into account that $ a \u0303 0 [ n 0 ] \u2192 \u2212 3 2 k B T n 0 $ leads to the dispersion relation

where $ v th = k B T m e $ 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, $ \Pi 0 ( \omega ) $, Eq. (28), without any additional terms related to the fermionic pressure.

## IV. GENERALIZED NON-LOCAL BOHM POTENTIAL IN LINEAR RESPONSE

In Sec. III, 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, $ \Pi QHD id ( k , \omega ) $, 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

where *μ*_{0} is a constant. The equilibrium configuration is now exposed to a weak external perturbation giving rise to $ n ( r , t ) = \u2009 n 0 ( r ) + n 1 ( r , t ) , \u2009 w ( r , t ) = w 1 ( r , t ) , \u2009 \phi ( r , t ) = \phi 0 ( r ) + \u2009 \phi 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}

where *w*_{1} determines the velocity perturbation via $ v 1 = \u2212 \u2207 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 wavenumber $ ( k , \omega ) $ space, which we denote by $F$

From Eqs. (69) and (70), immediately follows the QHD result for the inverse polarization function $ \Pi QHD \u2212 1 ( k , \omega ) = e \phi \u0303 1 / n \u0303 1 $

Neglecting the exchange-correlation contribution to the free energy and taking $ \Pi QHD = \Pi RPA $, it is straightforward to deduce from Eq. (71) that

where we used the definition (28) of the long-wavelength limit, $ \Pi 0 ( \omega ) $, 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 Sec. III. 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}

where the kernel *K* is symmetric with respect to permutation of $r$ and $ r \u2032 $. Using Eqs. (18) and (29), Eq. (73) can be written as

where $ K \u0303 ( [ 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 non-interacting free energy density functional.^{45,67–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 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 Sec. V.

## 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, $ \Pi QHD id \u2192 \Pi 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., Ref. 75 or the polarization function, within the formalism of local field corrections^{74}

with $ u \u0303 ( k ) = 4 \pi 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., $ \Pi QHD ( k , \omega ) \u2261 \Pi LFC ( k , \omega ) $, 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 Singwi-Tosi-Land-Sjoelander (STLS) approximation^{79–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} $ G ( k ) = A [ 1 \u2212 exp \u2009 ( \u2212 B k 2 ) ] $, where the coefficients *A* and *B* can be obtained using analytical fits for the exchange-correlation free energy per electron, *f*_{xc,} and the pair-distribution function of the uniform electron gas, *g*(*r*), at $ r \u2192 0 $, based on quantum Monte Carlo data.^{73} For the ground state, $ \theta \u2192 0 $, an accurate analytical formula of *G*(*k*) which has been fitted to quantum Monte Carlo data^{84} was presented by 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.^{87,89}

### 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 definition

where $ a \u0303 0 [ n 0 ] $ and $ a \u0303 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 $ a \u0303 0 0 [ n 0 ] $ and $ a \u0303 2 0 [ n 0 ] $ are the respective zero-frequency limits.

For the exchange-correlation term in the momentum equation (70), the relaxation time approximation gives rise to a friction force

where $ w \u0303 1 = i \omega k 2 n 0 n \u0303 1 $, and Eqs. (81) and (85) were used. If one retains also terms scaling as $ \u223c a \u0303 0 $, and $ \u223c a \u0303 2 $ in Eq. (85), the so-called hydrodynamic Drude model used in plasmonics^{35,43,95} is recovered. The electron collision frequency, *ν*, in the Mermin polarization function (82) is due to all scattering processes^{20} and must incorporate a contribution from electron-electron collisions^{96,97} as well as electron-ion (neutral) collisions.

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, *ν* = *ν* (*ω*).^{98–101} 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^{98} [for electron-ion (neutral) collisions], or semiclassical molecular dynamics simulations.^{102} Now we can write down the general form of the QHD momentum equation

where the notation $ n 0 \u2192 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 AND OUTLOOK

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 $ \alpha \xaf $ 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 approximations, based on the relation between the second-order functional derivative of the free energy and the inverse RPA polarization function, such as a density-dependent (or independent) kernel,^{103} and the so-called two-point and single-point functionals.^{104,105} However, any choice of an ansatz must be checked for consistency and numerical stability^{106} 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 exchange-correlation 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 \u2212 r \u2032 | $. This feature is known as the ultra-non-locality 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 time-dependent potential^{19,107} as the exchange-correlation kernel of the homogeneous system is a short-ranged function of $ | r \u2212 r \u2032 | $.

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)

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 $ | \u2207 n 0 \xb7 \u2207 w 1 n 0 \u2207 2 w 1 | \u226a 1 $. Reformulating this condition in terms of the velocity, $ | \u2207 n 0 \xb7 v | \u226a n 0 | \u2207 \xb7 v | $, 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 $ \u2212 \u2207 w ( r , t ) = v ( r , t ) \u2192 v ( r , t ) \u2212 1 / c \u2009 A ( r , t ) $, e.g., Ref. 108, 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 No. 57214224. This work has been supported by the Ministry of Education and Science of Kazakhstan under Grant No. 0263/PSF (2017).

### APPENDIX A: FUNCTIONAL DERIVATIVE OF HAMILTONIAN IN GRAND CANONICAL ENSEMBLE

For the partial derivative of a Hamiltonian with respect to a parameter *a*, the proof of the relation $ \u3008 \u2202 H \u2202 a \u3009 = \u2202 \Omega \u2202 a $ is given in Ref. 51, where $ \u3008 \u2026 \u3009 $ 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 $ \u3008 \delta H \delta n \u3009 = \delta E \delta 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 $ \u222b n 1 ( r ) d r = 0 $. The mean density *n*_{0} is understood as the smoothed density distribution with respect to microscopic density fluctuations. The semi-classical Hamiltonian can be expanded around *n*_{0} as

from where we see that

where *δn* = *δn*_{1}.

Now we take into account that, for a functional of the form $ F [ f ] = exp \u2009 [ \u222b f ( x ) g ( x ) d x ] $, the functional derivative reads

where *T* is in energy units (*k _{B}* = 1).

Finally, with the help of (A4), we have

### APPENDIX B: CONVERGENCE OF THE EXPANSION OF THE RPA POLARIZATION FUNCTION

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, $ \Pi RPA \u2212 1 ( 0 , 0 ) = 2 a \u0303 0 $, are shown, for different values of the degeneracy parameter. At a fixed wavenumber, *k*/2*k _{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*/2

*k*< 1, the dependence of $ \Pi RPA \u2212 1 ( k , 0 ) $ on

_{F}*θ*is non-monotonic. In this limit, with increase in the degeneracy parameter up to ∼0.75 the value of $ \Pi RPA \u2212 1 ( k , 0 ) / \Pi RPA \u2212 1 ( 0 , 0 ) $ increases. In contrast, at $ \theta \u2273 0.75 $, the increase in the degeneracy parameter leads to the decrease in the dimensionless $ \Pi RPA \u2212 1 ( k , 0 ) $.

Further, we discuss the long-wavelength limit. In Fig. 5, the convergence of the expansion, Eq. (27), for the static case, $ \Pi RPA \u2212 1 ( k , 0 ) = 2 \u2211 a l \u0303 k l $ with the coefficients $ a l \u0303 $ 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 $ \Pi RPA \u2212 1 ( 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 \u226a 2 k F $, taking into account the first order correction, $ 2 a 2 \u0303 k 2 $, gives already a good description of the static polarization function, at least in the case of the homogeneous electron gas. For $ \theta \u2273 1 $, the accurate interpolation of the inverse static polarization function in the RPA at *k* < 2*k _{F}* is provided by taking into account the second non-zero correction (

*l*= 4): $ \Pi RPA \u2212 1 ( k , 0 ) = 2 ( a \u0303 0 + a \u0303 2 k 2 + a \u0303 4 k 4 ) $.

### APPENDIX C: STIFFNESS THEOREM AT FINITE TEMPERATURE

In this appendix, the relation between the static inverse density response function and a change in the free energy due to the density perturbation by an external field is given. Suppose we have a unperturbed free energy *F*_{0}[*n*_{0}], with the mean density *n*_{0} = const. Applying a weak external field, *V*_{ext,} the electron density becomes $ n = n 0 + n \u0303 $, where $ | n \u0303 | / n 0 \u226a 1 $. The total *intrinsic free energy* of the perturbed system reads

where the expansion of the *intrinsic free energy* in terms of the density perturbation $ n \u0303 $ is used [cf. Eq. (89)]. As the result of conditions $ \u222b n \u0303 ( r ) d r = 0 $ and $ \delta F [ n ] \delta n | n 0 = const $, there is no contribution in Eq. (C1) due to the term linear in $ n \u0303 $. Additionally, the potential energy $ U [ n \u0303 ] $ in Eq. (C1) (terms in parentheses on the l.h.s.) does not explicitly depend on *n*_{0}, due to quasi-neutrality condition, i.e., we assume that the system has a positive neutralizing homogeneous background which remains unchanged by the external electric field. Therefore, it is taken into account that $ \delta U [ n \u0303 ] \delta n | n = n 0 \u2261 0 $. The one-to-one correspondence between the external potential and the equilibrium density was proven by Mermin.^{54} Further, we skip the integral with the kernel *L* and higher terms which are related to nonlinear response features.

From Eq. (C1), the minimization condition of *F*[*n*] with respect to $ n \u0303 $ allows us to identify

where it is taken into account that *F*_{0}[*n*_{0}] is constant, and the effective potential is defined as

Now, using the Fourier expansion $ K ( r \u2212 r \u2032 ) = 1 \Omega \u2211 k K \u0303 ( k ) \u2009 exp \u2009 [ \u2212 i k \xb7 ( r \u2212 r \u2032 ) ] $ and the convolution theorem, from Eq. (C2) we find

Finally, substituting Eqs. (C2) and (C4) into Eq. (C1) we arrive at the desired relation between the inverse density response function and a change in the free energy due to the density perturbation by the external field

where the density response function is defined as $ \chi = n \u0303 ( k ) / [ e V \u0303 ext ( k ) ] $ and/or in terms of the polarization function $ \Pi = n \u0303 ( k ) / [ e \phi \u0303 eff ( k ) ] $, and the Fourier transform of Coulomb potential, $ u \u0303 = 4 \pi e 2 / k 2 $, as

Equation (C5) has the meaning of a *stiffness theorem at finite temperature* for the case of the density perturbation by the external electric field. The term on r.h.s. of Eq. (C5) is referred to as *stiffness energy*, which proportional to the square of the amplitude of the density deviation from the original equilibrium value, $ | n \u0303 | 2 $, and to the *stiffness* of the system (represented by *χ*^{–1}).

The relation (C4) for the ground state was obtained in Ref. 21. For the ground state, $ F \u2192 E $, the relation (C5) was used to compute the static response function by quantum Monte Carlo in Ref. 84. For the discussion of higher terms in Eq. (C1), we refer to Refs. 68, 69, and 109. Note that *F*[*n*] and *F*_{0}[*n*_{0}] in Eq. (C1) are the exact total free energies after and before application of the external field, respectively.

Now we show that the results of this paper on the Bohm potential are consistent with relations (C4) and (C5). Indeed, in the static case, Eq. (75) agrees with Eq. (C4) (if one considers non-interacting electrons) except for the constant summand $ a \u0303 0 [ n 0 ] $ which does not contribute to the total free energy since the volume integral of $ n \u0303 $ vanishes. This constant term appeared because *F*_{0}[*n*] in Eqs. (17) and (73) is the LDA approximation to the already perturbed system's free energy, in contrast to *F*_{0}[*n*_{0}] in Eq. (C1).

## References

In Ref. 45, the notation *T*[*n*] for the noninteracting free energy density functional was used instead of *F*_{id}[*n*].

_{xc}(r) of the homogeneous electron gas

^{–}– e

^{–}collisions