Generalized hydrodynamics revisited

During the past decade a number of attempts to formulate a continuum description of complex states of matter have been proposed to circumvent more cumbersome many-body and simulation methods. Typically these have been quantum systems (e.g., electrons) and the resulting phenomenologies collectively often called"quantum hydrodynamics". However, there is extensive work from the past based in non-equilibrium statistical mechanics on the microscopic origins of macroscopic continuum dynamics that has not been exploited in this context. Although formally exact, its original target was the derivation of Navier-Stokes hydrodynamics for slowly varying states in space and time. The objective here is to revisit that work for the present interest in complex quantum states - possible strong degeneracy, strong coupling, and all space-time scales. The result is an exact representation of generalized hydrodynamics suitable for introducing controlled approximations for diverse specific cases, and for critiquing existing work.


I. INTRODUCTION
Traditional hydrodynamics for a simple fluid describes the dynamics of the average local conserved fields associated with symmetries of the Galilean group (number density n(r, t), energy density e(r, t), momentum density p(r, t) 1,2 . These fields are chosen since they dominate the behavior on large space and time scales, leading to a closed dynamics simpler than that of the full manybody degrees of freedom. The governing equations have their origin in the macroscopic conservation laws that follow from averages of the corresponding exact microscopic conservation laws for the operators representing these fields. The averages are defined with respect to a chosen state (ensemble) whose dynamics is governed by the Liouville-von Neumann equation. The generic forms of the macroscopic conservation laws, classical or quantum, do not depend on the specific state of interest except through the explicit forms for the average energy and momentum fluxes. The latter are "constitutive equations" relating those fluxes to the hydrodynamic variables for a closed description of their dynamics. This is analogous to thermodynamics whose formulation is general, but which requires an equation of state for each specific system. For variations on small space and time scales the Liouville-von Neumann equation can be solved to obtain the average fluxes in the form of Fourier's law and Newton's viscosity law 3,4 . The resulting hydrodynamic equations are the non-linear Navier-Stokes equations, characterized by a thermodynamic pressure and transport coefficients given formally by Green-Kubo time correlation functions 5 .
The objective here is to extend the continuum description to general complex states being studied now in various evolving fields of condensed matter physics and materials sciences. Examples of recent reviews include warm, dense matter 6 , high energy density physics 7,8 , thermoelectric transport 9 , quantum plasmas 10,11 , and electrons in graphene 12 . This entails accommodating conditions of strong coupling, strong quantum degeneracy, and/or shorter space and time scales. This idea has been explored extensively some years ago for calculating equilibrium time correlation functions 13,14 . Such functions obey the linear Navier-Stokes equations at low frequencies and long wavelengths 15,16 . These conditions are met for the differential cross section measured by laser light scattering, and the latter is a useful tool for measuring the hydrodynamic transport coefficients. However, that approach does not apply to the neutron scattering cross section which samples much smaller wavelengths and larger frequencies. To address this problem the hydrodynamic description is modified (generalized) to include these extended domains. Typically this is done by making thermodynamic derivatives wave vector dependent and transport coefficients frequency and wave vector dependent. An early review of the methods with references is given by Mountain 17 .
The scope of the objective here is far greater, to provide the basis for a continuum description of macroscopic dynamics under the most general conditions. It is a formally exact resolution of the macroscopic conservation laws into physically transparent components suitable for constructing approximate models in specific cases of interest. In this way much of the present phenomenology about "new" hydrodynamics can be critiqued and controlled. The analysis uses an application of nonequilibrium statistical mechanics developed long ago for linear response, but not limited to that case. The two earliest and most complete formulations are those of McLennan and collaborators 3,18-20 and of Zubarev and collaborators 4,21 . The present work can be viewed as a "revisiting" of that approach assuring its exact treatment of extreme quantum and non-local effects. It is hoped in this way an important largely forgotten tool will become available for new developments in the formulation of practical continuum dynamics models.
One motivation for this revisiting of generalized hydrodynamics is a growing current interest in a continuum description for the electron density in quantum devices (e.g.,nanomaterials) and other charged particle quantum systems (e.g., plasmonics). A very recent review of current activity in quantum hydrodynamics for plasmas has just appeared 10 . There is a much earlier history of quantum hydrodynamics associated with attempts to reinterpret the Schrodinger equation as an equivalent continuum description 22,23 , time dependent Thomas -Fermi models 24 , or from a macroscopic Hamiltonian 25 , each leading to different quantum potentials occuring in related quantum continuum equations. See Stanton and Murillo 26 for some resolution of these different quantum potentials. Previous quantum hydrodynamics have often focused on nearly ideal plasmas (weak coupling, weak spatial inhomogeneities). One interest here is for the quite different states of warm dense matter 6 , comprised of strongly coupled ions and electrons at solid state densities and a wide range of temperatures and spatial inhomogeneity. A direct approach by standard many-body methods is very challenging. Density functional theory (DFT) for equilibrium properties of such states has met with some significant success 27 . Applications of time dependent density functional theory (TDDFT) are still at an early stage in this context [28][29][30] . In principle the density and momentum conservation laws provide the basis for a closed description of the density if the momentum flux can be expressed as a functional of the density via TDDFT 28,31 . Extensions of this idea including energy conservation are also being considered 9 . Here, the closure of all conservation laws is provided by a representation of the fluxes in terms of a formal solution to the Liouville-von Neumann equation as a functional of the fields, non-local in both space and time.
In the next few sections the basis for an exact continuum description for the average density, internal energy density, and momentum density (equivalently, flow velocity) is reviewed. The corresponding conservation laws are not closed since the energy and momentum fluxes are not specified as functionals of the dependent variables. Instead they are defined as averages of specific operators over a solution to the Liouville-von Neumann equation. That solution is represented as a reference state plus its deviation. The reference state is taken to be local equilibrium (also known as the "relevant state" 21 ). This is a composite state representing an equilibrium system at each space and time point, constrained to have the same values for the average local conserved densities as those for the actual non-equilibrium state. If only the hydrodynamic fields were measured, the local equilibrium state would be indistinguishable from that for the actual system -it is the best choice within information entropy theory 21,[32][33][34] . This local equilibrium contribution to the fluxes gives results that are specific functionals of the hydrodynamic fields. If deviations from the local equilibrium fluxes are neglected, a closed set of equations for the fields is obtained. This is referred to as "perfect fluid hydrodynamics". The terminology is chosen since these equations have no entropy production. The perfect fluid hydrodynamics is local in time, but non-local in space. Hence there are no inherent limitations on space and time scales. Calculation of the local equilibrium energy and momentum fluxes is a problem closely related to that of DFT for the thermodynamics of a non-uniform system. Specific limits of these perfect fluid equations subsume all previous "quantum hydrodynamics", and important differences are noted. In the "local density approximation" where spatial non-locality is neglected these become the usual Euler equations of hydrodynamics extended to quantum states.
In Section IV the remainder of the energy and momentum fluxes, additions to the local equilibrium contributions, are considered. These are defined by averages over a special solution to the Liouville-von Neumann equation that is the deviation from the local equilibrium ensemble. It is not actually a solution but rather a representation of a solution in terms of the unknown macroscopic fields of interest. In this way averages of the fluxes are given as functionals of these fields, allowing the desired closure of the macroscopic conservation laws. The results are similar to those known for linear constitutive equations, with the equilibrium time correlation functions replaced by corresponding non-linear local equilibrium time correlation functions. They apply for both classical and quantum conditions, and no limitations are imposed on length and time scales. With the energy and momentum fluxes determined in this way as functionals of the hydrodynamic fields, the macroscopic conservation laws become the desired exact generalized hydrodynamic equations.
The derivation of the macroscopic local conservation laws in the next section is based on the fundamental microscopic Heisenberg dynamics for the system. This is contrasted with some earlier treatments starting from a macroscopic mean field Hamiltonian 11,25,31 . The associated approximate macroscopic Hamiltonian dynamics of the latter entails several limitations, most importantly an irrotational velocity field, a mean field free energy functional, and no dissipative contributions to the fluxes. The details of the exact microscopic derivation here are relegated to the Supplemental Material, so the next section is devoted to their averages, notational conventions, and their transformation from energy and momentum fields to the more usual internal energy density and flow velocity fields. Also, the energy and momentum fluxes are divided into a local equilibrium contribution and dissipative contribution. The interpretation of each is discussed. Finally, the objectives and results are reviewed in the last section with an outlook for future developments.
As noted above, the approach to macroscopic transport given here follows closely the earlier applications of nonequilibrium statistical mechanics 3,4,18-21 . The present work updates this, and provides easier access to that complex work in the context of complex states not experimentally accessible until recently. It retains validity for quantum effects, non-linearity, strong coupling and spatial inhomogeneity. Being a formally exact representation of a continuum description, it subsumes the phe-nomenological versions developed in recent years, providing a structure for their critique and extension. Some examples of this connection to earlier work are given in the text and conclusions section.
On conclusion of this work a new review of quantum hydrodynamics appeared 10 , with extensive critique of various approximations in the context of plasmas and related electronic systems. Consequently, no attempt will be made here to make contact with the rather complete discussion there of existing work. A primary premise there is that phenomenological quantum hydrodynamics has often been misused due to a lack of understanding of its validity conditions. The need for strong theoretical underpinnings of any approximation or application is stressed. The objective here is to provide an exact hydrodynamic formulation as the basis for introduction of controlled applications.
The presentation given here is all within the context of a one component fluid with simple constituents (e.g., charges, atoms, molecules) governed by non-relativistic quantum mechanics. However, the philosophy of constructing an exact formulation for the macroscopic fields prior to the introduction of approximations extends to diverse other systems mutatis mutandis with different fields, (e.g. polymers, mixtures, broken symmetry states) or different dynamics (e.g., relativistic, granular/active).

II. MACROSCOPIC LOCAL CONSERVATION LAWS
A. Microscopic local conservation laws For simplicity of notation only a one component system is considered here and below. A system of N identical particles (Bosons or Fermions) contained in a volume Ω enclosed by a surface has a Hamiltonian operator given by (1) Here U (|q α − q β |) is the interaction potential for the pair of particles at positions q α , q β , and v ext (q α ) is an external single particle potential. The boundary conditions are contained in the external potential by the condition v ext (q α ) → ∞ for q α ∈ . The state of the system is specified by a statistical density operator ρ, defined in terms of its components for each N particle Hilbert space, ρ N , and normalized according to The definition of ρ N is such that it includes a restriction to the appropriate subspace of symmetrized or antisymmetrized elements of the N particle Hilbert space. For a given observable represented by the operator A N in the N particle Hilbert space its expectation value in the state ρ is defined by The notation A; ρ instead of the simplified A will be used should any confusion arise. The Hamiltonian without the external potential has all the continuous symmetries of the Galilei group (time translations, spatial translation, rotations, boosts) as well as discrete symmetries of parity and time reversal. Associated with the continuous symmetries are the usual conservation laws for number, energy, linear momentum, and angular momentum. Here only point particles (point sources of force) are considered so the relevant conservation laws are those of number density, momentum density, and energy density. The operators representing these densities are, respectively where [, ] + denotes the anticommutator, [A, B] + = AB + BA. Here ∆(r) is a function localized about some small domain centered at r and normalized to unity dr ∆(r) = 1.
Often in the literature the ideal case of a delta function is considered, ∆(r) → δ(r). Exact microscopic local conservation laws follow directly from the Heisenberg dynamics for the densities in eqs. (4a) to (4c). The analysis is straightforward but lengthy. For completeness, the details are repeated in the Supplemental Material, with the results These operator equations are exact, with explicit expressions for the energy flux s(r, t) and momentum flux t ij (r, t) given in the Supplemental Material. A summation convention is adopted for repeated indices for Cartesian coordinates.

B. Macroscopic local conservation laws
The corresponding macroscopic averages of these microscopic equations (defined as in (3)) give the corresponding macroscopic local conservation laws.
The sources in (7b) and (7c) represent the average work done by the external force and the average external force density, respectively. The microscopic conservation laws are independent of the state of the system considered, while the macroscopic conservation laws are specific to the state. Still, they are quite general and apply for both mixed and pure states, classical and quantum. It remains to determine the fluxes such as to provide a closed set of equations for the average fields n(r, t), e(r, t), and p(r, t).
Instead of the average momentum density, it is traditional to consider the flow velocity field u(r, t) defined by p(r, t) ≡ mn(r, t)u(r, t) .
Also, it is useful to extract the purely convective parts of the energy density and fluxes of energy and momentum 3 e(r, t) = e 0 (r, t) The operators with a subscript 0 are the same as those without the subscript, except with all particle momentum operators replaced by p α → p α − m u(r, t). Hence the corresponding averages are those quantities in the local rest frame of the fluid at point r. For example, e 0 (r, t) is the average local internal energy density. The local conservation laws become where D t = ∂ t + u · ∇ is the material derivative. The dependent fields are now n(r, t), e 0 (r, t), and u(r, t). The utility of these exact equations depends on determination of the energy and momentum fluxes, s 0 (r, t) and t 0ij (r, t), in the local rest frame.

C. Analysis of energy and momentum fluxes
For a closed set of equations, the required local rest frame fluxes should be expressed as functionals of the fields so that the macroscopic conservation laws would become a closed, deterministic set of equations for the fields. That is the final objective of the presentation here. It is accomplished in two steps. The first is to identify the result for a "perfect fluid" without dissipation. Next the remainder responsible for dissipation is obtained from a formal solution to the Liouville-von Neumann equation.
Note that the time dependence can be shifted to the state ρ using the cyclic invariance of the trace s 0 (r, t) = s 0 (r); ρ(t) , t 0ij (r, t) = t 0ij (r); ρ(t) , which in turn requires the solution to the Liouville-von where L N is the Liouville operator defined by for each N , for any operator X. The solution is written as the sum of a reference state ρ ℓ N [y(t)] plus its deviation ∆ N (t) The reference state is entirely determined by a set of conjugate fields {y(t)} in one-to-one correspondence with the macroscopic conserved fields. This correspondence is defined by the requirements n ℓ (r|y(t)) ≡ n(r, t) , p ℓ (r|y(t)) ≡ mn(r, t)u(r, t) , where the superscript ℓ denotes a reference ensemble average, A ℓ = A; ρ ℓ . The left sides of these equations are functionals of the conjugate fields while the right sides are the fields of the local conservation laws. In this way the conjugate fields {y(t)} are functionals of the average conserved fields n(r, t), e 0 (r, t), u(r, t), and vice versa. The reference state therefore has the exact average values for the conserved fields by construction. A choice for ρ ℓ with these properties is the local equilibrium ensemble 3,21 . To motivate it, note that the grand ensemble for a system at rest is where by normalization The corresponding result for a system moving with velocity w is obtained by the replacement (transformation to the moving frame) H → H − w · P + 1 2 mN w 2 , where P is the total momentum operator. Now consider an assembly of equilibrium systems where β, ν, w vary at each point in space and time where The operators N, H, P have been replaced accordingly by the associated operators for their densities In this notation the local equilibrium ensemble becomes where summation over repeated indices is implied. The "conjugate fields" of this local equilibrium ensemble are the coefficients of the corresponding conserved field operators {y(r, t)} ↔ −ν(r, t) + β(r, t)mw 2 (r, t)/2, β(r, t), −β(r, t)w(r, t)} .
It is interesting to note that this local equilibrium ensemble is also the "best choice" in the sense that it maximizes the information entropy for the given values of the conservative fields [32][33][34] . Furthermore, the local equilibrium state describes a "perfect fluid" in the sense that there is no dissipation, as discussed in the next section. With ρ ℓ N specified, the left sides of Eqs. (16a)-(16c) can be calculated as functionals of β(t), ν(t), w(t). Inverting these equations then gives these auxiliary fields as functionals of the average conserved fields of interest, and consequently To simplify and clarify the notation the fields n (r, t) , e 0 (r, t) , and u (r, t) will be denoted {ζ (r, t)} ≡ {n (r, t) , e 0 (r, t) , u (r, t)} so that (24) becomes All local equilibrium averages therefore become functionals of the conjugate variables, or equivalently of the average conserved fields. To make these two representation clear the notation in the following will be X (r) ; ρ ℓ [y(t)] ≡ X ℓ (r |y(t)) ≡ X (r |ζ(t)) , (27) where the one-to-one relationship of (16a) -(16c) can be expressed as y(r,t) = y(r |ζ(t)) or its inverse ζ(r,t) = ζ(r |y(t)). (28) This functional relationship is local in time. As a specific example the local equilibrium fluxes are s 0 (r) ; ρ ℓ [y(t)] = s ℓ 0 (r |y(t)) = s 0 (r |ζ(t)) (29) Since the operators s 0 (r) and t 0ij (r) are given by the derivation in the Supplemental Material, the functionals s 0 (r|ζ(t)) and t 0ij (r|ζ(t)) are given as explicit local equilibrium averages, and can be taken as formally known. In particular, it can be shown (see Appendix A) that t 0ij (r|ζ(t)) = π ij (r|n(t), e 0 (t)).
Here, π ij (r|n(t), e 0 (t)) is the equilibrium pressure tensor for a non-uniform system as a functional of the local density and internal energy density, independent of the flow velocity. With these results the macroscopic conservation laws become D t n(r, t) + n(r, t)∇ · u(r, t) = 0 , (33a) D t e 0 (r, t) + e 0 (r, t)∇ · u(r, t) + π ij (r|ζ(t))∂ i u j (r, t) +t * 0ij (r, t)∂ i u j (r, t) + ∇ · s * 0 (r, t) = −n(r, t)u(r, t) · ∇v ext (r, t) , (33b) The fluxes have been separated explicitly The contributions to the fluxes s * 0 (r, t) and t * 0ij (r, t) are those from ∆ N (t) in (15) and are discussed in Section IV.
Equations (33a)-(33c) are the basis for all analysis of continuum descriptions described here. They are still exact and general, with the constitutive equations for s * 0 and t * 0ij to be determined for specific states of interest. They remain to be specified from a solution to the Liouville-von Neumann equation. Explicit expressions are well-known for states with small space and time variations. In that case, the fluxes are obtained to first order in the spatial gradients of the fields, with transport coefficients given by Green-Kubo time correlation functions. For the fluid phase these are Fourier's law for the energy flux and Newton's viscosity law for the momentum flux. The resulting continuum equations are the Navier-Stokes hydrodynamic equations. The exact general form for arbitrary states is discussed in Section IV.

