Heat capacity of liquids: A hydrodynamic approach

We study autocorrelation functions of energy, heat and entropy densities obtained by molecular dynamics simulations of supercritical Ar and compare them with the predictions of the hydrodynamic theory. It is shown that the predicted by the hydrodynamic theory single-exponential shape of the entropy density autocorrelation functions is perfectly reproduced for small wave numbers by the molecular dynamics simulations and permits the calculation of the wavenumber-dependent specific heat at constant pressure. The estimated wavenumber-dependent specific heats at constant volume and pressure, $C_{v}(k)$ and $C_{p}(k)$, are shown to be in the long-wavelength limit in good agreement with the macroscopic experimental values of $C_{v}$ and $C_{p}$ for the studied thermodynamic points of supercritical Ar.


Introduction
Thermodynamics of the liquid state of matter was successfully explained on the basis of thermodynamic perturbation theory (TDPT) more than 40 years ago [1][2][3][4]. Together with the integral equation theory of static structure of liquids and computer simulations, it is the basis of current exploration of physics and chemistry of real liquid systems. Modern versions of TDPT [5,6] were recently developed on the basis of treatment of the interaction interatomic potential as a sum of an improved reference system (hard-sphere-or soft-sphere-like repulsion plus short-range Yukawa-like attraction) and a perturbation due to the remaining tail of the interparticle interaction.
It is also possible to represent the thermodynamics of the liquid state via time-dependent processes because dynamic processes in condensed matter systems are tightly connected with their thermodynamics. However, in contrast to solids, where thermodynamics can be easily defined by contributions from static lattice and small vibrations around the lattice sites in terms of phonon excitations [7], in fluids such an approach is not applicable due to the dynamic disorder and the absence of stable local energy minima for atomic particles. Furthermore, in liquid state there exist various relaxation processes which are absent in the solid state. As an example one can consider a relaxation connected with diffusivity of local temperature. In liquid system without any external temperature gradient, this relaxation process is caused by local temperature gradients which emerge due to adiabatic propagation of the most longwavelength acoustic excitations. The acoustic modes, which propagate with small wavenumbers due to local conservation laws, in fact create adiabatically compressed and decompressed local regions in liquid with different local temperatures. Such a complex origin of collective dynamics in fluids which is completely different from the elastic sound propagation in solids is tightly connected with the phenomenon of viscoelasticity of liquids [8]. Furthermore, the acoustic excitations in liquid state are only of longitudinal polarization (transverse long-wavelength acoustic modes are not supported in liquids, there can exist only short-wavelength shear excitations with a propagation gap in the long-wavelength region) and have non-zero damping. The dispersion law of longitudinal and transverse collective excitations in liquids cannot be considered as the linear in wave numbers [9,10] as this is taken in the Debye model of the heat capacity of solids. All these facts must be reflected in the construction of liquid thermodynamics represented via the dynamic processes.
Heat capacity being one of the important quantities in thermodynamics of the liquid state of matter can be easily calculated from the temperature fluctuations in the liquid system [11]. This is by date the most simple way of estimation of the specific heat at constant volume C v for realistic liquids from computer simulations. Schofield [12,13] has shown how the generalized wavenumber-dependent thermodynamic quantities can be estimated from static correlators of different statistically independent dynamic variables [14]. The generalized wavenumber-dependent thermodynamic quantities were studied later on for pure [15][16][17] and binary liquids [18] in classical and ab initio [19] molecular dynamics simulations.
Perhaps the first application of the hydrodynamic theory to collective dynamics of liquids was reported in the paper by Landau and Placzek [20], while systematic theoretical studies of liquid dynamics with hydrodynamic equations have started since the seminal studies by Mountain [21]. In 1971 Cohen et al. derived analytical expressions for hydrodynamic time correlation functions for pure and binary liquids [22], which contained the so-called "non-Lorentzian corrections" and satisfied all the sum rules existing within the hydrodynamic treatment [23]. Molecular dynamics (MD) simulations being a perfect tool for exploration of time-dependent correlations in realistic liquids and solids enable one to check the predictions of the hydrodynamic theory.
In this study we report the behaviour of autocorrelation functions of energy, heat and entropy densities obtained by molecular dynamics simulations of supercritical Ar at two densities. By date there were no MD studies of the entropy density autocorrelation functions of pure liquids, for which the hydrodynamics predicts purely single-exponential decay with time in the long-wavelength region. Besides, the zero-time value of the entropy density and heat density autocorrelation functions should be connected with the wavenumber-dependent specific heats C p and C v , respectively. We will perform a check of the long-wavelength limits of the C p (k) and C v (k) with the macroscopic experimental values for these thermodynamic quantities [24]. The remaining part of the paper is organised as follows. In the next section we give the predictions of the hydrodynamic theory on the studied time correlation functions. In section 3 we supply details of the MD simulations, and section 4 reports obtained the results on the autocorrelation functions of energy, heat and entropy densities obtained by molecular dynamics simulations of supercritical Ar, their corresponding time-Fourier transforms, and the estimated dependencies for the generalized specific heats C p (k) and C v (k). The last section contains conclusions of this study.

