Pressure-energy correlations in liquids. III. Statistical mechanics and thermodynamics of liquids with hidden scale invariance

In a series of papers published last year,^1–5 we introduced the concept of strongly correlating liquids and demonstrated by computer simulations that this includes a large class of model liquids. Specifically, the fluctuations, which are in many cases strongly correlated, are those of the configurational parts of pressure and energy, i.e., the parts in addition to the ideal gas terms, coming from the interatomic forces. Recall that for any microscopic state, energy $E$ and pressure $p$ have contributions both from particle momenta and positions,

For a liquid with pair interactions, if $v ( r )$ is the pair potential and $r i j$ is the distance between particles $i$ and $j$⁠, we have

Strong $W$⁠, $U$ correlation, if present at all, is observed under conditions of fixed volume, as illustrated in Fig. 1(a). The degree of correlation is quantified by the standard correlation coefficient $R$⁠, defined^1,4 by

Here and henceforth, unless otherwise specified, angular brackets ⟨ ⟩ denote thermal $N V T$ ensemble averages and $Δ$ denotes deviation from the average value of the quantity in question. We call liquids with $R > 0.9$ strongly correlating. Another characteristic quantity is the “slope” $γ$⁠, which we here define^1,4 as the ratio of standard deviations,

In the limit of perfect correlation $( R → 1 )$⁠, $γ$ becomes equal to the slope of average $W$ as a function of average $U$ at fixed volume, $γ= ( ∂ W / ∂ U ) V$⁠.

In Paper I (Ref. 4) of this series, it was shown that strongly correlating liquids are typically those with van der Waals type attraction and steep repulsion, which in simulations are often modeled by combinations of one or more Lennard-Jones (LJ)-type potentials. Typical slope values for the latter are of the order of 6, depending slightly on the state point (in the limit of very high density or temperature, the slope converges slowly to 4).

It is worth noting that the class of strongly correlating liquids does not simply correspond to radially symmetric pair potentials. Firstly, two metallic systems with many-body potentials were found to be strongly correlating; we believe metallic systems in general are strongly correlating, although this needs to be confirmed. Also many molecular liquids are strongly correlating. In fact, any potential with an inverse power-law (IPL) dependence on distances (not necessarily based on pair interactions) is perfectly correlating. Secondly, there exist radially symmetric pair potentials that are not strongly correlating, for example, the Dzugutov system.^4,7 One reason for strong correlation not to hold in some molecular systems is the presence of Coulombic terms in the potential. By themselves these would give 100% correlation, but their combination with LJ forces typically leads to a weak correlation. This was detailed in Paper I, which presented results from simulations of 13 different model liquids. In our present understanding based on these simulations, liquids with two length scales in their potentials are rarely strongly correlating.

Paper II (Ref. 5) in this series analyzed the case of the standard single-component LJ liquid in detail. Building on the fact that IPL potentials $∝ r − n$ are perfectly correlating, the results of this analysis can be summarized as follows: (1) almost all of the fluctuations in $W$ and $U$ come from interparticle separations in the region of the first peak of the radial distribution function $g ( r )$⁠; (2) in this region, the LJ potential is approximated very well by the sum of an IPL with exponent $n ∼ 18$ and a linear term $B + C r$⁠; and (3) when volume is fixed, the parts of $W$ and $U$ that come from the linear term are almost constant. Our initial and simpler explanation of strong $W U$ correlations¹ was based on the dominance of close encounters, i.e., that it is only the nature of the repulsive part of the potential that matters for the strong correlations. This explanation, however, is only adequate at high pressure/density. It does not explain the requirement of fixed volume, nor the fact that strong correlation is observed even at zero pressure, as well as for the low-temperature/low-pressure (classical) crystal. To see that an explanation at the individual pair interaction level is generally inadequate, consider Fig. 1(b) which shows a scatter plot of single-particle energy and virial. These are sums over the pair interactions a given particle has with its neighbors; summing over all particles gives the total potential energy and virial, respectively. If strong correlation held at the level of single pair interactions, it would also hold at the particle level, but it clearly does not. This emphasizes that strong correlation is a collective effect, as detailed in Paper II.

In this paper, we elaborate on the statistical mechanics and thermodynamics of strongly correlating liquids and report results from computer simulations, showing that strong virial/potential energy correlations are present even in nonequilibrium processes. The purpose is to present a number of new results supplementing those of Paper II in order to broadly illuminate the properties of strongly correlating liquids. Together Papers II and III give a fairly complete characterization of the properties of a strongly correlating liquid at one state point, as well as at different state points with the same density. Paper IV⁸ in this series goes on to consider varying-density curves of “isomorphic” state points in the phase diagram, which are characterized by several invariants; such curves exist only for strongly correlating liquids. Paper IV also discusses how and why there are, in fact, other numbers close to the $γ$ of Eq. (6) that reflect a liquid’s hidden scale invariance. For any strongly correlating liquid, these numbers are close to each other, so for simplicity throughout this paper we shall only use one $γ$⁠, namely the one defined by Eq. (6) (in Paper IV the definition of $γ$ changes to a slighly smaller number).

