The Rosenfeld-Tarazona expression for liquids' specific heat: A numerical investigation of eighteen systems

We investigate the accuracy of the expression of Rosenfeld and Tarazona (RT) for the excess isochoric heat capacity, C_V^{ex} \propto T^{-2/5}, for eighteen model liquids. Previous investigations have reported no unifying features of breakdown for the RT expression. Here liquids with different stoichiometric composition, molecular topology, chemical interactions, degree of undercooling, and environment are investigated. We find that the RT expression is a better approximation for liquids with strong correlations between equilibrium fluctuations of virial and potential energy, i.e., Roskilde simple liquids [Ingebrigtsen et al., Phys. Rev. X 2, 011011 (2012)]. This observation holds even for molecular liquids under severe nanoscale confinement, the physics of which is completely different from the original RT bulk hard-sphere fluid arguments. The density dependence of the specific heat is predicted from the isomorph theory for Roskilde simple liquids, which in combination with the RT expression provides a complete description of the specific heat's density and temperature dependence.

Fundamental theories for the pressure (or density) and temperature dependence of thermodynamic quantities have gained renewed attention in the last decade. These theories can serve as a valuable input to equations of state 1,2 , but also as input to scaling strategies which relate key dimensionless transport coefficients to thermodynamic quantities, such as the excess entropy 3,4 (with respect to an ideal gas) or the excess isochoric heat capacity 5 . Predicting dynamical quantities from first principles is a challenging task. One such theory is mode-coupling theory 6 (MCT) which relates the dynamic density correlations of a fluid to its static structure. However, MCT becomes less reliable as the characterizing features of supercooled liquids become pronounced [7][8][9] . Extensions of MCT to more complex environments, such as nanoscale confinement or homogenous shear flow, have been slow to emerge [10][11][12][13][14] . Recently, however, other promising MCT-type-theories have started to appear for supercooled liquids 15 .
Alternatively, theories which relate dynamics to thermodynamics such as that of Adam and Gibbs 16,17 and Rosenfeld's excess entropy scaling 3,18 , relate to the dynamics also in the highly supercooled liquid regime 19,20 . Excess entropy scaling strategies, as proposed by Rosenfeld, have proven succesful in predicting the dynamics of not only single-component atomic fluids 3,4 , but also binary mixtures 4,19-21 , ionic substances 22,23 , small molecules [24][25][26][27][28] , and polymers 29 . In fact, excess entropy scaling strategies have been used as a reliable predictor even for the perplexing dynamics of nanoconfined liquids 5,30-35 which exhibit stratification and positiondependent relaxation processes; the applicability of excess entropy scaling strategies to these situations emphasizes the usefulness of the thermodynamic approach to dynamics.
To fully harness the power of predicting dynamics from thermodynamics, however, it is imperative to develop reliable theories for the pressure (or density) and temperature dependence of thermodynamic quantities. Rosenfeld and Tarazona 36 (RT) argued for a mathematically simple expression for the density and temperature dependence of the potential energy U for fluids. Their arguments are based on thermodynamic perturbation theory, using the fundamental-measure reference functional for hardspheres in combination with an expansion of the free energy around the η = 1 packing fraction. The arguments are involved and not easy to follow, but their expressions have found widespread application (Refs. 1,2,25,26,. From the potential energy one gains access to thermodynamic quantities such as the excess isochoric heat capacity C ex V = (∂U/∂T ) V and the excess entropy S ex via C ex V = T (∂S ex /∂T ) V . Both of these quantities enter the aforementioned strategies. For a long time only few studies focused on the heat capacity 58 . Recently, however, the heat capacity of liquids has started to receive more attention [59][60][61][62][63] . The RT expressions for the potential energy and excess isochoric heat capacity read where α(ρ) and β(ρ) are extensive functions of density ρ that relate to the specific system 36 . Several numerical investigations have tested the applicability of the RT expressions for various model liquids. These liquids span from simple atomic model fluids to liquids showing a wide range of structural, dynamical, and thermodynamical anomalies in their phase diagram. More specifically, the RT expressions have been investigated for single-component atomic fluids [36][37][38][39] , binary mixtures 2,40-47 , ionic substances 1,48-50 , hydrogenbonding liquids 25,51-55 , small molecules 26,56 , and sheared liquids 57 . These investigations showed that the RT expressions give a good approximation for a range of systems, but are less accurate when applied to systems known to not have strong virial/potential energy correlations 64 , such as the Dzugutov liquid and Gaussian core model, as well as for SiO 2 and BeF 2 in their anomalous regions. For SPC/E water different results for the applicability have been reported 25,51,54 .
The purpose of this paper is to provide insight into the conditions under which RT applies by investigating 18 different model systems possessing different stoichiometric composition, molecular topology, chemical interactions, degree of undercooling, and environment. We use GPU-optimized 65 NVT molecular dynamics computer simulations 66-69 (in total over 40000 GPU hours) to calculate the potential energy and excess isochoric heat capacity along a single isochore for each of these 18 model systems (for the single-component Lennard-Jones (SCLJ) liquid we also vary the density). Here and henceforth quantities are reported in dimensionless units, e.g., by setting σ = 1, ǫ = 1, etc. The heat capacity is calculated via Einstein's fluctuation formula C ex V = (∆U ) 2 /k B T 2 (within statistical error, the numerical derivative of the potential energy agrees with the value calculated from the fluctuations). Previous investigations have often calculated the excess entropy, but this is computationally very demanding and not necessary to test Eqs. (1) and (2). Table I presents the investigated model systems, which range from simple atomic fluids to molecules under severe nanoscale confinement. The densities represent typical liquid-state densities.   (3)) for the potential energy and excess isochoric heat capacity, respectively, for the isochore of density ρ. The virial/potential energy correlation coefficient R is given for the lowest temperature state point Tmin. The abbreviations used are: Kob-Andersen binary Lennard-Jones mixture (KABLJ); inverse power-law fluid with exponent n (IPL n); LJ polymer chain of length n (LJC n); Lewis-Wahnström o-terphenyl (OTP); singlecomponent Buckingham liquid (SCB); single-component LJ liquid (SCLJ); Wahnström binary LJ mixture (WABLJ). The "Nanoconfined dumbbell" is confined to a (smooth) slit-pore of width H = 8.13, corresponding to roughly 16 molecular lengths. Figures 1(a) and (b) show, respectively, the potential energy and excess isochoric heat capacity at constant density as a function of temperature for all investigated systems. In all cases, the excess isochoric heat capacity decreases with increasing temperature. The data points of Fig. 1 were generated by the following procedure.
1. First the system is cooled at constant density until one of the following happens: a) The system crystallizes; b) The pressure becomes negative; c) The relaxation time is of the order 10 5 time units. This happens at the temperature T min . The system is then equilibrated at T = T min ; in the case of crystallization or negative pressure, the temperature is increased slightly (and this new temperature defines T min ).
2. Next, the temperature is increased from T min up to T max = 3T min , probing state points along the isochore with a spacing of ∆T = (T max − T min )/7. A total of eight equilibrium state points are hereby generated for each isochore.
Turning now to the RT expressions, we show in Figs. 2(a) and (b) the coefficient of determination 81 D for the potential energy and excess isochoric heat capacity as a function of 1 − R (see below). For a generic quantity X, the coefficient of determination D X is defined by where f (X i ) is a function that provides the model values, and the average X is taken over a set of data points X = {X 1 , ..., X N }. In our case f (X i ) is given by fits to the data points in X using, respectively, U = A 0 T 3/5 + A 1 , and C ex V = 3/5A 2 T −2/5 , where A 0 , A 1 , and A 2 are constants. D X measures the proportion of variability in a data set that is accounted for by the statistical model 81 ; D X = 1 implies perfect account of the variability.
The virial/potential energy correlation coefficient R is and calculated from the canonical ensemble equilibrium fluctuations at T min . A new class of liquids was recently proposed, namely the class of strongly correlating liquids. These liquids are simple in the Roskilde sense of the term 82 and defined by R ≥ 0.90. Only inverse power-law fluids are perfectly correlating (R = 1), but many models 64 as well as experimental liquids 83 have been shown to belong to the class of Roskilde simple liquids. This class is believed to include most or all van der Waals and metallic liquids, whereas covalently, hydrogenbonding or strongly ionic or dipolar liquids are generally not Roskilde simple 64 . The latter reflects the fact that directional or competing interactions tend to destroy the strong virial/potential energy correlation.
By plotting D as function of 1 − R, we investigate whether a correlation between the "simple" expression of RT for the specific heat and "Roskilde simple liquids" exists. For many of the investigated systems that fall into this simple class such a correlation is not expected from the original RT bulk hard-sphere fluid derivation (see below).  Table I). The definition of the symbols is given in Fig. 1. (a) DU . (b) D C ex V . The Girifalco system gives a negative value for D C ex V and has for clarity of presentation been left out (see Table I). For both the potential energy and the excess isochoric heat capacity, RT is seen to deteriorate as R decreases below 0.90 (to the right of the red line); in particular see insets for the SCLJ liquid.
We observe from Fig. 2(a) that for all liquids D U gives a value close to 1, but RT provides a better approximation for liquids with R larger than 0.90 (to the left of the red line). A similar behavior is observed for D C ex V in Fig. 2(b), however, note the change of scale. The insets of both figures show for the SCLJ liquid how RT deteriorates as R decreases below 0.90. We conclude that the RT expressions work better for systems that are Roskilde simple at the state points in question.
Originally 36 , RT was argued from fundamentalmeasure thermodynamic perturbation theory for bulk hard-sphere fluids and via simulation shown to describe inverse power-law systems to a high degree of accuracy. Later investigations showed that RT is a good approximation also for LJ liquids. These systems are Roskilde simple, and a recently argued quasi-universality 84 for Roskilde simple single-component atomic systems implies this behavior.
We have shown that the key determining factor for RT is not whether systems are atomic or molecular (see results for dumbbell, OTP, and LJC), but rather the degree of strong correlation between virial and potential energy. This was shown to be the case even for severely nanoconfined molecular systems which exhibit a completely different physics from bulk hard-sphere fluids 5 and are thus not expected to satisfy the original RT arguments. The latter is, in particular, true also for the elongated nonspherical molecules studied here. The observed correlation between RT and Roskilde simple liquids is thus highly nontrivial.
As a further validation of the above viewpoint, we relate the function α(ρ) in the RT expression to h(ρ) for Roskilde simple liquids. For such a liquid, temperature separates 77 into a product of a function of excess entropy per particle and a function of density via T = f (s ex )h(ρ). Roskilde simple liquids are characterized by having isomorphs to a good approximation 85 . Isomorphs are curves in the thermodynamic phase diagram along which structure and dynamics in reduced units, as well as some thermodynamic quantities are invariant. Along an isomorph both C ex V and h(ρ)/T are invariant, and consequently one may write Since by the RT expression; C ex V = 3/5α(ρ)T −2/5 = 3/5 α(ρ) 5/2 /T 2/5 , it follows that h(ρ) = α(ρ) 5/2 or, equivalently, For a LJ system, it was shown in Refs. 77 and 86 that h(ρ) is given by in which γ 0 is calculated from the virial/potential energy fluctuations at ρ = 1 and T = 1 via γ 0 = ∆W ∆U / (∆U ) 2 .
Equation (6) is tested in Fig. 3 for the KABLJ and the repulsive LJ system (for which, respectively, γ 0 = 5.35 and γ 0 = 3.56). The repulsive LJ system is defined from v(r) = (r −12 + r −6 )/2 and has R above 99.9% in its entire phase diagram; γ varies from 2 at low density to 4 at high density. We determine α(ρ) for different densities by fitting Eq. (1) as a function of temperature for each isochore and system. h(ρ) is calculated analytically from Eq. (7). Figure 3 shows that α(ρ) as predicted by the isomorph theory to a very good approximation is given by h(ρ) 2/5 . A complete description is thus given via Eqs.
(2) and (7) for the density and temperature dependence of the specific heat, i.e., C ex V = (h(ρ)/T ) 2/5 . Scaling strategies which relate dynamics to thermodynamics have in the past proven useful to predict perplexing dynamical phenomena. We identified here the range of applicability for RT as the class of Roskilde simple liquids. By combining the RT expressions with the isomorph theory, we were able to provide also the full density and temperature dependence of the specific heat. The predictive power of the aforementioned scaling strategies for most van der Waals and metallic liquids is hereby significantly increased.  7)) for the KABLJ and repulsive LJ system (see text). The red and blue curves are proportional to h(ρ) 2/5 with the proportionality constant determined from the highest density state point (ρ = 4.00) for each system.