Hydrodynamic approach
Hydrodynamic approach to collective dynamics of liquids is valid only on macroscopic length and time scales where the most long-time relaxation processes in continuum are considered. The macroscopic equations in the hydrodynamic theory are in fact the local conservation laws for the densities of particles, of total momentum and of energy. For the case of longitudinal dynamics, the only relevant on macroscopic scale contributions to the hydrodynamic time correlation functions [8] are the ones coming from acoustic longitudinal propagating modes with the linear in wave numbers k dispersion law Here c s is the adiabatic speed of sound and the damping coefficient Γ = (D L +(γ−1)D T )/2 depends on the longitudinal kinematic viscosity D L , the thermal diffusivity D T and the ratio of specific heats γ = C p /C v .
The only relevant relaxing mode for the longitudinal dynamics on macroscopic scales is the one coming from thermal relaxation, i.e., connected with the diffusivity of local temperature D T .
The analytical expressions for the hydrodynamic time correlation functions, which we will study by MD simulations, have in general a three term form: an exponential function connected with thermal relaxing mode and two oscillating terms coming from the acoustic excitations [22]. The heat density autocorrelation function reads where the amplitudes of mode contributions are constrained by the lowest-order sum rule Note that the mode contributions A hh and B hh are constants in the long-wavelength limit, while the amplitide of the "non-Lorentzian" term D hh (k) linearly depends on k. However, outside the hydrodynamic region when one applies MD simulations the amplitudes A hh and B hh become k-dependent making possible estimation of the wave-number dependent specific heat. Similar is the expression for the energy density autocorrelation function: while for the entropy density autocorrelation function, the hydrodynamics predicts a single exponential form for pure liquids: Hence, the wavenumber-dependent specific heats C v (k) and C p (k) are connected with the static heatdensity and entropy-density autocorrelation functions, which in fact are the zeroth moments of the corresponding spectral functions S hh (k, ω) and S ss (k, ω). In the long-wavelength limit, the wavenumber dependent specific heats tend to their macroscopic values C v and C p . The advantage of the hydrodynamic approach lies in a possibility to estimate the contributions from relaxing and propagating modes to corresponding thermodynamic quantities [25].

Molecular dynamics simulations
We performed molecular dynamics simulations for supercritical Ar at temperature T = 280 K and two densities 750 kg/m 3 and 1464.5 kg/m 3 using a system of 2000 particles interacting via the Woon potential [26] with the parameters taken from [27]. The cut-off radius for the effective two-body potential was 12 Å. This potential permits a nice reproduction of the thermodynamic and dynamic quantities of supercritical Ar by MD simulations as it was shown in [28].
The simulations were performed in the microcanonical ensemble with the time step of 2 fs. The duration of production runs was of 360 000 time steps. Every sixth configuration was used for the sampling of dynamic variables of particle density, momentum density, energy density and first time derivative of the momentum density. The averages of static and time correlation functions over all possible directions of different wave vectors with the same wave number were performed.
During the MD simulations we sampled the spatial Fourier components of hydrodynamic variables of: particle density

13606-3
and energy density e(k, t ) = 1 where r i (t ), v i (t ) and ε i (t ) are the trajectories, velocities and single-particle energies of the i -th particle [14]. The single-particle energies were calculated as follows: where Φ(r ) is the effective pair potential used in the simulations.
The above dynamic variables were used to calculate the hydrodynamic time correlation functions as well as the heat density autocorrelation functions where the brackets mean the ensemble average and static averages read as follows f ne (k) = 〈n(k)e(−k)〉 and f nn (k) = 〈n(k)n(−k)〉.
In order to study the entropy density autocorrelation functions F ss (k, t ), we sampled the density of the longitudinal component of stress tensor σ L (k, t ) via [14] J where the overdot denotes the time derivative of the corresponding dynamic variable. Hence, the dynamic variable of entropy density can be sampled from MD simulations as [14] s with the corresponding static averages f eσ L (k) = 〈e(k)σ L (−k)〉 and f σ L n (k) = 〈σ L (k)n(−k)〉.