The organization of the paper is as follows. Section II begins with a discussion of the scaling properties of systems with IPL potentials, the natural starting point for a discussion of the hidden scaling properties of strongly correlating liquids. This is followed by a generalization to allow an extra term in the free energy depending on the volume only. Some, but not all, scaling properties of IPL systems are inherited by this generalization. Following this, in Sec. III we discuss the “extended inverse-power law” (eIPL) potential introduced in Paper II, which includes the above-mentioned linear term. We illustrate with simulation results the key property that the linear term contributes significantly to the virial and potential energy fluctuations when the volume may fluctuate, but little when it is fixed. Hence, it gives rise, approximately, to a volume-dependent term in the free energy of the type discussed in the previous subsection. This leads to an inherited approximate scaling property, which we refer to as “hidden scale invariance” since it is not immediately obvious from the intermolecular potential. The argument about how and why hidden scale invariance causes strong $W U$ correlations makes no assumption about equilibrium. To emphasize this point, Sec. IV presents results from nonequilibrium computer simulations of strongly correlating molecular liquids, in particular, aging following a temperature jump, and crystallization, both at constant volume. The property of strong correlation is shown to apply even in these out-of-equilibrium situations. Section V discusses ensemble dependence; it is shown here that the virial/potential energy correlation is always stronger in the $N V T$ ensemble than in the $N V E$ one. The last main section, Sec. VI, comprises of two topics under the heading “thermodynamics of strongly correlating liquids.” First, we discuss the relation of pressure-energy correlations to the thermodynamic Grüneisen parameter $γ G$⁠, showing that the slope $γ$ [Eq. (6)] is larger than $γ G$ by roughly a factor involving the ratio of excess (configurational) to total constant-volume specific heat. This ratio is around 2 for many simple liquids.⁹ The second part formulates the property of strong correlation in the energy-bond language known as “network thermodynamics.”¹⁰ 

The purpose of this section is to summarize the properties of IPL potentials and identify which of these properties are inherited by strongly correlating liquids and which are not.

IPL potentials—sometimes referred to as soft-sphere potentials—have been used in liquid state theory for many years as convenient model systems.^11–18 Such potentials have a number of simple properties. IPL potentials have, however, been considered unrealistic because their predicted equation of state is quite unrealistic and because they have no stable low-pressure liquid phase and no van der Waals loop, problems that derive from the fact that IPL potentials are purely repulsive. Moreover, the correct IPL exponent fitting the LJ liquid is around 18 (Papers I and II, Refs. 16, 19, and 20), not 12 as one might naively guess from the repulsive $r − 12$ term of the LJ potential; this may have confused people searching from an effective IPL description of the LJ liquid. A major point made in this series of papers is that, when interpreted correctly, IPL potentials are much more realistic than generally thought because they describe well a number of properties of strongly correlating liquids. For reference, we now briefly summarize the long established properties of IPL liquids.

Consider $N$ identical particles in volume $V$ interacting by a pair potential of the form $v ( r ) = A r − n$⁠; we make the pair assumption for simplicity but note that the argument below generalizes immediately to any potential that is an Euler homogeneous function of the position coordinates. From the standard partition function, the (Helmholtz) free energy $F$ is conveniently written^6,21 as the ideal gas term plus the nontrivial “excess” free energy, $F = F id + F ex$⁠. The first term is the free energy of an ideal gas at the same volume and temperature, $F id = N k B T ln ( ρ Λ 3 / e )$⁠, where $ρ = N ∕ V$ is the particle number density and $Λ = h ∕ 2 π m k B T$ is the thermal de Broglie wavelength. The excess free energy is given^6,21 by

Whenever $n > 3$⁠, this expression leads to a free energy with a well-defined extensive thermodynamic limit $( N → ∞ )$⁠.^11,12

It follows from Eq. (7) that the excess free energy of an IPL liquid is given as follows in terms of a function of density $ρ$ to the power $n ∕ 3$ over temperature $T , ϕ ( ρ n ∕ 3 ∕ T )$ (Klein’s theorem^11,12):

