Molecular hydrodynamic theory of the velocity autocorrelation function

The velocity autocorrelation function (VACF) encapsulates extensive information about a fluid's molecular-structural and hydrodynamic properties. We address the following fundamental question: How well can a purely hydrodynamic description recover the molecular features of a fluid as exhibited by the VACF? To this end, we formulate a bona fide hydrodynamic theory of the tagged-particle VACF for simple fluids. Our approach is distinguished from previous efforts in two key ways: collective hydrodynamic modes are modeled by \emph{linear} hydrodynamic equations; the fluid's static kinetic energy spectrum is identified as a necessary initial condition for the momentum current correlation. Our formulation leads to a natural physical interpretation of the hydrodynamic VACF as a superposition of quasinormal hydrodynamic modes weighted commensurately with the static kinetic energy spectrum, which appears to be essential to bridging continuum hydrodynamical behavior and discrete-particle kinetics. Our methodology yields VACF calculations quantitatively on par with existing approaches for liquid noble gases and alkali metals; moreover, our hydrodynamic model for the self-intermediate scattering function extends the applicable domain to low densities where the Schmidt number is of order unity, enabling calculations for gases and supercritical fluids.


INTRODUCTION
][9][10][11] Using simple physical arguments, they recognized that this protracted decay was due to delayed viscous momentum transport-a phenomenon known as hydrodynamic memory.
In this Communication, we formulate a theory of the VACF from a purely hydrodynamic standpoint that directly confronts these questions.Central to our framework is a general set of physical desiderata (Table I), which constrains a chosen hydrodynamic model so as to correctly reproduce both short-time molecular kinetics and collective dynamics across all timescales while minimizing ad hoc assumptions.Crucially, we find physical consistency requires that the equal-time velocity covariance be an integrable distribution over wavenumber, which we identify as an initial condition corresponding to the fluid's static kinetic energy spectrum.8][59] The resulting moment-equations naturally yield wavenumber-dependent current correlations that reproduce, among other phenomena, molecular-scale  viscoelasticity. 60,61By contrast, generalized hydrodynamics 1,3,62 begins with the Navier-Stokes equations, extending them to molecular scales by phenomenologically promoting transport coefficients to nonlocal quantities with frequency-and wavenumber-dependence.][64][65][66][67][68][69] The present theory offers important advantages, however.In particular, it implies a natural physical interpretation of the hydrodynamic VACF: a superposition of quasinormal hydrodynamic modes, each mode k having a wavelength 2π/k and weight u 2 k proportional to the equilibrium probability density of fluid kinetic energy.Moreover, unlike other approaches, the methodology enables realistic VACF calculations not only for simple liquids, but also gases using the same underlying framework.We apply the methodology to representative fluids-liquid rubidium and argon, and gaseous and supercritical argon-and find remarkable agreement with MD calculations. 28,31,34,70EORY OF THE HYDRODYNAMIC VACF Let u(x, t) be the Eulerian velocity field of a fluid derived as the velocity moment of an exact single-particle distribution function, f (x, v, t), and x α (t) be the position of constituent particle α of mass m.We then require the fluid velocity at x α (t) to be equal to the particle velocity: u[x α (t), t] = v α (t), where v α (t) ≡ ẋα (t).Expressing the Lagrangian velocity using Fourier components u k (t), the velocity covariance where in the last line we assume a spatially homogeneous system and decompose the ensemble average into a product of averages.Importantly, we assume a priori uncorrelated particle displacements and fluid velocities whose underlying correlations will be reproduced by applying appropriate physical constraints (cf.Table I).
Assuming an isotropic fluid, we decompose the velocity covariance, is the selfintermediate scattering function (SISF), the characteristic function for tagged-particle displacements.Altogether, this yields a general expression for the VACF where N ≡ ∞ 0 dkk 2 q(k) −1 is the normalization at t = 0.

Equilibrium energy partitioning
From Eq. (3), it is apparent that the customary statement of equipartition, q(k) = k B T /m, leads to a divergent integral.This is unsurprising, as a true white noise spectrum is unphysical.Thus, the fluid static kinetic energy spectrum, or spectral density for short, must contain a microscopic parameter, a, that attenuates large-wavenumber contributions.][73] The inverse Fourier transform of q(k) is formally identical to the "form factor," f (r), in the velocity-field method, wherein a is a priori treated as an effective molecular radius, which, in practice, is fixed near the mean intermolecular distance (V /N ) 1/3 = n −1/3 0 .However, by only requiring consistency with the desiderata in Table I-viz.reqs.( 2) and (5)-we deduce that a more generally represents a kinetic correlation length, a scale above which a continuum description applies.Indeed, we find the shape of q(k), which must depend on the detailed intermolecular potential, strongly influences VACF oscillations, starkly contrasting with the form factor's largely peripheral role in ensuring normalizability.Note that Eq. (4) implies q(k) ≥ 0, whereas f (k) is negative for a range of wavenumbers; conversely, q(r) = ∞ 0 dk k 2 q(k) sin(kr)/kr can oscillate radially about zero, thereby inducing oscillations in the VACF.