III. PERFECT FLUID HYDRODYNAMICS
In this section it is first shown that the contributions to dissipation (entropy production) are entirely due to the components of the fluxes s * 0 and t * 0ij . Hence, their neglect results in a "perfect fluid hydrodynamics." In the local density approximation (see below) they are the usual Euler-level hydrodynamics. More generally they extend the Euler equations to strong spatial non-locality.

A. Entropy production
The entropy and local entropy are defined here from the information entropy, maximized for the given average conserved fields [32][33][34] where is the index function for the local equilibrium ensemble in (19), so (37) becomes From Eq. (21) the operators for the conserved fields are denoted by ψ κ (r) (κ = 1, 2, 3) and the conjugate average fields are those of Eq. (23). They are defined in terms of the average conserved fields by the conditions eqs. (16a) to (16c). In this notation, Eq. (20) becomes The entropy production is defined by Next use (16a) -(16c), which define the conjugate fields in terms of the average conserved fields, to write Then using the macroscopic conservation laws, where The entropy production becomes where σ ext (t) is the entropy production due to the external force To interpret this further write γ iκ (r, t) = γ ℓ iκ (r|y(t)) + γ * iκ (r, t) and use the cyclic invariance of the trace to show that the contributions from γ ℓ iκ (r|y(t)) and the external force cancel (see Eq. (A17) of Appendix A) Consequently all of the entropy production is due to γ * iκ (r, t) All dissipation is associated with γ * iκ (r, t).

