In this third paper of the series, which started with Bailey et al [J. Chem. Phys. 129, 184507 (2008); ibid. 129, 184508 (2008)], we continue the development of the theoretical understanding of strongly correlating liquids—those whose instantaneous potential energy and virial are more than 90% correlated in their thermal equilibrium fluctuations at constant volume. The existence of such liquids was detailed in previous work, which identified them, based on computer simulations, as a large class of liquids, including van der Waals liquids but not, e.g., hydrogenbonded liquids. We here discuss the following: (1) the scaling properties of inverse powerlaw and extended inverse powerlaw potentials (the latter includes a linear term that “hides” the approximate scale invariance); (2) results from computer simulations of molecular models concerning outofequilibrium conditions; (3) ensemble dependence of the virial/potentialenergy correlation coefficient; (4) connection to the Grüneisen parameter; and (5) interpretation of strong correlations in terms of the energybond formalism.
I. INTRODUCTION
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 $\Delta $ 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” $\gamma $, which we here define^{1,4} as the ratio of standard deviations,
In the limit of perfect correlation $ ( R \u2192 1 ) $, $\gamma $ becomes equal to the slope of average $W$ as a function of average $U$ at fixed volume, $\gamma = ( \u2202 W / \u2202 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 LennardJones (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 manybody 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 powerlaw (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 singlecomponent LJ liquid in detail. Building on the fact that IPL potentials $ \u221d r \u2212 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 \u223c 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^{1} 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 lowtemperature/lowpressure (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 singleparticle 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^{8} in this series goes on to consider varyingdensity 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 $\gamma $ 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 $\gamma $, namely the one defined by Eq. (6) (in Paper IV the definition of $\gamma $ 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 inversepower law” (eIPL) potential introduced in Paper II, which includes the abovementioned 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 volumedependent 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 outofequilibrium 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 pressureenergy correlations to the thermodynamic Grüneisen parameter $ \gamma G $, showing that the slope $\gamma $ [Eq. (6)] is larger than $ \gamma G $ by roughly a factor involving the ratio of excess (configurational) to total constantvolume specific heat. This ratio is around 2 for many simple liquids.^{9} The second part formulates the property of strong correlation in the energybond language known as “network thermodynamics.”^{10}
II. PROPERTIES OF IPL SYSTEMS AND GENERALIZATIONS
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.
A. IPL potentials
IPL potentials—sometimes referred to as softsphere 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 lowpressure 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 \u2212 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 \u2212 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 ( \rho \Lambda 3 / e ) $, where $ \rho = N \u2215 V $ is the particle number density and $ \Lambda = h \u2215 2 \pi 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 welldefined extensive thermodynamic limit $ ( N \u2192 \u221e ) $.^{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 $\rho $ to the power $ n \u2215 3 $ over temperature $ T , \varphi ( \rho n \u2215 3 \u2215 T ) $ (Klein’s theorem^{11,12}):
Equation (8) implies that a number of derived quantities are also functions of $ \rho n \u2215 3 \u2215 T $. As important examples, recall the following standard identities: the excess entropy: $ S ex = \u2212 ( \u2202 F ex \u2215 \u2202 T ) V $, the potential energy: $ U = F ex + T S ex $, the virial $ W = \u2212 V ( \u2202 F ex \u2215 \u2202 V ) T $, the excess isothermal bulk modulus: $ K T ex = V ( \u2202 2 F ex \u2215 \u2202 V 2 ) T $, the excess isochoric specific heat per unit volume: $ c V ex = \u2212 ( T \u2215 V ) ( \u2202 2 F ex \u2215 \u2202 T 2 ) V $, the excess pressure coefficient: $ \beta V ex = ( 1 \u2215 V ) ( \u2202 W \u2215 \u2202 T ) V = \u2212 \u2202 2 F ex \u2215 \u2202 T \u2202 V $. These quantities are all functions of the single variable $ \rho n \u2215 3 \u2215 T $; more accurately, one has [where $ f 1 ( x ) = x \varphi \u2032 ( x ) \u2212 \varphi ( x ) $, $ f 2 ( x ) = x \varphi \u2032 ( x ) $, $ f 3 ( x ) = ( n \u2215 3 ) 2 x 2 \varphi \u2033 ( x ) + [ ( n \u2215 3 ) + ( n \u2215 3 ) 2 ] x \varphi \u2032 ( x ) $, and $ f 4 ( x ) = \u2212 x 2 \varphi \u2033 ( x ) $],
The functions $\varphi $ and $ f 1 , \u2026 , 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 $ \rho n \u2215 3 \u2215 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 , \u2026 , N ) $ is a solution to Newton’s equations. Then, it is straightforward to show that $ r i ( 1 ) ( t ) = \alpha r i ( \lambda t ) $ is also a solution if $ \alpha \u2212 ( n + 2 ) = \lambda 2 $. In particular, if $ r i ( t ) $ refers to equilibrium ($ N V E $ or $ N V T $) dynamics at a state point with density $ \rho 0 $ and temperature $ T 0 $, then $ r i ( 1 ) ( t ) = \alpha r i ( \lambda t ) $ refers to equilibrium dynamics at density $ \rho 1 = \rho 0 \u2215 \alpha 3 $ and temperature $ T 1 = T 0 \alpha 2 \lambda 2 $ (temperature scales as the meansquare velocity and velocities get a factor $ \alpha \lambda $). Using the above relation between $\alpha $ and $\lambda $, this implies that
This means that two states with different densities and temperatures but the same $ \rho n \u2215 3 \u2215 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 $ \tau A $ via $ \tau A = \u222b 0 \u221e \u27e8 A ( 0 ) A ( t ) \u27e9 d t \u2215 \u27e8 A 2 \u27e9 $, it follows that any two states with same $ \rho n \u2215 3 \u2215 T $ have the same “reduced” (dimensionless) relaxation time $ \tau \u0303 A $, if this quantity is defined by $ \tau \u0303 A = \tau A \u2215 t 0 $ where the characteristic time $ t 0 $ is defined by $ t 0 = \rho \u2212 1 \u2215 3 m \u2215 k B T $. Similarly, if one defines the reduced diffusion constant $ D \u0303 $ by $ D \u0303 = D \u2215 D 0 $, where $ D 0 = \rho \u2212 2 \u2215 3 \u2215 t 0 = \rho \u2212 1 \u2215 3 k B T \u2215 m $, then $ D \u0303 $ 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 $ \rho n \u2215 3 \u2215 T $ are identical provided that lengths are scaled by $ \rho \u2212 1 \u2215 3 $.
B. Inheritance of scaling properties by generalized IPL potentials
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 , \u2026 , 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 , \u2026 , r N ) \u2245 f ( V ) $—at least for states that carry Boltzmann weights of any significance. The approximate identity $ U diff ( r 1 , \u2026 , 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 $ \rho n \u2215 3 \u2215 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 ) \u2260 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 $ \rho n \u2215 3 \u2215 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.
III. LJ AS A GENERALIZED IPL POTENTIAL: THE eIPL POTENTIAL APPROXIMATION
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 singlepair 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 firstpeak 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 \u2215 6 \sigma $, 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 nearzero pressure, the best fit was obtained with $ n = 19.2 $ (while the exponent implied by the slope [Eq. (6)], $ \gamma = 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 firstpeak 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 powerlaw potential almost without affecting the fluctuations as long as the volume is held constant. The analysis of Paper II, which also included an indepth 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 \gamma $, where $\gamma $ 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 powerlaw 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 nearzero pressure state point used in Fig. 1 of Paper I, the exponent determined from $\gamma $ is $ n = 3 \gamma = 18.9 $ and the optimal value of $A$ is $ 1.3437 \u03f5 \sigma 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 \u2212 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 ) \u223c 10 \u2212 5 \u03f5 $, the cutoff is given by $ r c = \u2212 B \u2215 C $. For the fit shown in Fig. 3(b), $ r c = 1.61 \sigma $.
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 bondlength 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 singlecomponent system has a rather welldefined nearestneighbor distance, which, at densities where the pressure is not too negative, will be less than $ r c $. In a threedimensional liquid, however, the nearestneighbor 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 threedimensional (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 \u2192 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 \u2aa2 1 $, while the linear term gets multiplied by $ \u2212 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 $.
Ensemble .  Quantity .  LJ .  IPL .  diff .  lin .  rest . 

$ N V T $  $U$  0.0231  0.0225  0.0075  0.0085  0.0063 
$W$  0.1468  0.1402  0.0227  0.0301  0.0350  
$ N p T $  $U$  0.0484  0.0320  0.0665  0.0539  0.0144 
$W$  0.1704  0.1997  0.0589  0.0417  0.0430 
Ensemble .  Quantity .  LJ .  IPL .  diff .  lin .  rest . 

$ N V T $  $U$  0.0231  0.0225  0.0075  0.0085  0.0063 
$W$  0.1468  0.1402  0.0227  0.0301  0.0350  
$ N p T $  $U$  0.0484  0.0320  0.0665  0.0539  0.0144 
$W$  0.1704  0.1997  0.0589  0.0417  0.0430 
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 manybody 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 coarsegrained model of water using a shortrange manybody potential recently introduced by Molinero and Moore,^{27} 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^{28} and the Dzugutov potential.^{7} Likewise, the addition of Coulomb terms to a LJtype potential generally ruins strong correlations (Paper I).
IV. OUTOFEQUILIBRIUM DYNAMICS IN MOLECULAR MODELS
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.
A. Temperature downjump simulations of three molecular model liquids
Figure 6(a) shows the results for a temperature downjump at constant volume, starting and ending in equilibrium ($ N V T $ simulations).^{29} The system studied is an asymmetric dumbbell liquid consisting of two differentsized LJ particles glued together by a bond of fixed length, with parameters chosen to mimic toluene.^{2} 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 $; $ \gamma = 6.1 $). When the liquid is similarly equilibrated at $ 240 K $, the red blob appears. To test for correlation in an outofequilibrium 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 downjump. 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 orthoterphenyl (LW OTP) model, which consists of three LJ spheres at fixed length and angle with parameters optimized to mimic orthoterphenyl.^{30} This liquid is also strongly correlating ($ R = 0.91 $; $ \gamma = 7.6 $). The colors are as in Fig. 6(a): green gives a hightemperature equilibrium state $ ( T = 600 K ) $, red gives a lowtemperature 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 farfromequilibrium 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.^{31} 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, hightemperature equilibrium, red denotes the lowtemperature 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.
B. Pressure and energy monitored during crystallization of a supercooled liquid: The LW OTP model
A different faroutofequilibrium 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 $ \rho = 1.135 g \u2215 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 timeaveraged 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 $\gamma $, 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 $\gamma $—is noticeable because physical characteristics are rarely unaffected by a firstorder 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.
C. Glasses and inherent states
The above outofequilibrium 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 outofequilibrium 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 deepestdescent 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 outofequilibrium state, of course, and inherent states may be regarded as zerotemperature glasses. The plot shows, once again, that strong correlations are present also far from equilibrium.
V. ENSEMBLE DEPENDENCE OF THE CORRELATION COEFFICIENT
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,^{36} 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}
A. $ N V T $ versus $ N V E $
To compare the $ N V T $ and $ N V E $ ensembles, we take $F$ and $f$ as the energy $E$ and the inverse temperature $ \beta = 1 \u2215 ( k B T ) $, respectively, keeping the volume fixed. Noting that $ \u2202 \u2215 \u2202 \beta = \u2212 k B T 2 \u2202 \u2215 \u2202 T $ and $ ( \u2202 \beta \u2215 \u2202 E ) V = 1 \u2215 ( \u2212 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 $ \beta 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 $ \u27e8 \Delta U \Delta W \u27e9 N V T = k B T 2 V \beta V ex $ and $ \u27e8 ( \Delta U ) 2 \u27e9 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 \beta 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.
System .  $\rho $ .  $T$ .  $ R T $ .  $ C V ex \u2215N$ .  $ R E $ [Eq. (47)] .  $ R E $ (measured) . 

LJ  1.00  1.00  0.991  1.5  0.982  0.983 
LJ  1.00  0.80  0.991  1.7  0.981  0.981 
LJ  0.82  0.80  0.943  0.90  0.912  0.918 
LJ  0.82  0.67  0.949  1.3  0.909  0.904 
KABLJ  1.2  0.47  0.936  2.1  0.859  0.862 
System .  $\rho $ .  $T$ .  $ R T $ .  $ C V ex \u2215N$ .  $ R E $ [Eq. (47)] .  $ R E $ (measured) . 

LJ  1.00  1.00  0.991  1.5  0.982  0.983 
LJ  1.00  0.80  0.991  1.7  0.981  0.981 
LJ  0.82  0.80  0.943  0.90  0.912  0.918 
LJ  0.82  0.67  0.949  1.3  0.909  0.904 
KABLJ  1.2  0.47  0.936  2.1  0.859  0.862 
B. $ N V T $ versus $ N p T $
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 \u2261 U + p V $ and the volume $V$. We use again Eq. (37), but now take $F$ as $V$ and $f$ as $ p \beta $, 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.
VI. THERMODYNAMICS OF STRONGLY CORRELATING LIQUIDS
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.
A. Relation to the Grüneisen parameter
The Grüneisen parameter was originally introduced to characterize the volume dependence of normal modes of a crystal,^{38,39}
where $ \omega i $ is the frequency of the $ i th $ normal mode. By assuming that $ \gamma i $ is the same for all modes and denoting the common value by $ \gamma G $, one can derive the Mie–Grüneisen equation of state^{38} as
Here, $p$ is the pressure, $ u ( v ) $ with $ v = V \u2215 N $ is the “static” energy of the crystal per atom (the energy of the forcefree configuration about which vibrational motion occurs), and $ E vib $ is the vibrational energy. In general $ \gamma 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 $ ( \u2202 p \u2215 \u2202 T ) V = ( \gamma G \u2215 V ) ( \u2202 E vib \u2215 \u2202 T ) V = ( \gamma G \u2215 V ) ( \u2202 E \u2215 \u2202 T ) V $, i.e.,
This expression is the slope of the pressure versus energy curve at fixed volume, analogous to the $\gamma $ of Eq. (6) but for the presence of the kinetic terms (recall that for an IPL liquid $ \gamma = n \u2215 3 = ( \u2202 W \u2215 \u2202 U ) V $). If $ \alpha p $ is the coefficient of thermal expansion, $ K T $ is the isothermal bulk modulus, and $ c V = C V \u2215 V $ is the isochoric specific heat per unit volume, Eq. (55) implies via standard thermodynamic identities
This relation allows $ \gamma G $ to be determined from experimentally accessible quantities; in fact, Eq. (56) can be taken as a thermodynamic definition of $ \gamma G $.^{40}
There have been suggestions of how to connect the socalled density scaling exponent^{24,25}—the one controlling the relaxation time via the variable $ \rho \gamma \u2215 T $—with the Grüneisen parameter, notably by Roland and coworkers.^{41–43} In Ref. 41, equality of $ \gamma G $ and $\gamma $ was argued theoretically by reference to the Avramov model. More recently, Roland and Casalini showed^{43} that equality is not consistent with experimental results; rather $ \gamma G $ is smaller than $\gamma $ by a factor of order 3. This discrepancy was reconciled in the context of the entropy model for relaxation (by which the relaxation time is a unique function of the socalled configurational entropy $ S c $), by introducing an excess or corrected Grüneisen parameter defined via
Here, $ \Delta c V $ is the difference in the isochoric specific heats per unit volume between the liquid and the glass. Because $ \Delta c V $ is smaller than $ c V $, one has $ \gamma G corr > \gamma G $. By arguing from the experimental data that the nonconfigurational part of the entropy (associated with vibrations, equal to the entropy of the glass) is independent of volume and assuming that $ \Delta C V $ is constant, Roland and Casalini derive
Hence, $ \gamma G corr $ is the densityscaling exponent.
We take a different approach to connecting the Grüneisen parameter with the slope $\gamma $ (which provides a good estimate of the density scaling exponent, see Ref. 22 and Paper IV). Instead of splitting the entropy, we split the pressure into potential and kinetic parts and get from $ \gamma G = V ( \u2202 p \u2215 \u2202 E ) V $
Expressing the temperature derivatives in terms of fluctuations (Appendix B of Paper I) gives
In the limit of strong correlation, one can replace $ \u27e8 \Delta U \Delta W \u27e9 $ with $ \gamma \u27e8 ( \Delta U ) 2 \u27e9 $. Writing the resulting expression in terms of the excess (configurational) specific heat $ C V ex $ gives
In the harmonic approximation, good for many simple liquids close to their melting point,^{9} $ C V ex = ( 3 \u2215 2 ) N k B $ (while it is generally larger in the supercooled state). Thus, the term $ N k B $ in the numerator is expected to be roughly a factor of 10 smaller than the other term; we drop it and arrive at
This ratio is around onehalf in the harmonic approximation, otherwise larger.
B. Energybond formulation of the strongly correlation property
This section derives a general thermodynamic condition of the property of strong $ W U $ correlations, a condition which linearly constrains small variations in entropy, volume, temperature, and pressure [Eq. (72), below]. It is convenient to approach the problem from a general point of view. The energybond formalism provides an abstract description of the interactions between a system and its surroundings.^{10,44–49} An energy bond has an “effort” variable $ e ( t ) $ and a “flow” variable $ f ( t ) $, where $ e ( t ) f ( t ) $ is the free energy transferred into the system per unit time. The “displacement” $ q ( t ) $ is the timeintegrated flow, i.e., $ q \u0307 ( t ) = f ( t ) $. The energybond formalism is general, but we only discuss the linear case where it is most useful. Thus, we consider a system that is slightly perturbed from equilibrium. It is assumed that the underlying microscopic dynamics is described by a stochastic equation, i.e., inertial forces are ignored. This is believed to be a good approximation whenever the highly viscous phase is approached; thus the below results are mainly of relevance to viscous liquids.
Linearresponse theory is characterized by the fluctuationdissipation (FD) theorem, which in the energybond formalism is given as follows. Consider a situation with $n$ energy bonds and external control of the effort variables. If $ \u27e8 f i ( 0 ) f j ( t \u2032 ) \u27e9 0 $ is the equilibrium flow autocorrelation function, the average flow at time $t$ is given by
If the arbitrary additive constants of the displacements are chosen such that $ q i = 0 $ on average, the timeintegrated version of this is
If the flow variables are externally controlled, the FD theorem is
In most cases, efforts are invariant under time reversal and flows change sign. The Onsager reciprocity relation is $ \u27e8 f i ( 0 ) f j ( t ) \u27e9 0 = \u27e8 f j ( 0 ) f i ( t ) \u27e9 0 $ (or $ \u27e8 e i ( 0 ) e j ( t ) \u27e9 0 = \u27e8 e j ( 0 ) e i ( t ) \u27e9 0 $, depending on which variables are externally controlled and which are free to fluctuate). From the FD theorem, expressions for the frequencydependent response functions are easily derived. Consider, for instance, the compliances $ J i j ( \omega ) $, defined by $ J i j ( \omega ) = q i ( \omega ) \u2215 e j ( \omega ) $ for a periodic situation with infinitesimal perturbations around equilibrium, $ e ( t ) = Re [ e ( \omega ) exp ( i \omega t ) ] $, etc. For these quantities, the FD theorem implies
The case relevant to strongly correlating liquids is that of two energy bonds that are not independent, as we now proceed to show. The two energy bonds are those of standard thermodynamics, reflecting the fundamental relation $ d E = T d S \u2212 p d V $.^{50,51} These two energy bonds are given as follows. The thermal energy bond has temperature variation as the effort variable and entropy variation as the displacement variable ($ e 1 ( t ) = \delta T ( t ) $, $ q 1 ( t ) = \delta S ( t ) $), the mechanical energy bond has pressure variation as the effort variable and the negative volume variation as the displacement variable ($ e 2 ( t ) = \delta p ( t ) $, $ q 2 ( t ) = \u2212 \delta V ( t ) $). Usually, the two standard thermodynamic energy bonds are independent, but we are here interested in the case when they are not.
Treating the problem of two constrained energy bonds from a general perspective, we shall prove that the following four criteria are equivalent.
 The variables of the two energy bonds are linearly constrained as follows:(67)$ a q 1 ( t ) + b q 2 ( t ) = c e 1 ( t ) + d e 2 ( t ) . $

The system’s relaxing properties, i.e., its noninstantaneous responses, are described^{50} by a single variable $ \u03f5 ( t ) $ as follow:
(68)$ q 1 ( t ) = J 11 \u221e e 1 ( t ) + J 12 \u221e e 2 ( t ) + \gamma 1 \u03f5 ( t ) , $$ q 2 ( t ) = J 21 \u221e e 1 ( t ) + J 22 \u221e e 2 ( t ) + \gamma 2 \u03f5 ( t ) . $In these equations the $ J \u221e $’s are the compliances referring to the shorttime, nonrelaxing response (the highfrequency response). Note that $ J 12 \u221e = J 21 \u221e $ by the FD theorem.
 The relaxing parts of the correlation functions entering into Eq. (66) are proportional. More precisely, the correlation functions obey(69)$ \u27e8 q 1 ( 0 ) f 1 ( t ) \u27e9 0 \u221d \u27e8 q 1 ( 0 ) f 2 ( t ) \u27e9 0 \u221d \u27e8 q 2 ( 0 ) f 1 ( t ) \u27e9 0 \u221d \u27e8 q 2 ( 0 ) f 2 ( t ) \u27e9 0 ( t \u2260 0 ) $and(70)$ \u27e8 q 1 ( 0 ) f 1 ( t ) \u27e9 0 \u27e8 q 2 ( 0 ) f 2 ( t ) \u27e9 0 = \u27e8 q 1 ( 0 ) f 2 ( t ) \u27e9 0 \u27e8 q 2 ( 0 ) f 1 ( t ) \u27e9 0 ( t \u2260 0 ) . $
Note that by differentiation, Eq. (69) implies for $ t \u2260 0 $ that $ \u27e8 f 1 ( 0 ) f 1 ( t ) \u27e9 0 \u221d \u27e8 f 1 ( 0 ) f 2 ( t ) \u27e9 0 \u221d \u27e8 f 2 ( 0 ) f 1 ( t ) \u27e9 0 \u221d \u27e8 f 2 ( 0 ) f 2 ( t ) \u27e9 0 $.

The dynamic Prigogine–Defay ratio^{50} $ \Lambda ( \omega ) $ is unity at all frequencies (where double prime denotes the imaginary part),
Proof that $ 1 \u21d4 2 $. By elimination of the variable $\u03f5$ from Eq. (68), 2 implies 1. To prove the reverse implication, suppose that Eq. (67) applies and fix the dimensions such that the constants $c$ and $d$ are dimensionless. Define $ J 11 \u221e = ( 1 + c ) \u2215 a $, $ J 12 \u221e = J 21 \u221e = \u2212 1 \u2215 b $, and $ J 22 \u221e = ( d b + a ) \u2215 b 2 $. Introducing the variables $ \u03f5 1 = q 1 \u2212 J 11 \u221e e 1 \u2212 J 12 \u221e e 2 $ and $ \u03f5 2 = q 2 \u2212 J 21 \u221e e 1 \u2212 J 22 \u221e e 2 $, it follows that $ a \u03f5 1 + b \u03f5 2 = 0 $. This means that we are in the situation described by Eq. (68) with a common relaxing variable to the two energy bonds, $ \u03f5 ( t ) \u221d \u03f5 1 ( t ) \u221d \u03f5 2 ( t ) $, and symmetric shorttime compliances, $ J 12 \u221e = J 21 \u221e $.
Proof that $ 2 \u21d2 3 $. In terms of functional derivatives with respect to the efforts at an earlier time $ ( t \u2032 < t ) $ since $ \u03f5 ( t ) $ for small variations in the effort variables is linear in these, via the FD theorem timereversal invariance implies that $ \delta q 1 ( t ) \u2215 \delta e 2 ( t \u2032 ) = \delta q 2 ( t ) \u2215 \delta e 1 ( t \u2032 ) $. Thus, Eq. (68) implies that $ \delta \u03f5 ( t ) \u2215 \delta e 2 ( t \u2032 ) \u221d \delta \u03f5 ( t ) \u2215 \delta e 1 ( t \u2032 ) $. From this, Eqs. (69) and (70) now follow via the FD theorem and Eq. (68).
Proof that $ 3 \u21d2 4 $. According to the FD theorem, the compliance matrix imaginary parts are given by $ J i j \u2033 ( \omega ) = ( 1 \u2215 k B T ) \u222b 0 \u221e \u27e8 q i ( 0 ) f j ( t \u2032 ) \u27e9 0 sin ( \omega t \u2032 ) d t \u2032 $. In conjunction with Eqs. (69) and (70), this implies that $ \Lambda ( \omega ) = 1 $ at all frequencies.
Proof that $ 4 \u21d2 2 $. We refer to calculations of Ref. 50 which considered a system described by stochastic dynamics, i.e., with no inertial forces. Generalization of the arguments given there for the two standard thermodynamic energy bonds to the case of two arbitrary energy bonds proves the required implication.
This completes the proof of the equivalence of points (1)–(4). For the case where the two energy bonds are the two standard thermodynamic bonds, the constraint Eq. (67) translates into (changing here the sign of $b$)
How does this all relate to strong $ W U $ correlations in liquids? Via the equivalence of Eq. (72) to Eq. (68) and to unity dynamic Prigogine–Defay ratio [Eq. (71)], the results derived in Refs. 1–3 imply that Eq. (72) describes a 100% correlating liquid subjected to small perturbations from equilibrium. Generally, for any strongly correlating liquid Eq. (72) is obeyed with good accuracy. Thus, Eq. (72) gives the required thermodynamic formulation of the hidden scale invariance characterizing strongly correlating liquids.
Equation (72) implies that for strongly correlating liquids, the four thermodynamic variables, entropy, volume, temperature, and pressure, cannot vary independently. Referring to Eq. (68), it is clear that for certain simultaneous changes in the four thermodynamic variables, the relaxing part is left unchanged; this suggests that for such changes, the system is taken to a state where it is immediately in thermal equilibrium. This observation inspired the works leading to Paper IV where “isomorphs” are introduced. These are curves in the phase diagram along which several quantities are invariant, and along which jumps from equilibrium at one state point take the system to a new state that is instantaneously in thermal equilibrium.
Finally, we would like to draw attention to an analog of strongly correlating liquids. Consider a relaxing dielectric such as, e.g., a highly viscous dipolar liquid placed in a metal capacitor. This system’s interaction with its surroundings may be described by two energy bonds: one energy bond is defined by the capacitor charge (electronic plus induced) and the voltage across the capacitor, the other energy bond is the induced dielectric charge at the capacitor surface and a fictive electric field only coupling to the liquid’s dipoles. Because of Gauss’ law, these two energy bonds are not independent, but constrained by a linear displacementfield relation of the form Eq. (67). Thus, from the energybond formalism point of view, a strongly correlating liquid is analogous to the standard measuring cell used for probing $ \u03f5 ( \omega ) $ of dipolar viscous liquids, with the strong virial/potential energy correlations reflecting one of Maxwell’s four equations.
VII. CONCLUDING REMARKS
We have illuminated a number of features of strongly correlating liquids’ hidden scale invariance. The linear term in the eIPL potential, which hides the approximate scale invariance, contributes little to the thermal fluctuations at fixed volume; this is why strongly correlating liquids inherit a number of IPL properties. As shown in previous papers,^{1–5} the hidden scale invariance has important experimental consequences, including that of density scaling.^{22,57} The general physical picture is that van der Waals liquids and metallic liquids—because they are strongly correlating—are simpler than hydrogenbonding liquids, ionic liquids, and covalently bonded liquids, which are not strongly correlating.
Paper IV further investigates the consequences of a liquid being strongly correlating. This is done by defining socalled isomorphs in the liquid’s phase diagram and showing that a number of properties to a good approximation are invariant along an isomorph. The isomorph definition does not refer to $ W U $ correlations. Only strongly correlating liquids have isomorphs, however; this is because the existence of isomorphs is a direct consequence of the hidden scale invariance characterizing strongly correlating liquids.
ACKNOWLEDGMENTS
We thank Tage Christensen and Søren Toxværd for helpful input. The center for viscous liquid dynamics “Glass and Time” is sponsored by the Danish National Research Foundation (DNRF).
APPENDIX: CALCULATING THE $ H ex V $ CORRELATION COEFFICIENT IN THE $ N p T $ ENSEMBLE
Here, we provide the details of the calculation of the correlation coefficient between volume and excess enthalpy in the $ N p T $ ensemble. We apply Eq. (49) with $ A , B \u220a { U , V } $. First, we need the pressure derivatives at constant temperature of $ \u27e8 U \u27e9 $ and $ \u27e8 V \u27e9 $. Taking $U$ first, we have (noting that for simple averages like $ \u27e8 U \u27e9 $, it is not necessary to specify the ensemble because of equivalence of ensembles)
Note that $V$ in the derivative is without averaging signs since there it is a parameter of the relevant ensemble $ ( N V T ) $. The volume derivative of $ \u27e8 U \u27e9 $ is calculated as follows. The excess (configurational) partition function $ Z ( V , T ) $ is the integral (where $\beta =1/ k B T$)
Here, $\Gamma $ indexes points in configuration space and $ d \Gamma = d 3 N r \u2215 V N $. In the following, we use the configuration space identity $ \u2202 U \u2215 \u2202 V = \u2212 W \u2215 V $ valid for changes that scale the coordinates of the microscopic configurations, corresponding to fixed socalled reduced coordinates (compare Appendix A of Paper IV); constant temperature is implicit, as is the dependence of $U$ on $\Gamma $ and $V$,
Thus we have (adding the subscript $ N V T $ to the fluctuation expression since this is ensemble dependent)
The pressure dependence of $ \u27e8 V \u27e9 $ is given by
To keep the notation simple, averaging signs are henceforth omitted from simple averages such as $ \u27e8 V \u27e9 $, $ \u27e8 W \u27e9 $, etc. We can write expressions for the variances of $U$ and $V$ in the $ N p T $ ensemble using Eq. (49),
We need also the covariance
Now we have all we need to construct the $ H ex V $ correlation coefficient in the $ N p T $ ensemble where $ H ex =U+pV$. The covariance between $ H ex $ and $V$ is
where we have used $ p V = N k B T + W $. The variance of $ H ex $ is more tedious,
Again using $ p V = N k B T + W $ allows some simplication,
Now we can form the $ H ex V $ correlation coefficient,
where $ a = \u27e8 \Delta U \Delta W \u27e9 N V T + N ( k B T ) 2 $ and $ b 2 = K T V k B T \u27e8 ( \Delta U ) 2 \u27e9 N V T $.
REFERENCES
The system consisted of 512 asymmetric dumbbell molecules modeled as two LJ spheres connected by a rigid bond. The dumbbells were parametrized to mimic toluene. A large sphere (mimicking a phenyl group) was taken from the LewisWahnström OTP model (Ref. 30) with the parameters $mp=77.106u$, $\sigma p=0.4963nm$, and $\u03f5p=5.726kJ\u2215mol$. A small sphere (mimicking a methyl group) was taken from UAOPLS having $mm=15.035u$, $\sigma m=0.3910nm$, and $\u03f5m=0.66944kJ\u2215mol$. The bonds were kept rigid with a bond length of $d=0.29nm$. The volume was $V=77.27nm3$, giving an average pressure of approximately $1atm$. The temperature was held constant at $T=130K$ using the Nosé–Hoover thermostat. $NVT$ simulations were carried out using GROMACS software (Refs. 52 and 53) using the Nosé–Hoover thermostat (Refs. 54 and 55). Molecules were kept rigid using the LINCS algorithm (Ref. 56).
In particular, averages of quantities which are the sums of singleparticle functions (Ref. 6).