The hydrodynamic VACF
For an isotropic fluid, the VACF decomposes into longitudinal and transverse components with respect to wavevector k as ψ(t) = 1 3 ψ (t) + 2 3 ψ ⊥ (t); Eq. (3) becomes where, after changing notation from q(k) to q(ka) and setting N = N p a 3 , we have explicitly introduced a and a parameter p that controls the shape of q(ka).
Equation ( 5) and the desiderata in Table I constitute a general formulation of the hydrodynamic VACF.To implement the framework, one specifies the static spectrum, q(ka), and correlation functions, C ,⊥ (k, t) and F s (k, t).In what follows, we develop hydrodynamic models satisfying reqs.(1-7) in Table I to obtain solutions for C ,⊥ (k, t) and F s (k, t).We then introduce representative models for q(ka) that yield simple-fluid VACFs in good agreement with MD calculations.

HYDRODYNAMIC MODEL FORMULATION
To describe collective hydrodynamic modes, we develop a regularized variant of the well-known 10-moment equations, a relaxation system of partial differential equations (PDEs) that we call R10. 59When appropriately constrained by the desiderata, R10 reproduces accurate VACFs for simple liquids, as well as realistic VACFs for intermediate-density fluids and gases.Note that R10 is not necessarily the optimal (or unique) hydrodynamic model consistent with the desiderata; in principle, our framework is compatible with other methodologies, including generalized hydrodynamics 1,62,74 and the GENERIC formalism. 50,5110 may be obtained as a generalization of the BGK-Boltzmann equation using the moment method of hydrodynamics, where we consider two collisional relaxation rates, one each for the longitudinal and transverse components of the deviatoric stress tensor.We derive analytic expressions for longitudinal and transverse current correlations, which capture damped molecular-scale sound and elastic shear waves, respectively.In particular, the corresponding memory equations and kernels [req.(1)] are not assumed, as is commonly done, but rather implied by the R10 PDEs.
To describe the self -motion of tagged particle α, we propose a hydrodynamic model for the self-density, n s , motivated by the R10 equations for collective hydrodynamics variables (n, j, and S).We then derive analytic expressions for the SISF via F s (k, t) = n * s (k, 0)n s (k, t) .As there are no known exact solutions for F s , approximations have often leveraged the Gaussian assumption, which is exact in the small-and large-k limits and, conveniently, directly relates to the mean-square displacement (MSD), e.g., is the MSD of tagged particle α.For instance, a cumulant expansion of F s (k, t) with a Gaussian leading term directly reveals non-Gaussian effects, which are known to be relatively small for simple liquids. 1,75Unlike common approximations, our model for F s satisfies reqs.(1-6), which, along with our model for C ⊥ , extends the description to low Schmidt number, Sc ≡ ν/D s ∼ O(1), and reproduces the exponential decay range of gases-a feature not captured by simple diffusion (i.e., Fick's law) and other Gaussian models.

Regularized 10-moment model (R10)
Consider the following linear transport equations for an adiabatic equation of state: where ρ(x, t) is the mass density, j(x, t) is the mass current density, S(x, t) is the deviatoric stress tensor, and ε is the rate-of-strain tensor; v 0 ≡ k B T /m and v are thermal and longitudinal phase velocities, respectively.We define the Fourier transform of the collision frequency tensor Υ(x) as where k ≡ k/k, 1 is the unit tensor, and kk and 1 − kk are longitudinal (normal) and transverse (shear) projection operators; γ and υ are the respective collision frequencies.Similarly, for the stress diffusion tensor, D(x), with regularization diffusion coefficients D and D ⊥ .
Two points should be highlighted.First, the Fouriertransformed collision term, −Υ k • S k , is a relaxation approximation related to the BGK collision operator in the BGK-Boltzmann equation that gives rise to viscoelasticity: e.g., υ −1 is the transverse component's Maxwell relaxation time.Second, the Fourier-transformed diffusive regularization term, −k 2 D k • S k , extends the description to higher Knudsen (Kn) number 50,55 -essential for molecular-scale fluid flows where Kn ∼ O(1) and conventional Navier-Stokes fails. 76,77Importantly, regularization implements moment closure when constrained by the Green-Kubo relation for self-diffusion [req.(7)], as well as eliminating spatial structure below physically meaningful scales [req.(4)].

Relaxation limit of the stress components
Working in Fourier space, the transverse projection of Eqs. ( 7) and ( 8) yields while the longitudinal projection of Eqs. ( 6) to ( 8) yields We then make the key assumption that the longitudinal stress relaxation rate, γ + k 2 D , is faster than any other timescale, which circumvents solving a third-order differential equation for C and proves to be a good approximation; Eq. ( 15) relaxes to where Here, γ is determined in the small-k limit from available data for the bulk viscosity, µ ≡ v 2 0 /γ, and D is treated as a free parameter; the intuitive choice D → 4µ/3 yields good results for present calculations.

Modal current correlation functions
Using Eqs. ( 11) and ( 12) for the transverse component, and the relaxation approximation for the longitudinal component, Eqs. ( 13), ( 14) and ( 16), we derive ordinary differential equations (ODEs) in time for the modal correlation functions using j k (t) = ρ 0 u k (t), ρ 0 is the equilibrium mass density.For transverse current correlations, C ⊥ , we combine Eqs. ( 11) and ( 12) to obtain a second-order ODE for the transverse current where 18) has a memory equation form, req.(1), with kernel Similarly, substitution of the relaxation limit, Eq. ( 16), into Eq.( 14) eventually yields a secondorder ODE for the longitudinal current correlations where With the initial conditions C ,⊥ (k, 0) = 1 and Ċ ,⊥ (k, 0) = 0 [reqs.(2-3)], solutions to Eqs. ( 18) and ( 19) are readily found: where Equa- tions (20) and (21) have the flexibility to satisfy all requirements in Table I; the corresponding dynamic spectral densities are provided in Appendix A.