B. Hydrodynamics without dissipation
A special case of interest is conditions for which the irreversible energy and momentum fluxes can be neglected, s * 0i → 0 and t * 0ij → 0. This is never strictly true but can be a reasonable approximation for certain flows.
The resulting continuum equations are referred to as the perfect fluid equations. These equations have no unknown components (except the external forces to be chosen) and hence comprise a closed set of equations for the fields, justifying the terminology "hydrodynamics." The explicit functional dependence of π ij (r|ζ(t)) on the fields requires evaluation of the local equilibrium average, Eq. (32). The simplest approximation is the "local density approximation" whereby the functionals are replaced by functions of the fields at the point of interest r. The local equilibrium ensemble defined as a functional of the conjugate fields then becomes the equilibrium grand canonical ensemble, defined as a function of the fields at r, Similarly the relation of these conjugate fields to the average conserved fields, (21), becomes ψ κ (r); ρ e (β(r, t), ν(r, t), w(r, t)) = ψ κ (r, t).
The local equilibrium averages on the left side then become equilibrium averages. Then (52) expresses them as functions of the non-equilibrium conserved fields. The associated pressure tensor is The problem then reduces to a determination of the equilibrium pressure for the system being considered. This is a two step process. First calculate the pressure as a function of β(r, t), ν(r, t) and then use the equilibrium relations, β(r, t) = β e (n(r, t), e 0 (r, t)), ν(r, t) = ν e (n(r, t), e 0 (r, t)) to express the pressure as a function of n(r, t), e 0 (r, t). This is the usual procedure in the derivation of Navier-Stokes hydrodynamics for states with small spatial gradients and hence the local density approximation is justified. More generally, it is an uncontrolled assumption that must be removed for states with strong spatial inhomogeneities.