Equation (8) implies that a number of derived quantities are also functions of $ρ n ∕ 3 ∕ T$⁠. As important examples, recall the following standard identities: the excess entropy: $S ex = − ( ∂ F ex ∕ ∂ T ) V$⁠, the potential energy: $U = F ex + T S ex$⁠, the virial $W = − V ( ∂ F ex ∕ ∂ V ) T$⁠, the excess isothermal bulk modulus: $K T ex = V ( ∂ 2 F ex ∕ ∂ V 2 ) T$⁠, the excess isochoric specific heat per unit volume: $c V ex = − ( T ∕ V ) ( ∂ 2 F ex ∕ ∂ T 2 ) V$⁠, the excess pressure coefficient: $β V ex = ( 1 ∕ V ) ( ∂ W ∕ ∂ T ) V = − ∂ 2 F ex ∕ ∂ T ∂ V$⁠. These quantities are all functions of the single variable $ρ n ∕ 3 ∕ T$⁠; more accurately, one has [where $f 1 ( x ) = x ϕ ′ ( x ) − ϕ ( x )$⁠, $f 2 ( x ) = x ϕ ′ ( x )$⁠, $f 3 ( x ) = ( n ∕ 3 ) 2 x 2 ϕ ″ ( x ) + [ ( n ∕ 3 ) + ( n ∕ 3 ) 2 ] x ϕ ′ ( x )$⁠, and $f 4 ( x ) = − x 2 ϕ ″ ( x )$],

The functions $ϕ$ and $f 1 , … , f 4$ all depend on $n$⁠, but, for simplicity of notation, we have not indicated this explicitly. Dividing across by the dimensional factors on the right hand side (for example, $k B T$ in the case of potential energy and virial), one arrives at dimensionless forms of the excess entropy, potential energy, etc., which are functions of $ρ n ∕ 3 ∕ T$ only.

Turning now to the dynamics, consider the standard molecular dynamics (MD) case where the equations of motion are Newton’s equations. Suppose $r i ( t ) ( i = 1 , … , N )$ is a solution to Newton’s equations. Then, it is straightforward to show that $r i ( 1 ) ( t ) = α r i ( λ t )$ is also a solution if $α − ( n + 2 ) = λ 2$⁠. In particular, if $r i ( t )$ refers to equilibrium (⁠$N V E$ or $N V T$⁠) dynamics at a state point with density $ρ 0$ and temperature $T 0$⁠, then $r i ( 1 ) ( t ) = α r i ( λ t )$ refers to equilibrium dynamics at density $ρ 1 = ρ 0 ∕ α 3$ and temperature $T 1 = T 0 α 2 λ 2$ (temperature scales as the mean-square velocity and velocities get a factor $α λ$⁠). Using the above relation between $α$ and $λ$⁠, this implies that

This means that two states with different densities and temperatures but the same $ρ n ∕ 3 ∕ T$ have dynamics that scale into one another by simple scalings of space and time. In particular, if for any quantity $A$ one defines the relaxation time $τ A$ via $τ A = ∫ 0 ∞ ⟨ A ( 0 ) A ( t ) ⟩ d t ∕ ⟨ A 2 ⟩$⁠, it follows that any two states with same $ρ n ∕ 3 ∕ T$ have the same “reduced” (dimensionless) relaxation time $τ ̃ A$⁠, if this quantity is defined by $τ ̃ A = τ A ∕ t 0$ where the characteristic time $t 0$ is defined by $t 0 = ρ − 1 ∕ 3 m ∕ k B T$⁠. Similarly, if one defines the reduced diffusion constant $D ̃$ by $D ̃ = D ∕ D 0$⁠, where $D 0 = ρ − 2 ∕ 3 ∕ t 0 = ρ − 1 ∕ 3 k B T ∕ m$⁠, then $D ̃$ is the same for the two states. Summarizing,

Finally, we note that it follows from the above scaling property that

Thus, the Boltzmann factors of the two microscopic configurations are the same. Consequently, the scaling of the dynamics holds also for stochastic dynamics. This observation, in a generalized form, is the starting point of Paper IV in this series. By the same argument, the structure of states with the same $ρ n ∕ 3 ∕ T$ are identical provided that lengths are scaled by $ρ − 1 ∕ 3$⁠.

This section discusses how a potential may have a number of IPL properties to a good approximation, thus justifying the term “hidden scale invariance.” Consider a general potential between particles $i$ and $j$⁠, rewriting it (as can always be done) as a sum of an IPL potential plus the difference, denoted as “diff,”

For any microscopic configuration $( r 1 , … , r N )$⁠, the potential energy is then the sum of an IPL term and a “diff” term, and the excess free energy is given by

We now investigate consequences of the crucial assumption (see Paper II and Sec. III) that $U diff$ to a good approximation is only a function of volume: $U diff ( r 1 , … , r N ) ≅ f ( V )$—at least for states that carry Boltzmann weights of any significance. The approximate identity $U diff ( r 1 , … , r N ) = f ( V )$ means that the second exponential can be moved outside the integral, and we get