Self-intermediate scattering function
The regularized relaxation models represented by Eqs.(11) and (12) and Eqs. ( 13) to (15) for the collective hydrodynamic variables suggest the following model where n s (k, t) is the (self-)density of a non-interacting collection of tagged particles (i.e., test particles 75 ), j s (k, t) the self-current, γ s the Brownian collision frequency corresponding to (Stokes) friction, and D s the self-diffusion coefficient; here, we take D s → v 2 0 /γ s , which enforces req.( 4) and ensures the positivity of F s for all k and t.
Equations ( 22) and ( 23) are a natural generalization of Fick's Law of self-diffusion, 44,45 which is recovered in the full relaxation limit (γ s → ∞, k → 0) where Eqs.(22)  and ( 23) reduce to a diffusion equation for n s (k, t).The full ODE corresponding to Eqs. ( 22) and ( 23) is 24) the roots of which factor nicely to give Note that density fluctuations decay exponentially for large k.To see that this model is reasonable, consider Eq. ( 24) in memory equation form where The Markovian solution, obtained by freezing F s (k, t ′ ) at the upper limit t, is which, when D s = 0, yields the formal result for conventional Langevin dynamics. 1Also note that when D s → 0 in Eq. ( 25), F s (k, t) → 1 and Eq. ( 27) would not satisfy req.(4).It is unsurprising that D s > 0 is necessary for a physically meaningful SISF.

IMPLEMENTATION OF THE FRAMEWORK
The VACF is calculated by evaluating Eq. ( 5) using the transverse and longitudinal current correlations, Eq. ( 20) and Eq. ( 21), and SISF, Eq. ( 25), along with the static spectrum, q(ka), discussed subsequently; importantly, all VACF calculations use these equations (and parameters).However, determining gas parameters (Appendix D) is considerably more complicated since the short-time sum rule [req.(5)] only applies to dense fluids, F s ≈ 1 cannot be assumed for req.(7), and empirical data for transport parameters is limited.A suitable kinetic theory, such as Enskog theory (used here) or modifications thereof, [78][79][80] can be used with our framework to derive sum rules connecting molecular parameters to macroscopic quantities, 10,81 as well as estimating transport coefficients in lieu of empirical inputs.