C. Pressure tensor
The pressure tensor and other local equilibrium averages simplify by transforming to a local rest frame to eliminate their dependence on the flow field. The local equilibrium ensemble depends on the velocity field only where G([q α ]) is the generation of the transformation. Therefore In particular the average momentum density can be written Comparison to the condition (16c) shows Let A 0 (r; u (r, t)) be an operator in the local rest frame; i.e. depends on the particle momenta through p ′ α = p α − mu (q α , t) . Then its local equilibrium average is independent of the velocity field In all of the following the averages of interest are for operators in the local rest frame and therefore the flow velocity is taken to be zero.
The pressure tensor is defined as the local equilibrium average of the microscopic momentum flux in the local rest frame t 0ij (r) and (61) In all of the following the choice is made.
Separate the pressure tensor into its diagonal and traceless parts The first term on the right will be referred to as the mechanical pressure The second term of (63) is the traceless part of the pressure tensor. Its volume integral must vanish since there is no external vector from which a nondiagonal tensor can be formed. Otherwise the offdiagonal contributions to π ℓ ij (r|β(t), ν(t)) are non-zero, in general.

Local pressures
The trace of the momentum flux operator follows from (60) where e K 0 (r) is the kinetic energy density operator The mechanical pressure is therefore where the notation introduced in Eq. (27) has been extended for the local equilibrium average of products by defining and use has been made of the identity To interpret the result in Eq. (68) it is useful to introduce a local pressure associated with the virial equation such that Eq. (72) is recognized as the virial equation associated with the local equilibrium ensemble. Then it follows that the mechanical pressure has the same volume integral as the virial pressure To obtain this use has been made of A third local pressure, the thermodynamic pressure π T (r|β(t), ν(t)), can be identified from the partition function Q [β(t), ν(t)] to define the thermodynamics of a local equilibrium system is extensive this can also be expressed as The volume derivative in (77) is evaluated in Appendix B with the result that the thermodynamic local pressure is the same as the local virial pressure In summary the volume integrals of all three local pressures are the same but the local mechanical pressure differs from the virial and thermodynamic local pressures In summary, there are several "reasonable" definitions for a local pressure that all agree with the thermodynamicglobal pressure but differ otherwise. In the present context the correct local pressure for the hydrodynamic equations is that obtained by direct evaluation of (68). A definition of the thermodynamic local pressure and its relationship to an associated pressure tensor coupling thermal and mechanical properties is given by Refs. 35 and 36.