From this follows directly that

This implies that

While the systems under consideration here have the same excess entropy as their “hidden” IPL systems, several other quantities have contributions from the $f ( V )$ term. These quantities do not obey IPL scaling. In contrast, the scaling behavior for dynamics and structure is inherited. To see this, consider two state points with the same $ρ n ∕ 3 ∕ T$⁠. For the pure IPL system $[ f ( V ) = 0 ]$⁠, the two state points have the same dynamics and structure as argued in the previous section. Letting $f ( V ) ≠ 0$ simply shifts the energy surface, which changes neither the dynamics nor the structure.^22,23 This scenario—scaling of the dynamics, but not the equation of state—is exactly what is experimentally observed for a large number of viscous liquids. For example, in van der Waals liquids relaxation times are found to be a function of $ρ n ∕ 3 ∕ T$ (using $n$ as an empirical parameter),^24,25 but the scaling does not apply to the (excess) pressure with the exponent determined from the scaling of relaxation time, as required for IPL scaling.^24,26

In Sec. III we provide numerical evidence that there are indeed systems that to a good approximation fulfill the assumption introduced above that $U diff$ is a function of volume only.

In this section we examine the extent to which the LJ potential may be approximated by an eIPL potential [including a linear term, Eq. (33), below] by considering the fluctuations at a particular state point of the LJ liquid. We choose a state point whose pressure is near zero because here it is particularly clear that single-pair effects are insufficient to explain the strong $W U$ correlation (Fig. 1).

The analysis of Paper II took its starting point in assuming that the approximating IPL should match the potential closely at a particular value of the interparticle separation. An important conclusion of the analysis was, however, that the success of the IPL approximation is derived not from its behavior near any particular $r$-value, but rather from the fact that the difference from the real potential is close to linear over the whole first-peak region. In this section we reexamine the idea that the fluctuations are well described by an IPL potential and the argument for why the difference term almost does not fluctuate. We show explicitly that the latter contributes little to the fluctuations at constant volume, but significantly when the volume is allowed to fluctuate as in the $N p T$ ensemble. This demonstrates that the LJ potential is of the type considered in the previous subsection.

We wish to determine to what extent the LJ potential

can be matched, for the purpose of describing fluctuations of potential energy and virial at fixed volume, by an IPL

Here, $B$ indicates an optional constant. To start with, how should the exponent $n$ and the coefficients $A$ and $B$ be chosen? An obvious choice, followed in the first part of Paper II, is to require that the two potentials, $v LJ$ and $v IPL$⁠, should agree as much as possible around a particular value of $r$⁠, denoted $r 0$⁠. Given $r 0$⁠, if we require the functions and their first two derivatives to match at $r 0$⁠, this determines all three parameters $A$⁠, $n$⁠, and $B$⁠. The exponent $n$ is given (Paper II) by

Here, the notation $n ( 1 ) ( r 0 )$ refers to one kind of $r$-dependent effective IPL exponent, based on the ratio of the second and first derivatives. For $v IPL ( r )$ this simply returns $n$⁠. Otherwise, it gives a local matching of the $v LJ ( r )$ to $v IPL ( r )$⁠. This leaves effectively one parameter to vary, namely $r 0$⁠, which must be less than the location of the minimum of the LJ potential $r m = 2 1 ∕ 6 σ$⁠, where $n ( 1 )$ diverges. The parameter $r 0$ may be chosen to optimize the match of the fluctuations in the total energy and virial. For an $N V T$ simulation at $T = 80 K$ and near-zero pressure, the best fit was obtained with $n = 19.2$ (while the exponent implied by the slope [Eq. (6)], $γ = 6.3$⁠, was slightly smaller, 18.9).

Later in Paper II, it was demonstrated that there is no particular reason why the potentials should match close at a particular value of $r$ since the fluctuations have contributions from the whole first-peak region, including beyond the potential minimum. The reason that any kind of matching is possible over this region—where $v LJ$ clearly does not resemble a power law—is that a linear term may be added to the power-law potential almost without affecting the fluctuations as long as the volume is held constant. The analysis of Paper II, which also included an in-depth treatment of the perfect LJ (fcc) crystal which is also strongly correlating, showed that the more relevant $r$-dependent effective exponent is the higher order $n ( 2 )$ defined by

This also returns $n$ for $v IPL ( r )$⁠, but since it does not involve the first derivative, it returns $n$ even if a linear function of $r$ is added to the potential as in the eIPL potential defined by

This potential fits the LJ potential very well around its minimum (Paper II) and thus includes some of its attractive part.