Static spectrum
Pending a first-principles derivation of q(ka), we consider representative two-parameter models: a symmetric generalized Gaussian with p ≥ 1, and generalized Lorentzian with p ≥ 6, where the normalizations depend on shape parameter p via N G ≡ N G p a 3 and N L ≡ N L p a 3 .Equations ( 28) and ( 29) approach the equipartition spectrum for small a; softer intermolecular potentials correspond to sharper spectral cutoffs (larger p), which amplify oscillations in the VACF [cf.3][84][85] Note that normalization [req.(2)] yields N G p = p/[8 1/p Γ(3/p)], whereas solutions for N L p are unavailable for general p, so Eq. ( 28) is preferred for analytic manipulation.
But the short-time sum rule [req.(5)] generally yields a ≤ a 0 (cf.Appendix B), implying that the molecularscale parameters υ = v 2 0 /ν and γ s = v 2 0 /D s .Assigning consistent values for a, υ, and γ s therefore requires a more careful procedure to concomitantly satisfy reqs.( 5) and (6).Importantly, treating the products a 2 υ and a 2 γ s in Eq. ( 31) as invariant quantities preserves the asymptotic form.

Parameter determination at lower densities
Given the paucity of numerical and experimental data for dilute monatomic fluids, we determine inputs using Enskog theory supplemented by heat capacity data from the National Institute of Standards and Technology (NIST) Chemistry WebBook. 88For brevity, procedural details are provided in Appendix D. II.Calculations for argon-like gaseous and supercritical fluids use expressions for the Enskog viscosity and diffusion coefficients combined with reqs.( 4) and (7) to derive reasonable values for the key parameters D ⊥ , υ, γ s , and a (cf.Appendix D).All calculations in presented in the article and its appendices can be performed with the Mathematica notebook provided in the supplementary material.

VACF calculations for liquid rubidium and argon use values listed in Table
Figure 1 shows the VACF for liquid rubidium using a generalized Lorentzian spectrum with p = 14, which should be compared to the velocity-field results (Ref.46, Fig. 1; Ref. 70, Fig. 2).The oscillatory VACF behavior of liquid alkali metals observed in MD simulations 46,47,70,96 requires a relatively sharp cutoff (large p), which corresponds to the relatively soft repulsive core of the intermolecular potential as compared to liquid argon. 21,82,83igure 2 shows the VACF for liquid argon near the triple point using a Gaussian (p = 2) spectrum (cf.velocityfield result: Ref. 47, Fig. 2).We remark that the plateau often seen in liquid argon(-like) VACFs from MD 25-28,92,97-99 can be reasonably reproduced in our formulation using, e.g., a generalized Lorentzian spectrum with p ≈ 7.5, v 2 ≈ 3.0v 2 0 , and adjusting D s upward by ∼10% so as to increase oscillations primarily in the longitudinal (but not transverse) component.
Figure 3 shows VACFs for argon-like gaseous and supercritical fluids.The exponential and diffusive (t −3/2 ) decay ranges are a direct consequence of the explicit handoff in Eq. ( 25) and Eq.(D3).These prominent features were observed in MD calculations for hardspheres 31 and a supercritical Lennard-Jones (LJ) fluid, 34 as well as immersed-particle fluctuating hydrodynamics simulations at low Sc 43,100,101 and analytic calculations for Basset-Boussinesq-Oseen (BBO) dynamics 12,102-105 with general slip boundary conditions. 106In particular, the dips at early times (Fig. 3, t 10 ps) come from the longitudinal current, which appears to capture the effect of strongly damped sound waves.Similar features are clearly evident in MD calculations: see Ref. 31, Fig.S1 and also Appendix E for an indirect comparison with Ref. 34, Fig. 1, which shows excellent quantitative agreement.Analytic calculations for "BBO particles" 107 also exhibit qualitative similarities (Ref.106, Figs. 2 and 3).