Relation of pressure gradient to free energy gradient
It is often chosen to use a free energy density functional rather than the pressure to characterize the local equilibrium contribution to the momentum flux. To see how this is done, define the Legendre transform It follows from the definition of Q ℓ that the first functional derivatives are (82) Here F [β(t), n (t)] is the dimensionless free energy (in the uniform β limit F → βF ). Consider the special case The Finally, equating (84) and (85) and taking δr → 0 (86) This is the commonly used expression by others for the case of constant β(r,t). In our case it is important to remember there is another contribution to ∇π T (r |β(t), ν (t)) due to its variation with respect to β(r,t). This can be treated in the same way by introducing a double Legendre transformation from β(t), ν (t) to e 0 (t), n (t). The result is D. Relation to previous work As noted at the end of the Introduction above a new review of quantum hydrodynamics for plasmas has just appeared 10 . It provides an overview of much of the previous work on continuum descriptions for electrons, with extensive references. In most cases the phenomenology does not originate with traditional hydrodynamics, but is "reinvented". Consequently, the brief discussion here is motivated by recovering typical examples from the present exact analysis showing how some common limitations (weak inhomogeneities, weak coupling) can be removed in a controlled fashion.
Most previous work is based on the continuity equation and the momentum equation, without reference to the energy conservation law and without dissipation (i.e., special cases of the perfect fluid equations). Exceptions are references 8,9 and Ref. [37], although the latter is restricted to pure states. The implicit assumption in neglecting the energy equation is that the pressure tensor is independent of the internal energy. Then it is possible to show for the perfect fluid that the inverse temperature β is spatially uniform and dissipation is weak 3 . The perfect fluid equations (50a) -(50c) reduce to With the condition of uniform temperature these equations are still general, following from the underlying conservation laws for mass and momentum. In contrast, most earlier work is restricted to pure states or phenomenological macroscopic (coarse grained) Hamiltonian dynamics. Important restrictions to that work are explicit at the outset: 1) there is no dissipation, 2) the velocity field is irrotational (∇ × u(r, t) = 0), and 3) there are no off-diagonal elements to the pressure tensor. The origins and applications of continuum equations such as (88a) and (88b) have been reviewed extensively in references 11,38 in addition to the review noted above 10 . Instead, the focus is limited to results obtained from a phenomenological macroscopic Hamiltonian dynamics, due to Bloch 25 and developed primarily by others, e.g. Refs. 38-40. A Hamiltonian depending on the macroscopic fields n(r) and φ(r), where the flow velocity is defined in terms of this scalar potential by u(r) = ∇φ(r), is Here βF [n] is the dimensionless equilibrium free energy functional defined above, and v ext (r) is the given external single particle potential. Then the usual Hamilton's equations for the variables n(r) and φ(r) give the forms (88a) and (88b) with (Note that the free energy here includes the Hartree contribution, whereas that of reference 11 has extracted it explicitly). This result is in fact the same as the diagonal part of the pressure tensor in (63) where the thermodynamic pressure and free energies are related by (86) ∇π(r|n(t)) = n(r)∇ δF [n] δn(r) .
In most of the literature on generalized hydrodynamics for electrons (88b) therefore has the form There is considerable analysis of F [n(t)] within recent finite temperature DFT 6 . It is typically written as where F s , F H , F xc are the non-interacting, Hartree, and exchange-correlation contributions respectively. The parametrization of each term by the constant temperature β −1 has been left implicit. The most common case considered is weakly coupled electrons, for which F xc is neglected. Determination of the non-interacting contribution as an explicit functional of the density is still a formidable challenge, but for weakly inhomogeneous states it can be calculated from a gradient expansion. To second order in the density gradient it is 41 Here F TF [β, n] is the Thomas-Fermi free energy 42 . The coefficient of the square gradient contribution is known exactly for the special case of the uniform electron gas. More generally it is known in terms of the density response function for the corresponding uniform equilibrium system.
As a special case, the results for the uniform electron gas at zero temperature are where v H (r) is the Hartree potential, and where the Thomas-Fermi pressure π TF (β, n(r, t)) and Bohm potential v B (r) 23 are The appearance of the Bohm potential in this context has led to considerable confusion 11,26 . A result similar to (96) follows from an exact transformation of the Schroedinger equation (Madelung transform 9,22 ) except without the Thomas-Fermi pressure and without the factor of 1/9 for the Bohm potential. Its validity is strictly related to a pure state and does not extend to the mixed state ensembles considered here 38 . The approximation (96) is reasonable as long as the density of the underlying system stays close to homogeneity.

E. Application to strong coupling and strong inhomogeneities
The perfect fluid equations may have important applications to the complex electron -ion systems of warm, dense matter. Although neglecting dissipation, they describe the dominant convective dynamics without restric-tion to length and time scales. To incorporate the strong coupling all three contributions of the free energy functional in (93) must be included. Furthermore, since electrons in the vicinity of ions will have strong inhomogeneities, the limitation to gradient expansions must be relaxed. The range of temperatures for the electrons should extend from near zero to well above the Fermi temperature. This ambitious scope of applicability has been addressed successfully and two main advances are now available. An accurate representation of the equilibrium free energy for the uniform electron gas across the entire temperature/density plane has been developed from recent quantum Monte Carlo simulations 43,44 . This provides the essential local density approximation that is necessary to assure the uniform limit of any approximate functional for real electron-ion systems 45 . Functional development for strong inhomogeneities has included a generalized gradient approximation for both the noninteracting free energy 46 and the exchange-correlation free energy 27 . With these developments the perfect fluid hydrodynamics can be applied to warm dense matter conditions.