There are several possible ways of choosing the “best” eIPL to match the real potential. These will give slightly different exponents and coefficients $A$⁠, $B$⁠, and $C$⁠. We do not investigate them here; rather the purpose is to validate the basic idea of the eIPL approximation. Therefore, we choose a simple matching scheme, whereby we match the fluctuations to those of the IPL potential, without including a linear term, in order to determine $n$ and $A$⁠. For simplicity, we take the exponent directly from the observed fluctuations: $n = 3 γ$⁠, where $γ$ is defined in Eq. (6). To fix the coefficient $A$⁠, agreement with the potential energy and virial fluctuations is optimized by proceeding as follows. For a given configuration generated in a LJ MD simulation, we calculate the LJ potential energy $U LJ$ and the power-law potential energy $U IPL$⁠, similarly the corresponding virials $W LJ$ and $W IPL$⁠. The difference quantities $U diff$ and $W diff$ are defined as

A perfect match of the fluctuations would mean that $U diff$ and $W diff$ have zero variance. Therefore, we choose $A$ to minimize the sum of the relative diff variances,

For the near-zero pressure state point used in Fig. 1 of Paper I, the exponent determined from $γ$ is $n = 3 γ = 18.9$ and the optimal value of $A$ is $1.3437 ϵ σ n$⁠. Before examining the difference potential, what do we get if we simulate with the matched IPL potential? Figure 2 shows the radial distributions $g ( r )$ obtained for the above state point and another with a higher density and temperature, for three potentials: LJ, the repulsive $r − 12$ term of the LJ potential, and the optimal IPL potential with $n = 18.9$⁠. We used the same IPL potential at both state points (i.e., we did not adjust $A$ and $n$ to match the second state point). The first thing to note is that the $n = 18.9$ potential gives a structure much closer to that of the LJ potential than does the repulsive $n = 12$ term alone, in particular, the latter system has crystallized at the higher density and temperature. The second point is that there is still a difference between the LJ and the $n = 18.9$ IPL, present in both state points. The first peak in the LJ system is slightly higher and narrower, although its position is barely altered. Thus, the real potential gives a slight increase in order—however, the difference in the coordination number is less than 0.1 (integrating to the first minimum after the peak).

Figure 3(a) shows the LJ potential, the IPL potential with parameters optimized as described above, their difference, and the radial distribution function. As shown in Fig. 3(b), the main part of the difference potential $v diff ( r )$ is nearly linear. Thus, a good approximation to the real potential is the eIPL potential of Eq. (33) for $r$ less than a cutoff $r c$ and zero otherwise. Neglecting the small value of $v IPL ( r c ) ∼ 10 − 5 ϵ$⁠, the cutoff is given by $r c = − B ∕ C$⁠. For the fit shown in Fig. 3(b), $r c = 1.61 σ$⁠.

What are the implications of the linear term? A linear term in the pair potential contributes a term proportional to the sum of all bond lengths to the total potential energy. It was shown in Paper II that at constant volume, this sum is a constant in one dimension, and it was argued that it is approximately constant in three dimensions. The difference is because of two things. First, in three dimensions, there are contributions to bond-length changes from transverse components of relative displacements between the two particles defining a bond; second, within the eIPL the potential is only linear up to $r c$⁠, which means that argument that the linear term contributes little to the fluctuations depends on all bond lengths being less than $r c$⁠. In one dimension, at moderate temperatures, a single-component system has a rather well-defined nearest-neighbor distance, which, at densities where the pressure is not too negative, will be less than $r c$⁠. In a three-dimensional liquid, however, the nearest-neighbor distance is not as well defined—the radial distribution function does not go to zero after the first peak. Therefore, there will always be fluctuations at $r c$ as the lengths of bonds fluctuate back and forth across $r c$⁠, so the sum of bond lengths, which are less than $r c$⁠, will fluctuate. In Paper II, it was shown that for a three-dimensional (classical) crystal at low temperature—where this is not an issue because $g ( r )$ does go to zero after the peak—the correlation coefficient $R$ becomes very high, over 99.5%, as $T → 0$ (but not 100%).

We can check directly the effect of adding the linear term to the IPL. Figure 4 shows scatter plots of IPL (black) and eIPL (red) energy and virial plotted against the true LJ values. These were calculated for a set of configurations drawn from an $N V T$ simulation using the true (LJ) potential. Including the linear term makes little difference. It somewhat improves the match to the energies, though not to the virials (possibly due to the discontinuity in the pair virial at $r c$⁠). The $W U$ correlation coefficient of the eIPL potential energy and virial is 0.917 [compare to the true (LJ) value of 0.938 and the pure IPL value of 1.0].