DISCUSSION
The present theory reflects an eclectic synthesis of ideas dispersed throughout the literature along with several new concepts.We have organized the theoretical formulation and framework in this Communication so as to highlight important results.In summary, we: parameters that recover the short-time VACF behavior while preserving the correct long-time t −3/2 decay.It is worth mentioning that one can derive telegrapher's equations from R10 (for D ⊥ or D s = 0).Trachenko 108 derived a telegrapher's equation as a continuum liquid dynamics model.0][111] As pointed out by Khrapak, 112 Zwanzig's speculative model of molecular self-motion in a liquid 113 -where collective rearrangements correspond to configurational transitions between metastable equilibria-is consistent with the persistent random walk picture.From our hydrodynamic standpoint, this persistence originates from the finite relaxation time of molecular-scale stresses.
Regularization, however, is essential-especially at low densities.Importantly, Eq. (25) shows that exponential decay occurs when k 2 D s > γ s , which implies kλ E > 1, where λ E ∼ v 0 /γ s is the Enskog mean free path.More precisely, when a < k −1 < λ E , the expected exponential decay range of a gas emerges from the contributions of large-k modes.Also, note that the slowest relaxation rate of C ⊥ (k, t) is given by the smaller of the two exponents in Eq. ( 20), which, when expanded for large k, gives −v 2 0 /D ⊥ when D ⊥ > 0 (cf.Appendix D).Thus, the regularization coefficients D ⊥ and D s , the latter of which arises from our hydrodynamic model for the self-density, are necessary for obtaining the expected exponential decay of the dilute gas VACF.
The exponential decay of the SISF (and VACF) at low densities also hints at a deeper connection to VACFs for large Brownian particles, which also exhibit exponential decay. 114][119][120][121][122][123][124][125] We remark that the non-Gaussian behavior of the VACF at intermediate times, which switches from damped oscillations to pure exponential decay, coincides with the emergent exponential ranges in F s and C ⊥ .It is therefore possible that the hydrodynamic VACF formulation can lend dynamical insight into the liquid-vapor phase transition not otherwise available through other methods.Realizing a fully capable VACF theory would, however, require a means (i.e., a sum rule 10,81 ) to determine a as a continuous function of density (and temperature) through the phase transition without relying solely on the short-time sum rule involving ω E [req.(5)], which is ill-defined at lower densities.][131][132][133][134][135][136] Extensions to the present work include deriving q(ka) from first principles, as well as (numerically) solving the full zero-frequency constraint for D ⊥ [req.(7)] at all densities and third-order ODE for C (k, t) (i.e., beyond the relaxation limit).However, we expect the present theory to be valuable well beyond VACF calculations.For example, Alder and Wainwright 6 and, more recently, Han et al. 30 and Lesnicki and Vuilleumier, 29 have convincingly demonstrated that collective motions in discreteparticle fluids are well represented by hydrodynamic flow fields down to the single-particle scale.This suggests that stochastically driven R10 equations, 77 which gen-eralize the Landau-Lifschitz Navier-Stokes equations, 137 would represent a viable molecular-hydrodynamic model capable of reproducing molecular-scale flows.
Indeed, the efficacy of our VACF formulation rests largely on q(ka).Even apparently subtle differences in the shape of q(ka) substantially alter the character of the VACF, e.g., the oscillatory nature for liquid rubidium or "plateau" for liquid argon (data not shown).The details of the VACF depend on the manner in which the fluid's distribution of kinetic energy transitions from (approximate) equipartition at long wavelengths to zero below the molecular scale.The kinetic energy distribution must, in turn, depend on the intermolecular potential, which also determines the radial distribution function, g(r), or static structure factor, S(k).However, whereas g(r) or S(k) characterizes static density correlations (zeroth velocitymoment), q(k) characterizes static momentum correlations (first velocity-moment).Given its pivotal role in capturing molecular-scale behavior, it is thus our belief that the static spectrum is key to describing how continuum hydrodynamic modes emerge from the molecular scale and bridging the continuum and discrete-particle perspectives.
where τ c = ω −1 E is the timescale characterizing the initial (de)correlation of the VACF and the Einstein frequency, ω E , is formally defined via where g(r) is the radial distribution function and φ(r) is the intermolecular potential.For large Sc-assumed to be the case for liquids-the short-time expansion of the modal correlation function dominates that from the selfintermediate scattering function.That is, when ν ≫ D s , we have v 2 0 /υ ≫ v 2 0 /γ s and then γ s ≫ υ.Thus, in the liquid state, F s (k, t) ≈ 1 over the timescale τ c and the initial decay of the VACF (and the location of the first zero-crossing) is controlled by the modal current correlation functions, C ⊥ and C .The integrand of the hydrodynamic VACF formula, Eq. ( 5), can thus be approximated as and so the short-time VACF is approximately B4) Equating the quadratic terms in Eq. (B1) and Eq.(B4), we obtain which is the second frequency-moment condition (i.e., fsum rule). 3o obtain Eq. ( 33), explicitly evaluate Eq. (B5) for the generalized Gaussian spectrum The scale a can now be deduced from the correlation time τ c through the second wavenumber moment of the static spectrum, q(k): where v 2 a ≡ v 2 + 2v 2 0 /3.Carrying out the integral analytically and rearranging to solve for a yields a = which agrees with Eq. (33).For the specific case of a pure Gaussian static spectrum, a = √ 3v a /ω E .In the liquid state, it is seen that the main decay timescale is controlled by the molecular scale a via Appendix C: Zero-frequency constraint for liquids [req.(7)] The Green-Kubo relation for self-diffusion is In the liquid regime, it is reasonable to take C (k, t) → 0 and F s (k, t) → 1, which allows us to evaluate Eq. (C1) analytically using Eq. ( 5) to obtain where χ(p) is a coefficient that depends on the sharpness of the spectral cutoff.After rearranging, we obtain which, as shown below, yields analytic constraints for D ⊥ when D s is known (e.g., experimentally) and the static spectrum is represented by either the generalized Gaussian or Lorentzian forms [Eqs.( 28) and ( 29)].
For the transverse current, Eq. (C1) becomes Carrying out the time integration first yields so that the k-integral that remains to be evaluated is For a generalized Gaussian spectrum, Eq. (C6) becomes and, for a generalized Lorentzian, Inspection of Eqs.(C7) and (C8) reveals that the transverse current correlation uniquely contributes a pdependent coefficient of the a 2 υ term, χ(p), unlike the longitudinal current correlation, as well as depending on the functional form of q(ka) = q p (ka): Thus, D ⊥ can be readily computed using Eq.(C7) or Eq.(C8) with Eq. (C3).arises only from the product F s (k, t)C ⊥ (k, t).Then, one might define an effective decay rate, Ω 0 , to be the sum of the SISF and transverse current collision rates, i.e., Ω 0 = γ s + υ.It will be shown below that this assumption is not unreasonable and can be formally justified by enforcing req.( 7) using large-k approximations for the correlation functions.Specifically, the exponential range should arise for values of k where kv 0 > τ and Ω 0 = τ −1 [req.(4)].Indeed, unlike in the liquid case, the SISF cannot be treated as approximately constant (i.e., F s = 1) over the time interval with the dominant contribution to self-diffusion, which precludes an exact analytic integration over k.However, for dilute monatomic gases, the dominant contribution to D s generally comes from the large-k forms of C ⊥ (k, t) and F s (k, t), in the Green-Kubo relation allowing an approximate k integration.In this scenario, regularization plays a critical role, as seen by considering the slowest relaxation rate of C ⊥ (k, t)-the smaller of the two exponents in Eq. ( 20), which, for large k, is This relaxation rate must be finite as k → ∞ [req.( 4)], which, evidently, necessitates the inclusion of the regularization diffusion coefficient, D ⊥ .A similar argument holds for the exponential range of F s (k, t), where D s is the corresponding regularization coefficient.
Similarly, the longitudinal relaxation rate for largek is ∼v 2 /µ; however, the bulk viscosity of a dilute monatomic gas is typically very small compared to its shear viscosity, [138][139][140] making its contribution to the diffusivity small as well.To see this analytically, consider the VACF with a Gaussian static spectrum (p = 2) and bulk viscosity set to zero: the longitudinal contribution is the rapidly decaying form ψ (t) ∼ (1 − v 2 t 2 /a 2 ) exp −v 2 t 2 /2a 2 , which integrates exactly to zero.Thus, while C can affect the short-time VACF structure, only C ⊥ and F s dictate the exponential range and overall diffusivity under dilute conditions.
Zero-frequency constraint [req.(7)].Excluding the longitudinal component, the Green-Kubo relation for large k becomes [req.(7)] Note the appearance of the additional parameter γ s due to the treatment of finite Sc.To uniquely determine the parameters, we now require a separate relation between D ⊥ and υ or another way to determine γ s .At present, it seems reasonable to assume D ⊥ → v 2 0 /υ, which is consistent with the assumption that D s → v 2 0 /γ s , which was used in Eqs. ( 22) to ( 25)-the equations for the selfdensity and SISF.We remark that when D ⊥ = v 2 0 /υ, the roots of Eq. ( 18) factor nicely to give matching the neat form of Eq. ( 25) and preserving the positivity of the VACF in accordance with what is expected for a dilute gas.
Setting D ⊥ = v 2 0 /υ in Eq. (D2), we now have Clearly, we cannot directly set D s = v 2 0 /γ s on the lefthand side of Eq. (D4), as there would be no positive solution for υ.We instead assume, as in the liquid case, that γ s represents the rescaled (self-)collision frequency, while γ s0 is its corresponding base value determined by D s ≡ v 2 0 /γ s0 .With this assumption, Eq. (D4) becomes Parameter rescaling [req.(6)].Finally, we require that the Schmidt number be invariant after rescaling: Sc ≡ ν/D s = γ s0 /υ 0 = γ s /υ, which fixes the ratio of the rescaled collision frequencies to the ratio of their base values.Substituting υ = γ s /Sc into Eq.(D5) and solving for γ s leads to the condition where the second equality is obtained by recalling that a is determined by preserving the long-time asymptotic form of the VACF [req.(6)], which implies Thus, γ s , υ, and γ may be determined, respectively, via γ s0 = v 2 0 /D s , υ 0 = v 2 0 /ν, and γ 0 = v 2 0 /µ given known inputs for the transport coefficients (D s , ν, and µ).
Parameters for gaseous and supercritical argon.Given the sparsity of numerical and experimental data for transport coefficients of dilute monatomic gases, the VACFs in Fig. 3 were produced using transport coefficients calculated from Enskog theory, 11,141,142 supplemented by isothermal data for argon (at 130 K and 300 K) from the NIST Chemistry WebBook.For the transport coefficients, we took 0 /υ 0 , and µ → ν bE = v 2 0 /γ 0 , where are the Enskog diffusion coefficient, (kinematic) shear viscosity, and (kinematic) bulk viscosity, respectively. 142n Eqs.(D8) to (D10), n 0 is the equilibrium number density, σ the hard-sphere diameter, and b ≡ (2/3)πσ 3 is the second virial coefficient for a hard-sphere fluid.D 0 and ν 0 are the values of the self-diffusion and shear viscosity coefficients in the limit of zero density, and g(σ) is the radial distribution function evaluated at the hard-sphere point of contact.Analytic approximations for g(σ) can be obtained from the Percus-Yevick equation, scaled particle theory, [143][144][145] or the Carnahan-Starling equation of state; [146][147][148] we used the Carnahan-Starling approximation for the 3D pair distribution: where the packing fraction φ ≡ (π/6)n 0 σ 3 .For all calculations in Fig. 3, we used a Gaussian static spectrum (p = 2) primarily as proof-of-principle, though one should expect higher values of p to apply only at the highest densities; e.g., in the liquid state, where propagating longitudinal and shear wave modes may be significant for softer intermolecular potentials.We also took σ = 3.405 Å (argon Lennard-Jones diameter), m = 40 Da for the (atomic) mass, and set D ⊥ = v 2 0 /υ (as discussed above) and D = 4µ/3 (as in the liquid case).We determined v 2 from the adiabatic index, i.e., v 2 /v 2 0 = C P /C V , using heat capacity values obtained from the NIST Chemistry WebBook for each combination of temperature and density, T and ρ 0 , where ρ 0 = mn 0 .In Ref. 34, MD simulations are performed for a supercritical Lennard-Jones (LJ) fluid at a reduced density n * 0 = n 0 σ 3 = 0.5 and reduced temperature T * = k B T /ǫ = 1.5, where σ and ǫ are the LJ diameter and energy, respectively.To compare with the results of Ref. 34, we took m = 40 Da, σ = 3.405 Å, and ǫ 120 K, so that the (dimensional) mass density and temperature are, respectively, ρ = mn 0 = 841.3kg/m 3 and T = 180 K; note that, for this temperature-density combination, P ≈ 19.2 MPa.For v , we again use the adiabatic index computed from heat capacities obtained from the NIST Chemistry WebBook, SRD 69 for isothermal properties of argon (i.e, v 2 /v 2 0 = C P /C V ).To obtain base values for the collision frequencies, i.e., υ 0 , γ 0 , etc., we used the same Enskog diffusion and viscosity coefficient formulas given in Appendix D. However, we chose to treat this comparatively dense supercritical (SC) state slightly differently, as the packing fraction is φ ≈ 0.262-more than half of argon's triple point (TP) density (φ ≈ 0.445).In particular, we employed a hybrid treatment (detailed below) wherein the molecular scale a was computed via the second-moment condition [req.(5)  and Eq.(B7)] using simple physical arguments to obtain an estimated Einstein frequency of ω E = 4. ω E , is physically relevant for SC argon at the given density because the surface-to-surface spacing is (still) substantially smaller than the particle diameter; i.e., we assume that particles cannot easily "break through the cage" formed by its neighbors at this (packing) density.Second, to arrive at an estimate for ω E , we assume that the larger interparticle spacing in the SC case leads to an increased oscillation period, ∆τ , over that of TP argon, where ∆τ = τ SC − τ TP , and τ SC ≡ 2π/ω SC E and τ TP ≡ 2π/ω TP E are the respective oscillation periods.Finally, to estimate the increase in oscillation period, we make the simple physical assumption that ∆τ is an intervening ballistic interval arising from the increased surface-to-surface distance, an additional separation of d SC − d TP ≈ 0.695 Å over the TP case.This leads to ∆τ ≈ 2 d SC − d TP /v b , where the factor of 2 accounts for two ballistic traversals cage oscillation and we take the ballistic speed to v b ≈ √ 2k B T /m = √ 2v 0 (i.e., the most probable speed of a Maxwell-Boltzmann distribution).
Putting everything together, we arrive at a rough estimate for the oscillation period for the dense SC fluid: for dilute fluids (orange) and explicitly setting a → a 0 (green).Note that the blue curve reproduces the rather subtle "plateau" region (t * ≈ 0.15-0.6)due to the soundlike mode originating from the longitudinal component, as well as the nuanced transition from exponential-like to diffusive decay for t * 1.