Results and discussion
Energy density autocorrelation functions A typical feature is an increasing amplitude and a decreasing frequency of the damped oscillations in the shape of F ee (k, t ) and F hh (k, t ) towards smaller wave numbers, which is typical of the contribution coming from acoustic excitations. The relative strength of the relaxing and oscillating contributions in F hh (k, t ) was estimated in the long-wavelength limit in [25] as the inverse of (γ − 1), i.e., the inverse of the Landau-Placzek ratio known for the density-density time correlation functions. Hence, for the case of the two thermodynamic points of the supercritical Ar studied here and having the experimental ratio of specific heats γ of 2.30 and 1.55 [24] for the low-and high-density studied systems, one should observe much smaller relaxing contribution for the low-density state. Indeed, in figure 2 (a) the oscillations in the shape of F hh (k, t ) even make the heat density autocorrelation function negative close to ωτ ∼ 1 for the smallest possible in these simulations wave number, while in figure 2 (b), due to the much stronger contribution from the relaxing thermal mode, the F hh (k, t ) is positive for all times. It is important that the dimensionless functions F hh (k, t ) have the zero-time values tending to the macroscopic values of C v due to (2.2). Indeed, in figure 2 (a) the zero-time values tend to the macroscopic value of C v = 1.775 [24], while for the high-density state they tend to the value C v = 2.34. The time-Fourier transformed energy density (e-e), heat energy density (h-h) and entropy density (s-s) autocorrelation functions, S i i (k, ω), i = e, h, s are shown for the two studied densities of the supercritical Ar in figure 4. The S ee (k, ω) and S hh (k, ω) have the typical of the dynamic structure factors shape with the central peak due to thermal relaxation and the side peak due to propagating acoustic excitations. The ratio of the integral intensities of the central and two side peaks known as the Landau-Placzektype ratio was obtained in [25] for the case of S hh (k, ω) and is equal to (γ − 1) −1 . The side peak of the S ee (k, ω) and S hh (k, ω) is the consequence of the corresponding oscillating contributions to the F ee (k, t ) and F hh (k, t ), respectively. It is remarkable that the calculated spectral function S ss (k, ω) shows only the one-peak shape with the central peak due to thermal relaxation (line connected star symbols in figure 4).  The zero-time values of the time correlation functions F hh (k, t ) and F ss (k, t ) permit estimation of the wavenumber-dependent specific heats C v (k) and C p (k), which in the long-wavelength limit tend to their macroscopic values. In figure 5 we show the wavenumber dependencies of C v (k) and C p (k) for the two densities of the supercritical Ar at temperature 280 K. The asterisks at k = 0 show the values of macro-scopic values of C v obtained from the standard expression via the thermal fluctuations in the fluid system [11]. One can see that the reported wave-number dependencies tend quite well to the experimental values for the supercritical Ar taken from the NIST database [24]: C v =1.775k B and C p =4.097k B for the density of 750 kg/m 3 and C v =2.34k B and C p =3.62k B for the density of 1464.5 kg/m 3 .

Conclusions
In the MD simulations we calculated the autocorrelation functions of energy density, heat density and energy density and compared them for the smallest available wave numbers with the predictions of the hydrodynamic theory. Our results give evidence of a nice correspondence of the hydrodynamic description of time correlation functions and their zero-time values, which made it possible to estimate the wavenumber-dependent specific heats C v (k) and C p (k), with the MD-based estimation of their macroscopic values and experimental values for specific heats of supercritical Ar taken from the NIST database.
These results are very important from the viewpoint of analysing the contributions from thermal relaxation and collective excitations to the specific heats of liquids. Our results and the hydrodynamic approach make evidence of the definitive role of the thermal relaxation in calculation of the specific heats of fluids when it is estimated via dynamic processes. These results will urge a further understanding of the dynamics of the supercritical state of matter, which shows many interesting features observed recently in the scattering experiments and MD simulations on supercritical Argon [29,30] and Oxygen [31].