It is instructive to repeat the above for configurations drawn from an $N p T$ simulation at the same state point, i.e., with pressure chosen as the average pressure of the corresponding $N V T$ simulation. The results are shown in Fig. 5. Here, it is clear that the IPL potential represents the potential energy fluctuations very poorly (black points), while adding the linear term makes a substantial difference (red points). As in the $N V T$ case, the linear term affects the virial fluctuations much less. This is presumably because when taking the derivative to form the virial, the IPL term gets multiplied by $n ⪢ 1$⁠, while the linear term gets multiplied by $− 1$ and is thus reduced considerably in significance (cf. the insets of Fig. 4).

The size of the variances of the different terms are compared in Table I. We do not make a detailed analysis of the variance (taking into account cross correlations, etc). In the $N V T$ case, the IPL contributions are of similar size to the full (LJ) fluctuations—naturally since we explicitly optimized this—and the diff contributions are small compared to the IPL ones. In the $N p T$ case, on the other hand, the diff contributions to the fluctuations of $U$ are more than double the IPL ones; this is not the case for the diff contributions to $W$⁠, though they are still a larger fraction of the total than in the $N V T$ case. These numbers are consistent with Fig. 5. The diff contributions to the energy must be larger than the IPL ones because the latter is negatively correlated with the true energy. The fact that the variance of $U diff$ is smaller than the sum of those of $U lin$ and $U rest$⁠, in the $N V T$ case, indicates that the latter two are negatively correlated. This is presumably due to bond lengths around $r c$ which alternately are counted as part of $U lin$ and as part of $U rest$ as they fluctuate. This effect is less noticeable in the $N p T$ case; there we see clearly that fluctuations in $U lin$ account for most of those in $U rest$⁠.

Based on the above, we can now answer the question: Is it possible to predict whether or not a liquid is strongly correlating by inspection of its potential (i.e., without simulating virial and potential energy fluctuations)? For liquids with particles interacting by pair potentials, the answer is in the affirmative: the liquid is strongly correlating if the potential around the first peak of the structure factor (the typical interparticle distance) may be fitted well by an eIPL and if the potential involves only one characteristic length. For more general potentials, the situation is more complex. Thus, it is possible to construct many-body potentials with angular dependencies which scale like IPL pair potentials; these have 100% $W U$ correlation because this property follows whenever the potential is an Euler homogeneous function. In most realistic cases, however, systems with angular dependencies are not expected to be strongly correlating. An example is the coarse-grained model of water using a short-range many-body potential recently introduced by Molinero and Moore,²⁷ a model that reproduces water’s properties with surprising accuracy. This potential is not strongly correlating because close to water’s density maximum the virial/potential energy correlation coefficient $R$ must be close to zero (Paper I). Examples of potentials with two length scales, that for this reason are not strongly correlating, are the Jagla potential²⁸ and the Dzugutov potential.⁷ Likewise, the addition of Coulomb terms to a LJ-type potential generally ruins strong correlations (Paper I).

According to the eIPL explanation detailed in Sec. III, strong $W U$ correlations characterize all configurations of the LJ liquid at a given volume. This means that the correlations should be there also under nonequilibrium conditions if the volume is kept constant. In this section we present numerical evidence that this prediction is indeed fulfilled, even for molecular liquids, provided that they are strongly correlating in their equilibrium $W U$ fluctuations.

Figure 6(a) shows the results for a temperature down-jump at constant volume, starting and ending in equilibrium (⁠$N V T$ simulations).²⁹ The system studied is an asymmetric dumbbell liquid consisting of two different-sized LJ particles glued together by a bond of fixed length, with parameters chosen to mimic toluene.² The system was first equilibrated at $300 K$⁠. The green ellipse is the scatter plot of simultaneous values of $U$ and $W$ in equilibrium at $T = 300 K$⁠. The strong $W U$ correlation is revealed by the elongation of the ellipse (⁠$R = 0.97$⁠; $γ = 6.1$⁠). When the liquid is similarly equilibrated at $240 K$⁠, the red blob appears. To test for correlation in an out-of-equilibrium situation, we changed temperature abruptly from the $300 K$ equilibrium situation to $240 K$⁠. The blue points show how virial and potential energy evolve following the temperature down-jump. Clearly, strong $W U$ correlations are also present during the aging toward equilibrium. Figure 6(b) shows the same phenomenon for the Lewis–Wahnström ortho-terphenyl (LW OTP) model, which consists of three LJ spheres at fixed length and angle with parameters optimized to mimic ortho-terphenyl.³⁰ This liquid is also strongly correlating (⁠$R = 0.91$⁠; $γ = 7.6$⁠). The colors are as in Fig. 6(a): green gives a high-temperature equilibrium state $( T = 600 K )$⁠, red gives a low-temperature equilibrium state $( T = 450 K )$⁠, and the blue points show the aging toward equilibrium after changing temperature from $600 to 450 K$⁠. The picture is the same as in Fig. 6(a). The blue points follow the dashed line. Thus, virial and potential energy correlate strongly also for far-from-equilibrium states. Figure 6(c) plots $W ( t )$ and $U ( t )$ after the temperature jump for the data in Fig. 6(a) for the asymmetric dumbbell liquid. $W ( t )$ and $U ( t )$ follow each other closely on the picosecond time scale as well as in their slow, overall drift to equilibrium.