( 7 )
long-time diffusivity b Ds = N ∞ 0 dt ∞ 0 dkk 2 q(ka)C ⊥ (k, t)Fs(k, t) Zero-frequency sum rule from Green-Kubo relation.f a General requirements that constrain the functional forms of correlation functions.b Specific requirements that constrain the values of correlation function parameters.c Requires D ⊥ , Ds = 0 and k-independence of C ⊥ and Fs for large k.See Appendix D. d Sum rule constraint is valid when Einstein frequency, ωE, is well defined, e.g, for liquids.See Appendix B. e Implies the invariant forms a 2 υ and a 2 γs.See Eqs.(32) to (33) and Eqs.(D6) and (D7).f Liquids: C ≈ 0, Fs ≈ 1 yields analytic solution for D ⊥ .Gases: C ≈ 0 yields constraint for γs, D ⊥ .See Appendices C and D.

( 1 )
Presented a new derivation and interpretation of the hydrodynamic VACF formulation, Eqs.(2) to (5).(2) Established a core set of physical desiderata whose constraints are sufficient to recover realistic VACFs.(3) Identified q(ka) as the initial condition of the velocity covariance function that characterizes molecularscale kinetic fluctuations and reproduces the VACF over all timescales by judiciously superimposing each hydrodynamic mode.(4) Proposed linear PDEs (R10) to model hydrodynamic collective modes, where regularization and the zerofrequency sum rule effect moment closure; yields analytic solutions for C ⊥ and C with sufficiently rich structure to resolve subtle details in the VACF.(5) Proposed a new hydrodynamic model for the selfdensity, leading to a viable analytic form of the selfintermediate scattering function, F s (k, t), for all densities; captures exponential decay at low densities.(6) Described the (re-)scaling of hydrodynamic model