F. Linear modes
The simplest application of the perfect fluid equations is their linearization about an equilibrium reference state, n(r, t) = n (0) (r, t) + n (1) (r, t) + .., Choose u (0) (r, t) = 0 and the external force equal to zero. The solution sought is therefore the response to an initial perturbation. To lowest order the equations give uniform, time independent quantities n (0) (r, t) = n (0) , e 0 ). To next order they are 0 ) is the enthalpy density. The pressure tensor contribution is ∂ j π ij (r|n(t), e 0 (t)) = ∂ j dr ′ c n (|r − r ′ |)δ ij n (1) (r ′ , t) The pressure derivatives are evaluated at the uniform equilibrium state and are equilibrium response functions The linear equations (102a) -(102c) are most easily solved in a Fourier representation, with the notation The velocity field is written in terms of its longitudinal component ( k · u (1) ) and two transverse components ( e 1 · u (1) , e 2 · u (1) ). Transverse modes decouple from the others and have no time dependence. The remaining three variables have coupled dynamics with three eigenvalues λ(k) = (k (h c e0 (k) + n c n (k)), −k (h c e0 (k) + n c n (k)), 0). (106) The first two are the generalization of long wavelength sound modes to arbitrary length scales, while the last is a constant heat mode (since there is no dissipation). For Coulomb interactions they become generalizations of the plasmon modes. This is the simplest example of "generalized hydrodynamics." A more interesting example would be linearization around the stationary state imposed by an external force, such as that of a configuration of ions or a confining force. These will be discussed elsewhere. All of the perfect fluid hydrodynamics above neglects dissipation. Two attempts to improve this have been described. The first 11 makes a phenomenological modification of the coefficient a 2 in the gradient expansion of the free energy, Eq.(94), to make it non-local in space and time. This new dependence is then determined by the requirement that the linearized equations reproduce chosen representations for the linear response, e.g. random phase approximation. In the present context it is seen that this is attempting to capture effects that properly are contained in the dissipative fluxes, through an ad hoc modification of the local equilibrium free energy. A more consistent extension of the perfect fluid equations is to include an approximation to the irreversible fluxes, e.g. using those known in the Navier-Stokes limit for small space and time variations 7 . The result is a mixed representation of perfect fluid effects on all length scales but weak dissipation on long length scales. Still, it is a reasonable first attempt at including all relevant physical effects.
Several remarks are appropriate at this point.
1. These expressions for the irreversible fluxes are exact. Although complex, the difficult many body problem has been reduced to calculating local equilibrium time correlation functions. This is a formidable problem, extending that already present for linear response expressed in terms of corresponding equilibrium time correlation functions. The fields of the local equilibrium ensemble for the correlation functions are β (r,t) and ν (r,t) only, since the latter are in the local rest frame.
2. They are linear integral equations for the fluxes since they appear also on the right sides of (116) and (117) as well. This complication can be eliminated by introducing an appropriate projection operator, at the price of a more complex generator for the dynamics than the Liouville operator. This is described briefly in Appendix C.
3. These expressions provide the desired closure of the formal hydrodynamic equations (33a) -(33c), giving the irreversible fluxes as functionals of the local conserved fields. This is done indirectly, with (116) and (117) expressed first as functionals of the conjugate fields of the local equilibrium ensemble. These conjugate fields are directly related to the conserved fields through (16a) -(16c). The choice to represent the hydrodynamic equations in terms of the conserved fields or the conjugate fields is a matter of convenience for a given physical state.
4. The fluxes are non-local in space, just as the pressure tensor functional, and hence can describe states with strong spatial inhomogeneity. This is an extension beyond the usual Navier-Stokes hydrodynamics that is limited to small spatial gradients. 5. In contrast to the perfect fluid hydrodynamics, the irreversible fluxes are also non-local in time. Hence they describe all time scales, and also hysteresis (memory effects).
6. The response to small initial perturbations or external forces about equilibrium give the exact hydrodynamic response. Since the hydrodynamic equations describe all space and time scales, the exact response functions are obtained from this description. Of course, they are given in terms of the corresponding limits for the local equilibrium time correlation functions. This is the generalized linear hydrodynamics for equilibrium time correlation functions studied some time ago 13,14 .
7. The equilibrium response functions obtained from the hydrodynamic equations provide an interesting new exact representation. The generalized hydrodynamic form allows a direct cross-over to the small wavevector, low frequency limit -the Navier-Stokes limit. As discussed by Kadanoff and Martin 16 this is a singular limit not directly obtained by standard many body methods. 8. A linearization about local equilibrium for the pressure tensor and irreversible fluxes yields the non-linear Navier-Stokes hydrodynamics with the Green -Kubo expressions for transport coefficients. This is more general than linear response about equilibrium, in that the pressure and transport coefficients are local functions of the hydrodynamic fields at the space and time point of interest (in density functional theory this is the "local density approximation", extended to all fields).

V. SUMMARY AND DISCUSSION
A very general class of problems across many fields is currently being addressed via a continuum description of a few relevant macroscopic fields. For example, time dependent density functional theory focuses solely on the space and time dependent average electron number density. Properties of interest are presumed to be expressed as functionals of that density. A broader class of macroscopic fields, as considered here, include the energy density and momentum density which together with the number density have a special feature: they are averages of the set of exact locally conserved quantities. Historically, this property has provided the basic framework for phenomenological macroscopic balance equations for number, energy, and momentum. For special states, the phenomenology can be justified and supplemented by materials properties from the underlying statistical mechanics.
The present context is an interest in this approach for non-traditional systems (e.g., degenerate electrons) with the introduction of new phenomenology and applications whose validity and context is not yet fully explored. The objective here has been to suggest an alternative approach, that of starting from an exact continuum formulation and building more controlled approximations from it. For example, the Euler equations of hydrodynamics is a well-understood and useful approximation to the full Navier-Stokes equations; here, the perfect fluid hydrodynamics of Section III is understood as the analogous approximation (no dissipation) to the exact hydrodynamics, extended to arbitrary space and time scales. Furthermore, the perfect fluid equations are seen to provide the exact short time behavior of the system.
The analysis here is quite general within the limit of a single component non-relativistic fluid. Extensions to multi-component systems and inclusion of relativistic effects (e.g., graphene) follow directly. The primary results obtained are the exact balance equations (121) -(123) for the average local conserved fields, and the exact "constitutive equations" for the fluxes (116) and (117). The balance equations require a pressure tensor as a functional of the density and energy density in the local equilibrium state. This is an extension of the corresponding problem of density functional theory to include nonuniform energy density (or temperature) as well as number density. This formal extension has been developed but is still awaiting application 9,47 . The connection of the local hydrodynamic pressure to thermodynamics and an associated pressure tensor are discussed above and elsewhere 35,36,48 . Recent practical forms for strong coupling and strong spatial inhomogeneities have been developed recently 27,46 . Further developments for conditions of warm, dense matter remain a forefront computational challenge.
The constitutive equations for the dissipative energy and momentum fluxes are non-local functionals of the fields with respect to both space and time. While quite complex, their exact representation parallels closely that for Navier-Stokes hydrodynamics -linear in the spatial gradients of the conjugate fields with coefficients given by equilibrium time correlation functions that are local in space and time functions of the fields. These are Fourier's law and Newton's viscosity law. The many-body challenge has been compressed to calculating or modelling the time correlation functions. The exact forms here, (116) and (117), have a similar structure generalized to nonlocality and non-linearity in the conjugate fields, characterized by local equilibrium time correlation functions. The many-body challenge is now calculating or modelling these new correlation functions. The advantage of all the formal analysis behind these equations is to embed this limited (but difficult!) computation within a structure that assures the correct physics of the underlying conservation laws. Approximations made in the computation are expected to have quantitative rather than qualitative consequences.
There are several directions for first exploitation of these results. As noted above, an interesting application would be to the linear dynamics of a strongly coupled, inhomogeneous system of electrons in an external potential of a frozen configuration of neutralizing ions (a model for warm, dense matter) using the ideal fluid approximation (no dissipation). This requires calculating thermodynamic derivatives for the inhomogeneous system -time independent two point correlation functions. Such dynamics can be studied for comparison by existing methods of Born-Oppenheimer molecular dynamics using DFT generated forces 49 .
Another direction to explore is a local density approximation for spatially smooth states. This entail replacing all spatial functional dependence of the conserved and conjugate fields by their values at the same external field point. The pressure tensor and irreversible flux correlation functions then become corresponding equilibrium properties for which methods for evaluation are available. This can be done without a corresponding localization in time, extending the description to all time scales. One method available for such a calculation including strong coupling is the Kubo-Greenwood method within standard DFT 50,51 .
A new challenge presented by the results here is a better understanding of the local equilibrium state, both its structure and fluctuations. This is a direct extension of the corresponding equilibrium problem to non-uniform equilibrium states. There are many opportunities for application of systematic formal methods, models, and novel simulation techniques.
The functional derivative of a functional is defined from its first order differential variation The dependence of y on t has been suppressed to simplify notation at this point. Next, define so that U (x) obeys the equation Integrating from x = 0 to x = 1, and using the fact that U (0) = 1 gives and so to linear order in δy The tilde above an operator is defined in (115). The functional derivative of any local equilibrium average is therefore Now, restoring the dependence of y on t, the time derivative of ρ ℓ [y(t)] can be calculated directly Finally, consider the Liouville operator acting on the local equilibrium ensemble which obeys the equation Integrating from 0 to 1 gives or = dry α (r, t) ∂ i γ iα (r|y(t)) + γ ℓ iα (r|y(t)) e −η[y] − f iα (r|y(t)) + f where the microscopic conservation law corresponding to (43) - (45) has been used, and γ iα (r|y(t)) = 1 0 e −zη γ iα (r) − γ ℓ iα (r|y(t)) e zη dz (A16) One final simplification follows from the average of (A15), since X ℓ = 0 for any operator X. Then (A15) becomes