What happens when the same simulation scheme is applied to a liquid that is not strongly correlating? An example is SPC water, where the hydrogen bonds are mimicked by Coulomb interactions.³¹ Figure 7 shows results of simulations of SPC water at two different densities, corresponding to (a) low pressure and (b) very high pressure . In the first case, virial and potential energy are virtually uncorrelated (⁠$R = 0.05$⁠, $T = 260 K$⁠); in the second case, correlations are somewhat stronger (⁠$R = 0.34$⁠, $T = 260 K$⁠), though still weak. As in Fig. 6, green denotes the initial, high-temperature equilibrium, red denotes the low-temperature equilibrium, and blue denotes the aging toward equilibrium. Clearly, for this system, $W$ and $U$ are not closely linked to one another during the relaxation toward equilibrium.

A different far-out-of-equilibrium situation is that of crystallization of a supercooled liquid monitored at fixed volume and temperature. To the best of our knowledge, crystallization of the LW OTP model has not been reported before, but Fig. 8 shows that for simulations over microseconds the supercooled liquid crystallizes at $T = 375 K$ and $ρ = 1.135 g ∕ cm 3$⁠. In the crystal, the LJ spheres are arranged in a slightly deformed bcc lattice where the molecules have otherwise random orientation. The crystal is shown in the inset of Fig. 8(a). Figure 8(a) shows how time-averaged pressure and energy develop during crystallization. Contributions to pressure and energy from the momentum degrees of freedom are virtually constant after averaging over $1 ns$⁠, so strong $W U$ correlations manifest themselves in strong averaged pressure/averaged energy correlations. Clearly, averaged pressure and averaged energy follow each other closely before, during, and after the crystallization. This confirms the above finding, as well as those of Paper I, that strong correlations apply also for the crystalline phase of a strongly correlating liquid. Note that the slope $γ$⁠, the proportionality constant between virial and potential energy fluctuations, is virtually unaffected by the crystallization. The persistence of strong virial/potential energy correlations during crystallization—and the insignificant change in $γ$—is noticeable because physical characteristics are rarely unaffected by a first-order phase transition. These simulations show that the property of strong virial/potential energy correlations pertains to the intermolecular potential, not to the particular microscopic configurations under study. Figure 8(b) shows the radial distribution functions for the liquid and crystalline phases.

The above out-of-equilibrium simulations show that the property of strong $W U$ correlation is not confined to thermal equilibrium. This is consistent with the eIPL description given in the previous section, but was shown here to apply even for molecular models. It appears that strongly correlating liquids have a particularly simple potential energy surface. These results have significance for any out-of-equilibrium situation. Consider the potential energy landscape picture of viscous liquid dynamics^32–35 according to which each configuration has an underlying inherent state defined via a deepest-descent quench, a state that contains most information relevant to the slow dynamics, which may be regarded as jumps between different inherent states.^32–35 Figure 9 shows a $W U$ plot of the asymmetric dumbbell model in different situations at same density: equilibrium states (upper right) and their corresponding inherent states (lower left, one quench per temperature), and glasses at different temperatures in between. A glass is an out-of-equilibrium state, of course, and inherent states may be regarded as zero-temperature glasses. The plot shows, once again, that strong correlations are present also far from equilibrium.

Most of the simulations of Paper I were done in the $N V T$ ensemble. An obvious question is how the correlations differ between ensembles. This section develops the necessary theory needed to answer this question and compares first the $N V T$ and $N V E$ ensembles, then the $N V T$ and $N p T$ ensembles.

It is well known that although simple thermodynamic averages are independent of ensemble,³⁶ fluctuations are generally ensemble dependent. Consider two ensembles, one with extensive variable $F$ held fixed and one with its conjugate intensive variable $f$ held fixed (defined so their product is dimensionless). The other parameters defining the ensembles are the same. The covariance of observables $A$ and $B$ in the two ensembles are then related as^6,37