Figure C. 1
compares the magnitude of χ(p) for the generalized Gaussian (blue) and generalized Lorentzian (orange) forms for the static spectrum.It is seen, for instance, that a generalized Gaussian with p = 4 gives roughly the same contribution to the overall diffusivity as the generalized Lorentzian with p ≈ 6.
75 ps −1 .The resulting VACF-Fig.E.1, blue-is in striking agreement with the MD calculations shown in Fig. 1 of Ref. 34.Our estimate for the Einstein frequency is based on the following calculations and physical reasoning.First, at φ ≈ 0.262, the mean interparticle separation (center-to-center) is d SC 0 ∼ n −1/3 0 ≈ 4.29 Å, so the contact (surface-to-surface) distance, d SC , may be estimated as d SC ≈ d SC 0 − σ ≈ 0.885 Å; likewise, d TP ≈ 0.190 Å at the TP.It is within reason to assume the Einstein frequency,

Table I .
Desiderata and constraints for physically admissible correlation functions.

Table II .
Key parameters.Left two columns: relations and constraints for molecular parameters.Right three columns: input data for liquid rubidium/argon calculations in Figs. 1 and 2. /Ds Ds 2.85 • 10 −9 m 2 /s 91 1.75 • 10 −9 m 2 /s 92 Finally, we use a generalized Gaussian spectrum with p ≈ 1.5 for which Eq. (B7) gives a ≈ 2.28 Å and yields excellent agreement (Fig. E.1, blue) throughout the timescales sampled by MD.Note that time is expressed in LJ reduced units, t * ≡ t/τ LJ , where τ LJ ≡ mσ 2 /ǫ is the characteristic LJ timescale; for argon, τ LJ ≈ 2.156 ps. Figure E.1 also shows VACFs computed using two different (larger) values of a: computing a from Eq. (D7)