Local equilibrium averages
Consider a local equilibrium average of an operator A(r) in its local rest frame (function of the relative momenta p α − mu (q α )). Also note that ρ ℓ N also depends on the flow velocity only through the local momenta. Therefore, a unitary transformation (see (54), (55), and Supplemental Material) removes the velocity dependence; for simplicity of notation the dependence on y(t) is suppressed in this subsection.
Now, let T denote the anti-unitary time reversal operator (see Supplemental Material), with the property T p α T −1 = −p α , T q α T −1 = q α . Suppose A(r) has the property T A(r)T −1 = −A(r), then its local equilibrium average for u = 0 vanishes The first equality is non-trivial since the cyclic invariance of the trace does not hold for anti-unitary operators; it is proved in Section S3 of the Supplemental Material. Also use has been made of T ρ ℓ N T −1 = ρ ℓ N . As an example, consider the choice A(r) → s 0 (r), the energy flux in the local rest frame. According to (A20) its local equilibrium average must vanish, As a second example consider the matrix g αβ (r, r ′ ) = ψ α (r) ψ β (r ′ ) ℓ in the local rest frame. The conserved densities are odd or even under time reversal operation. Hence the matrix elements for densities with opposite signs must vanish.

Appendix B: Thermodynamic pressure
The objective here is to identify the local thermodynamic pressure given in (77) The volume derivative is taken at constant β, ν. To make explicit the volume dependence of Q [β(t), ν(t)] assume a cubic volume with V = L 3 and define the unitary operator U L with the properties 52 Recall that the flow velocity dependence has been transformed away. The explicit form for the exponent is The derivative with respect to L at constant β, ν is The desired volume derivative is, inverting back the scale transformation Finally, comparison with (77) gives (r| y(t)) Comparison with (71) shows this is the same as the local virial pressure π T (r|β(t), ν(t)) = π v (r|β(t), ν(t))

Appendix C: Solution to Liouville-von Neumann equation
The Liouville-von Neumann equation is Look for solutions of the form gives (C4) For simplicity of notation the case of a time independent external force has been chosen. More generally, the solution operator must be changed everywhere below as follows Using (A9) and (A18) this becomes and g −1 αβ (r, r ′ |y(t ′ )) is the inverse of g αβ (r ′ , r|y(t ′ )) given by (A8) dr ′′ g ασ (r, r ′′ |y(t ′ ))g −1 σβ (r ′′ , r ′ |y(t ′ )) = δ αβ δ(r − r ′ ) (C9) Next, using the macroscopic conservation law 44 to eliminate ∂ t ′ ψ α (r, t ′ ) leads to The contribution from the local equilibrium flux can be written .
The last line of (C14) is the component of the sources f α (r) orthogonal to the local conserved densities. The energy and momentum sources are both proportional to conserved densities and hence this term vanishes, leaving (C16)

Projection operator representation
It is useful to note some properties of the operators occurring in (C14). First, it is observed that the operators ψ β (r) and Ψ α (r|y) form a biorthogonal set in the sense ψ α (r)Ψ β (r ′ |y(t)) ℓ y(t) = dr ′′ ψ α (r) ψ σ (r ′′ |y(t)) ℓ y(t) Also Φ iα (r|y) is orthogonal to this set Define a related projection operator P t for trace operators Y Then it is straightforward to show and consequently, from (C14) The equation of motion for ∆ N from (C14) is Finally, integrating this with the initial condition ∆ N (0) = 0 gives the desired result The evolution operator U (t, t ′ ) gives the modified Liouville-von Neumann dynamics In contrast to (C14) there is now no longer an explicit dependence on γ * iα (r, t) in (C24).

S1. MICROSCOPIC CONSERVATION LAWS
The quantum mechanical microscopic conservation laws for the number, momentum, and energy density operators are given in Refs. S1 and S2. However, the details of the derivation are not given and the results do not include an external force. For completeness the general derivation is given here.
The time dependence of an operator A(t) which depends on the position and momenta of the system particles is given by where the Hamiltonian of interest is of the form Associated with the continuous symmetries are the usual conservation laws for particle number, linear momentum, energy, and angular momentum. Here only point particles are considered so the relevant conservation laws are those of particle number, momentum, and energy. The associated local densities are defined by where [A, B] + = AB + BA is the anticommutator. The symmetrized products ensure the Hermitian nature of these local density operators. In Eq. (S2) v ext (q α , t) is an external potential. As this external potential has not been included in the definition of e(r) in Eq. (S5), e(r) is referred to as the "intrinsic" energy density.
The following identities will be used below:

A. Number conservation
Local conservation laws follow exactly from the Hamiltonian dynamics. The simplest is the conservation of number density. In all derivations below the Greek letters, α, β, γ, · · · index the particles while Latin indices i, j, k, · · · denote Cartesian coordinates. Unless noted, Einstein summation is assumed for repeated indices.
Beginning with the time derivative of the number density, Using Eq. (S7), Inserting the definition of the momentum density in Eq. (S4) gives the final result,

B. Momentum conservation
The time derivative of p(r, t) is (S14) The two commutators appearing in the above equation will be treated separately. The former is The other commutator is where and Using Eqs. (S16) and (S17) in Eq. (S15) yields (S20) Define the kinetic part of the total momentum flux to be The force density arising from the external potential is defined by The second term on the right side requires some care. First, interchanging the dummy labels α, β and using Newton's third law gives N α =β=1 Taking half of the sum of these two equivalent expressions gives Consider a parameter λ defining an arbitrary path from x(λ β ) = q β to x(λ α ) = q α . The following identity then holds, Thus the interaction force term becomes where the potential part of the momentum flux density is Using Eq. (S29) in Eq. (S24) gives Defining the total momentum flux density as the sum of its kinetic and potential parts, t jk = t K jk + t P jk , gives a compact form for the conservation of momentum equation, (S32)

C. Energy Conservation
The conservation of energy equation also follows from the time derivative of the local energy density, using Eqs.
The following two sections deal separately with each commutator appearing in Eq. (S34). The first commutator can be written as the sum of three terms, N α=1 1 2

Combining results
Using Eqs. (S35) and (S38) in Eq. (S34), The terms in the third and fourth lines can be combined, The third term can further be rewritten as where use has been made of Eq. (S28) in the next to last line and dummy indices r, k have been swapped in the last line. The energy equation is therefore Moving the three derivative terms to the left side, this can be rewritten as ∂ ∂t e(r, t) + ∇ · s(r, t) = w(r, t) with the energy flux and the work done by the external potential This work term can be simplified as follows, This yields the final form for the work term as

S2. UNITARY TRANSFORMATION TO REST FRAME MOMENTA
Operators and phase functions of interest include those in the local reference frame. For simplicity consideration here is restricted to single particle functions, where u (q) is the average flow velocity at the particle position q. In the classical case a change of variables in the momentum average to gives the simpler form which is independent of the velocity field. The objective here is to find a unitary operator that does the same for quantum mechanical operators, The generator of the classical canonical transformation in Eq. (S49) obeys S3 Using the identity generator q · p ′ a new generator for the deviation is introduced In the following G(q, p ′ ) will be referred to as the generator rather than F (q, p ′ ). Note that for this generator Now consider the associated infinitesimal transformation, mu i (q) → ǫmu i (q), G(q, p ′ ) → ǫG(q, p ′ ). To first order in ǫ, where {, } denotes the Poisson bracket. The Dirac -Weyl quantization of this procedure is given by although some limitations exist. According to this quantization procedure, the classical infinitesimal transformation corresponds to the finite transformation To demonstrate that this transformation does in fact give the desired result, transform a quantity A(q ′ , r ′ ), Since G = G (q) it commutes with q and so Next define The following identity was used above: Integrating (S63) from 0 to 1 gives Therefore (S60) gives the desired result In the text, most final results are expressed as local equilibrium correlation functions in the rest frame, e.g.
Recall that the subscript zero denotes the chosen operator with all momenta replaced by their local rest frame momenta p α → p α − mu (q α ). The local equilibrium ensemble also has this form. Then by cyclic invariance of the trace and such correlation functions are independent of the local flow velocity.

S3. TIME REVERSAL SYMMETRY
The objective of this section is to prove a property of local equilibrium correlation functions for operators that have a definite sign with respect to changes in the sign of the particle momentum operator. In the classical case this is accomplished by introduction of the parity operator for the momentum variables only. In the quantum case the parity operator changes both the momentum and position operators and the analysis no longer holds. Instead, an operator that changes the momentum sign but not the position operator sign is considered.
As in the classical case the quantum time-reversal operator T generates a reversal of motion operator in the sense It follows (e.g., coordinate representation for p) that T i = −iT , and hence T is an anti-unitary operator where c 1 and c 2 are arbitrary complex numbers. In addition it is required that T be norm-preserving where |α and |β are members of a complete basis set {α}, and It then follows that T is anti-unitary S4 Consider again the correlation function of (S68) It is shown here that this has the property Thus if both A 0 (r) and B 0 (r ′ ) have a definite sign under time reversal, the correlation function vanishes when the signs are different. The proof of (S76) is complicated by the anti-unitary nature of T which does not obey the cyclic invariance of the trace. The proof is as follows. Define where A is an arbitrary linear operator. Then where the second equality follows from the anti-unitarity of T , (S74). Writing γ|α and α| γ explictly in terms of A As a special case, the diagonal elements for a self-adjoint operator are Now the correlation function in (S75) can be expressed as The last line follows from (S80). Since |α and | α are in one-to-one correspondence the latter is a complete basis set. Furthermore T ρ ℓ T −1 = ρ ℓ and (S83) becomes (S76).

S4. TRANSFORMATION TO LOCAL REST FRAME
Operator functions of the particle positions and momenta typically refer to a fixed laboratory frame. They define both convective motion and motion relative to each particle's local average velocity field. Denote such an operator by A({q, p}), where {q, p} denotes the N particle positions and momenta. It is convenient to extract from this the corresponding operator representing only motion relative to the average flow A 0 ({q, p}) ≡ A({q, p − mu (q)}). They are related by the generator of eq.(S67), Operators of primary interest are the conserved densities and associated fluxes, {ψ(r)} = {n(r), e(r), p(r)}, (S85) The corresponding rest frame operators are easily found to be and the rest frame fluxes are A concise form for these transformations to the local rest frame can be described in matrix representation, and the flux vector transforms to the rest frame as where the matrix A(u) is given by The second equality is a simplified notation, whereby u is a three-component column vector, u T is its transpose, 0 T is the transposed zero vector, and I is the 3 × 3 identity matrix. To obtain (S90) use has been made of the fact that It is clear that the results (S89) and (S90) apply as well for the average conserved densities and fluxes ψ 0α (r) = A αβ (u)ψ β (r) (S94) and the flux vector transforms to the rest frame as since the average is a linear operation and commutes with A.