To compare the $N V T$ and $N V E$ ensembles, we take $F$ and $f$ as the energy $E$ and the inverse temperature $β = 1 ∕ ( k B T )$⁠, respectively, keeping the volume fixed. Noting that $∂ ∕ ∂ β = − k B T 2 ∂ ∕ ∂ T$ and $( ∂ β ∕ ∂ E ) V = 1 ∕ ( − k B T 2 C V )$⁠, where $C V$ is the extensive isochoric specific heat, $C V = V c V$⁠, the covariance of $U$ and $W$ is given by

Here, we have introduced the excess parts of the isochoric specific heat and pressure coefficient, $C V ex$ and $β V ex$⁠, respectively (note the change in notation from Paper II where the superscript “conf” was used—“ex” is, however, more standard in liquid state theory). These two quantities are given by (for monatomic liquids)

For the variances, one has

The above implies that the $N V E$$W U$-correlation coefficient $R E$ (in the following subscripts $E$ or $T$ indicate the $N V E$ or $N V T$ ensemble, respectively) is given by

We wish to express the right side in terms of the $N V T$ coefficient $R T$⁠. The definition of $R T$ implies that

Inserting this into Eq. (44) and making use of the fluctuation relations $⟨ Δ U Δ W ⟩ N V T = k B T 2 V β V ex$ and $⟨ ( Δ U ) 2 ⟩ N V T = k B T 2 C V ex$ (see, e.g., Appendix B of Paper I) gives

After squaring and cancelling factors of $k B T 2$ and $V β V ex$⁠, we get an expression relating $R E$ to $R T$⁠,

To get a feel for the relation, divide by $R T 2$⁠. This yields

First one notes that when $R T = 1$⁠, the denominator on the right side becomes equal to the numerator and $R E = 1$⁠. That is, the property of perfect correlation is independent of (fixed volume) ensemble. When $R T < 1$⁠, the denominator becomes greater than the numerator, and so $R E 2 < R T 2$⁠. That is, the correlation coefficient is smaller in the $N V E$ ensemble than in the $N V T$ one. How can we understand this? Consider the set of $W U$ points sampled by the system during an $N V E$ trajectory; this is an elongated blob in the $W U$ diagram. Changing the energy will cause the blob to move along a line almost parallel with the long axis of the blob (Paper I). Switching to the $N V T$ ensemble is equivalent to superposing several of these collinear blobs on top of each other—the result is necessarily longer, but not wider. This corresponds to higher correlation.

Simulation data confirming relation (47) are presented in Table II.

In the $N p T$ ensemble, where volume is allowed to fluctuate, we must consider different variables. The natural variables to correlate are the excess enthalpy $H ex ≡ U + p V$ and the volume $V$⁠. We use again Eq. (37), but now take $F$ as $V$ and $f$ as $p β$⁠, the pressure times inverse temperature, keeping temperature constant. Equation (37) becomes

The details of the calculation, which are somewhat tedious, are given in the Appendix. The result for the $H ex V$ correlation coefficient is

where

and

Notice that $R H ex V , N p T$ is strictly less than unity—even for perfectly correlating liquids (that is, with perfect $W U$ correlations in the $N V T$ and $N V E$ ensembles). For the LJ simulation of Fig. 1(b), $R H ex V , N p T = 0.86$ (the $N V T$$W U$ correlation coefficient for the same state point is 0.94). Unlike the situation when comparing the $N V T$ and $N V E$ ensembles, there does not seem to be a simple relation between the two correlation coefficients. It seems likely, though, that the $H ex V$ correlation in the $N p T$ ensemble is generally smaller than the $W U$ correlation in the $N V T$ ensemble; thus, the property of strong correlation is less evident in the $N p T$ ensemble.

The property of strong virial/potential energy correlation does not just refer to microscopic properties that are only accessible in simulation, it also has consequences for the liquid’s thermodynamics. The first subsection below relates the slope [Eq. (6)] to the Grüneisen parameter, the second subsection shows how to give a general thermodynamic formulation of the property of strong correlations.

The Grüneisen parameter was originally introduced to characterize the volume dependence of normal modes of a crystal,^38,39

where $ω i$ is the frequency of the $i th$ normal mode. By assuming that $γ i$ is the same for all modes and denoting the common value by $γ G$⁠, one can derive the Mie–Grüneisen equation of state³⁸ as

Here, $p$ is the pressure, $u ( v )$ with $v = V ∕ N$ is the “static” energy of the crystal per atom (the energy of the force-free configuration about which vibrational motion occurs), and $E vib$ is the vibrational energy. In general $γ G$ depends on volume, but this dependence is typically small enough that it can be neglected. From Eq. (54) it follows that, if $E$ is the total, thermally averaged internal energy, one has $( ∂ p ∕ ∂ T ) V = ( γ G ∕ V ) ( ∂ E vib ∕ ∂ T ) V = ( γ G ∕ V ) ( ∂ E ∕ ∂ T ) V$⁠, i.e.,