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., hydrogen-bonded liquids. We here discuss the following: (1) the scaling properties of inverse power-law and extended inverse power-law potentials (the latter includes a linear term that “hides” the approximate scale invariance); (2) results from computer simulations of molecular models concerning out-of-equilibrium conditions; (3) ensemble dependence of the virial/potential-energy correlation coefficient; (4) connection to the Grüneisen parameter; and (5) interpretation of strong correlations in terms of the energy-bond formalism.

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,

(1)
Here, K and U are the kinetic and potential energies, respectively, and T ( p 1 , , p N ) is the “kinetic temperature,” proportional to the kinetic energy per particle.6 The configurational contribution to pressure is the virial W , which is defined6 by

(2)

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

(3)
(4)

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 , defined1,4 by

(5)

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 define1,4 as the ratio of standard deviations,

(6)

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 correlations1 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 IV8 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.9 The second part formulates the property of strong correlation in the energy-bond language known as “network thermodynamics.”10 

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 written6,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 given6,21 by

(7)

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 theorem11,12):

(8)

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 ) ],

(9)
(10)
(11)
(12)
(13)
(14)

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

(15)

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,

(16)
(17)

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

(18)

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,”

(19)

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

(20)

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

(21)

From this follows directly that

(22)

This implies that

(23)
(24)
(25)
(26)
(27)
(28)

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).

FIG. 1.

(a) Scatter plot of total virial and potential energy (in LJ units) for the standard single-component LJ liquid at T = 80 K (argon units) and near-zero pressure, simulated at constant volume (density ρ = 34.6 mol l , argon units, left panel) and constant pressure ( 1.5 MPa , argon units, right panel). (b) Scatter plot of single-particle virial and potential energy for the same simulation as in the left panel of (a). The single-particle W U correlation is much weaker, R = 0.63 , showing that collective effects are crucial for the correlation.

FIG. 1.

(a) Scatter plot of total virial and potential energy (in LJ units) for the standard single-component LJ liquid at T = 80 K (argon units) and near-zero pressure, simulated at constant volume (density ρ = 34.6 mol l , argon units, left panel) and constant pressure ( 1.5 MPa , argon units, right panel). (b) Scatter plot of single-particle virial and potential energy for the same simulation as in the left panel of (a). The single-particle W U correlation is much weaker, R = 0.63 , showing that collective effects are crucial for the correlation.

Close modal

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

(29)

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

(30)

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

(31)

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

(32)

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

(33)

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

(34)
(35)

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,

(36)

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).

FIG. 2.

Comparison of g ( r ) for simulations using the LJ potential and two IPL potentials: the r 12 repulsive term in v LJ ( r ) and the IPL potential that optimizes the agreement in the fluctuations of potential energy and virial by minimizing Eq. (36). The left panel shows these at a density of 0.82 and a temperature of 0.67 (LJ units), the right one at a density of 0.90 and a temperature of 0.80 (where the r 12 potential leads to crystallization).

FIG. 2.

Comparison of g ( r ) for simulations using the LJ potential and two IPL potentials: the r 12 repulsive term in v LJ ( r ) and the IPL potential that optimizes the agreement in the fluctuations of potential energy and virial by minimizing Eq. (36). The left panel shows these at a density of 0.82 and a temperature of 0.67 (LJ units), the right one at a density of 0.90 and a temperature of 0.80 (where the r 12 potential leads to crystallization).

Close modal

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 σ .

FIG. 3.

(a) Illustration of the difference between the LJ potential v LJ ( r ) , the empirically matched IPL potential v IPL ( r ) with A = 1.3437 ϵ σ n and n = 18.9 , and their difference v diff ( r ) . (b) Linear fit, v lin ( r ) = min ( 0 , 3.6635 + 2.2756 r σ ) , to v diff between 0.95 σ and 1.5 σ , and the remainder v rest ( r ) v diff ( r ) v lin ( r ) (full black curve).

FIG. 3.

(a) Illustration of the difference between the LJ potential v LJ ( r ) , the empirically matched IPL potential v IPL ( r ) with A = 1.3437 ϵ σ n and n = 18.9 , and their difference v diff ( r ) . (b) Linear fit, v lin ( r ) = min ( 0 , 3.6635 + 2.2756 r σ ) , to v diff between 0.95 σ and 1.5 σ , and the remainder v rest ( r ) v diff ( r ) v lin ( r ) (full black curve).

Close modal

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].

FIG. 4.

Effect on fixed-volume fluctuations of adding a linear term to the IPL potential. The linear term is that shown in Fig. 3(b). Configurations were generated by an N V T simulation using the LJ potential, and the different determinations of energy (LJ, IPL, and eIPL) and virial were computed on these configurations. The dashed lines indicate a perfect match. Including the linear term when computing the energy improves the match to the true (LJ) fluctuations (the correlation coefficient goes from 0.950 to 0.970), while it reduces the match to the virial (the correlation coefficient goes from 0.987 to 0.971, which is probably related to the fact that the pair virial is discontinuous at r c —we find that smoothing the linear part around r c restores the match somewhat). The insets show the pair potentials and virials: brown dashed lines, LJ; black lines, IPL; and red lines, eIPL. The overall conclusion from Fig. 4 is that the addition of the linear term induces little change in the fluctuations.

FIG. 4.

Effect on fixed-volume fluctuations of adding a linear term to the IPL potential. The linear term is that shown in Fig. 3(b). Configurations were generated by an N V T simulation using the LJ potential, and the different determinations of energy (LJ, IPL, and eIPL) and virial were computed on these configurations. The dashed lines indicate a perfect match. Including the linear term when computing the energy improves the match to the true (LJ) fluctuations (the correlation coefficient goes from 0.950 to 0.970), while it reduces the match to the virial (the correlation coefficient goes from 0.987 to 0.971, which is probably related to the fact that the pair virial is discontinuous at r c —we find that smoothing the linear part around r c restores the match somewhat). The insets show the pair potentials and virials: brown dashed lines, LJ; black lines, IPL; and red lines, eIPL. The overall conclusion from Fig. 4 is that the addition of the linear term induces little change in the fluctuations.

Close modal

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).

FIG. 5.

Comparison of potential energy calculated using v LJ ( r ) , v IPL ( r ) (black points), and v eIPL ( r ) (red points) for configurations drawn from an N p T simulation using the LJ potential for the same state as in Fig. 4. The potential energy, in particular, is very poorly represented by the inverse power-law contribution when the volume is allowed to fluctuate. In fact, the correlation between U IPL and U LJ is not only weak, it is negative. Including the linear term makes a huge difference here, yielding a correlation coefficient of 0.977 between U eIPL and U LJ (changed from 0.201 ). The slope is somewhat less than unity, indicating that there are significant contributions from pair distances beyond r c (i.e., from the “rest” part of the potential). The linear term affects the virial fluctuations much less presumably because the derivative of the potential is dominated by the IPL term.

FIG. 5.

Comparison of potential energy calculated using v LJ ( r ) , v IPL ( r ) (black points), and v eIPL ( r ) (red points) for configurations drawn from an N p T simulation using the LJ potential for the same state as in Fig. 4. The potential energy, in particular, is very poorly represented by the inverse power-law contribution when the volume is allowed to fluctuate. In fact, the correlation between U IPL and U LJ is not only weak, it is negative. Including the linear term makes a huge difference here, yielding a correlation coefficient of 0.977 between U eIPL and U LJ (changed from 0.201 ). The slope is somewhat less than unity, indicating that there are significant contributions from pair distances beyond r c (i.e., from the “rest” part of the potential). The linear term affects the virial fluctuations much less presumably because the derivative of the potential is dominated by the IPL term.

Close modal

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 .

Table I.

Variances of potential energy U and virial W , and of various contributions to U and W , of two different ensembles at the LJ state point given by ρ = 0.82 and T = 0.67 (dimensionless units).

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 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,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 potential28 and the Dzugutov potential.7 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).29 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.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 ; γ = 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.30 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.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, 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.

FIG. 6.

Computer simulations of virial and potential energy during the aging of two strongly correlating molecular liquids following temperature down-jumps at constant volume ( N V T simulations). (a) The asymmetric dumbbell model at a density of ρ = 1.109 g ml . The liquid was first equilibrated at T = 300 K . Here, simultaneous values of virial and potential energy are plotted, producing the green ellipse, the elongation of which directly reflects the strong W U correlation. Temperature was then changed to T = 240 K where the red ellipse marks the equilibrium fluctuations. The aging process itself is given by the blue points. These points follow the line defined by the two equilibrium simulations, showing that virial and potential energy correlate also out of equilibrium. (b) Similar temperature down-jump simulation of the LW OTP system (Ref. 30). Again, green marks the high-temperature equilibrium ( T = 600 K ) , red the low-temperature equilibrium ( T = 450 K ) , and blue the aging toward equilibrium. In both (a) and (b), the slope of the dashed line is not precisely the number γ of Eq. (6) because the liquids are not perfectly correlating; the line slope is Δ U Δ W ( Δ U ) 2 (see Paper I, Appendix B), a number that is close to γ whenever the liquid is strongly correlating. (c) Virial and potential energy for the asymmetric dumbbell model as functions of time after the temperature jump of (a); in the lower subfigure, data were averaged over 10 ps . Virial and potential energy clearly correlate closely, both on short and long time scales.

FIG. 6.

Computer simulations of virial and potential energy during the aging of two strongly correlating molecular liquids following temperature down-jumps at constant volume ( N V T simulations). (a) The asymmetric dumbbell model at a density of ρ = 1.109 g ml . The liquid was first equilibrated at T = 300 K . Here, simultaneous values of virial and potential energy are plotted, producing the green ellipse, the elongation of which directly reflects the strong W U correlation. Temperature was then changed to T = 240 K where the red ellipse marks the equilibrium fluctuations. The aging process itself is given by the blue points. These points follow the line defined by the two equilibrium simulations, showing that virial and potential energy correlate also out of equilibrium. (b) Similar temperature down-jump simulation of the LW OTP system (Ref. 30). Again, green marks the high-temperature equilibrium ( T = 600 K ) , red the low-temperature equilibrium ( T = 450 K ) , and blue the aging toward equilibrium. In both (a) and (b), the slope of the dashed line is not precisely the number γ of Eq. (6) because the liquids are not perfectly correlating; the line slope is Δ U Δ W ( Δ U ) 2 (see Paper I, Appendix B), a number that is close to γ whenever the liquid is strongly correlating. (c) Virial and potential energy for the asymmetric dumbbell model as functions of time after the temperature jump of (a); in the lower subfigure, data were averaged over 10 ps . Virial and potential energy clearly correlate closely, both on short and long time scales.

Close modal
FIG. 7.

Virial vs potential energy after a temperature down-jump at constant volume applied to SPC water, which is not strongly correlating (colors as in Fig. 6). (a) SPC water at 1 atm equilibrated at T = 260 K , subsequently subjected to an isochoric temperature down jump to T = 160 K . Clearly, W and U are not strongly correlated during the aging process. (b) Same procedure starting from a 2 GPa state point.

FIG. 7.

Virial vs potential energy after a temperature down-jump at constant volume applied to SPC water, which is not strongly correlating (colors as in Fig. 6). (a) SPC water at 1 atm equilibrated at T = 260 K , subsequently subjected to an isochoric temperature down jump to T = 160 K . Clearly, W and U are not strongly correlated during the aging process. (b) Same procedure starting from a 2 GPa state point.

Close modal

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.

FIG. 8.

Crystallization of the supercooled LW OTP liquid where each molecule consists of three LJ spheres with fixed bond lengths and angles (Ref. 30). (a) Pressure (right) and energy (left) monitored as functions of time during crystallization at constant volume. Both quantities were averaged over 1 ns ; on this time scale the pressure/energy fluctuations directly reflect the virial/potential energy fluctuations (Paper II). The horizontal dashed lines indicate the liquid (upper line) and crystal (lower line), the averages of which were obtained from the simulation by averaging over times 0–2 and 5 10 μ s , respectively. Both liquid and crystal show strong correlations, and the correlations are also present during the crystallization. Inset: crystal structure from the simulation. (b) Radial distribution functions of liquid and crystalline phases. The two spikes present in both phases come from the fixed bond lengths.

FIG. 8.

Crystallization of the supercooled LW OTP liquid where each molecule consists of three LJ spheres with fixed bond lengths and angles (Ref. 30). (a) Pressure (right) and energy (left) monitored as functions of time during crystallization at constant volume. Both quantities were averaged over 1 ns ; on this time scale the pressure/energy fluctuations directly reflect the virial/potential energy fluctuations (Paper II). The horizontal dashed lines indicate the liquid (upper line) and crystal (lower line), the averages of which were obtained from the simulation by averaging over times 0–2 and 5 10 μ s , respectively. Both liquid and crystal show strong correlations, and the correlations are also present during the crystallization. Inset: crystal structure from the simulation. (b) Radial distribution functions of liquid and crystalline phases. The two spikes present in both phases come from the fixed bond lengths.

Close modal

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 dynamics32–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.

FIG. 9.

W U plot for the asymmetric dumbbell model for various states at the same volume. The upper right corner shows data for simultaneous values of virial and potential energy for four equilibrium simulations ( T = 240 350 K ) . When quenching each of these to zero temperature in order to identify the inherent states, the crosses are arrived at. The intermediate points are glasses prepared by different cooling rates: out-of-equilibrium systems generated by cooling in 1 ns from 240 K to the temperature in question. This plot confirms the findings of Figs. 6 and 8 that strong virial/potential energy correlations are not limited to thermal equilibrium situations.

FIG. 9.

W U plot for the asymmetric dumbbell model for various states at the same volume. The upper right corner shows data for simultaneous values of virial and potential energy for four equilibrium simulations ( T = 240 350 K ) . When quenching each of these to zero temperature in order to identify the inherent states, the crosses are arrived at. The intermediate points are glasses prepared by different cooling rates: out-of-equilibrium systems generated by cooling in 1 ns from 240 K to the temperature in question. This plot confirms the findings of Figs. 6 and 8 that strong virial/potential energy correlations are not limited to thermal equilibrium situations.

Close modal

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 as6,37

(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

(38)
(39)

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)

(40)
(41)

For the variances, one has

(42)
(43)

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

(44)

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

(45)

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

(46)

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

(47)

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

(48)

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.

Table II.

Check of relation (47) between R T and R E for the LJ and the Kob–Andersen binary Lennard-Jones (KABLJ) liquids. The units for ρ and T are defined in terms of the length and energy parameters σ and ϵ for the interactions of the large particles. The excess isochoric heat capacity was calculated from the potential energy fluctuations in the N V T ensemble.

System ρ T R T C V ex N 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 ρ T R T C V ex N 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 

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

(49)

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

(50)

where

(51)

and

(52)

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

(53)

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 state38 as

(54)

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.,

(55)

This expression is the slope of the pressure versus energy curve at fixed volume, analogous to the γ of Eq. (6) but for the presence of the kinetic terms (recall that for an IPL liquid γ = n 3 = ( W U ) V ). If α p is the coefficient of thermal expansion, K T is the isothermal bulk modulus, and c V = C V V is the isochoric specific heat per unit volume, Eq. (55) implies via standard thermodynamic identities

(56)

This relation allows γ G to be determined from experimentally accessible quantities; in fact, Eq. (56) can be taken as a thermodynamic definition of γ G .40 

There have been suggestions of how to connect the so-called density scaling exponent24,25—the one controlling the relaxation time via the variable ρ γ T —with the Grüneisen parameter, notably by Roland and co-workers.41–43 In Ref. 41, equality of γ G and γ was argued theoretically by reference to the Avramov model. More recently, Roland and Casalini showed43 that equality is not consistent with experimental results; rather γ G is smaller than γ 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 so-called configurational entropy S c ), by introducing an excess or corrected Grüneisen parameter defined via

(57)

Here, Δ c V is the difference in the isochoric specific heats per unit volume between the liquid and the glass. Because Δ c V is smaller than c V , one has γ G corr > γ 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 Δ C V is constant, Roland and Casalini derive

(58)

Hence, γ G corr is the density-scaling exponent.

We take a different approach to connecting the Grüneisen parameter with the slope γ (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 γ G = V ( p E ) V

(59)

Expressing the temperature derivatives in terms of fluctuations (Appendix B of Paper I) gives

(60)

In the limit of strong correlation, one can replace Δ U Δ W with γ ( Δ U ) 2 . Writing the resulting expression in terms of the excess (configurational) specific heat C V ex gives

(61)

In the harmonic approximation, good for many simple liquids close to their melting point,9  C V ex = ( 3 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

(62)

This ratio is around one-half in the harmonic approximation, otherwise larger.

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 energy-bond 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 time-integrated flow, i.e., q ̇ ( t ) = f ( t ) . The energy-bond 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.

Linear-response theory is characterized by the fluctuation-dissipation (FD) theorem, which in the energy-bond formalism is given as follows. Consider a situation with n energy bonds and external control of the effort variables. If f i ( 0 ) f j ( t ) 0 is the equilibrium flow autocorrelation function, the average flow at time t is given by

(63)

If the arbitrary additive constants of the displacements are chosen such that q i = 0 on average, the time-integrated version of this is

(64)

If the flow variables are externally controlled, the FD theorem is

(65)

In most cases, efforts are invariant under time reversal and flows change sign. The Onsager reciprocity relation is f i ( 0 ) f j ( t ) 0 = f j ( 0 ) f i ( t ) 0 (or e i ( 0 ) e j ( t ) 0 = e j ( 0 ) e i ( t ) 0 , depending on which variables are externally controlled and which are free to fluctuate). From the FD theorem, expressions for the frequency-dependent response functions are easily derived. Consider, for instance, the compliances J i j ( ω ) , defined by J i j ( ω ) = q i ( ω ) e j ( ω ) for a periodic situation with infinitesimal perturbations around equilibrium, e ( t ) = Re [ e ( ω ) exp ( i ω t ) ] , etc. For these quantities, the FD theorem implies

(66)

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 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 ) = δ T ( t ) , q 1 ( t ) = δ 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 ) = δ p ( t ) , q 2 ( t ) = δ 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.

  1. The variables of the two energy bonds are linearly constrained as follows:
    (67)
  2. The system’s relaxing properties, i.e., its noninstantaneous responses, are described50 by a single variable ϵ ( t ) as follow:

    (68)

    In these equations the J ’s are the compliances referring to the short-time, nonrelaxing response (the high-frequency response). Note that J 12 = J 21 by the FD theorem.

  3. The relaxing parts of the correlation functions entering into Eq. (66) are proportional. More precisely, the correlation functions obey
    (69)
    and
    (70)

    Note that by differentiation, Eq. (69) implies for t 0 that f 1 ( 0 ) f 1 ( t ) 0 f 1 ( 0 ) f 2 ( t ) 0 f 2 ( 0 ) f 1 ( t ) 0 f 2 ( 0 ) f 2 ( t ) 0 .

  4. The dynamic Prigogine–Defay ratio50 Λ ( ω ) is unity at all frequencies (where double prime denotes the imaginary part),

(71)

Proof that 1 2 . By elimination of the variable ϵ 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 = ( 1 + c ) a , J 12 = J 21 = 1 b , and J 22 = ( d b + a ) b 2 . Introducing the variables ϵ 1 = q 1 J 11 e 1 J 12 e 2 and ϵ 2 = q 2 J 21 e 1 J 22 e 2 , it follows that a ϵ 1 + b ϵ 2 = 0 . This means that we are in the situation described by Eq. (68) with a common relaxing variable to the two energy bonds, ϵ ( t ) ϵ 1 ( t ) ϵ 2 ( t ) , and symmetric short-time compliances, J 12 = J 21 .

Proof that 2 3 . In terms of functional derivatives with respect to the efforts at an earlier time ( t < t ) since ϵ ( t ) for small variations in the effort variables is linear in these, via the FD theorem time-reversal invariance implies that δ q 1 ( t ) δ e 2 ( t ) = δ q 2 ( t ) δ e 1 ( t ) . Thus, Eq. (68) implies that δ ϵ ( t ) δ e 2 ( t ) δ ϵ ( t ) δ e 1 ( t ) . From this, Eqs. (69) and (70) now follow via the FD theorem and Eq. (68).

Proof that 3 4 . According to the FD theorem, the compliance matrix imaginary parts are given by J i j ( ω ) = ( 1 k B T ) 0 q i ( 0 ) f j ( t ) 0 sin ( ω t ) d t . In conjunction with Eqs. (69) and (70), this implies that Λ ( ω ) = 1 at all frequencies.

Proof that 4 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 )

(72)

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 displacement-field relation of the form Eq. (67). Thus, from the energy-bond formalism point of view, a strongly correlating liquid is analogous to the standard measuring cell used for probing ϵ ( ω ) of dipolar viscous liquids, with the strong virial/potential energy correlations reflecting one of Maxwell’s four equations.

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 hydrogen-bonding 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 so-called 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.

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).

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 { U , V } . First, we need the pressure derivatives at constant temperature of U and V . Taking U first, we have (noting that for simple averages like U , it is not necessary to specify the ensemble because of equivalence of ensembles)

(A1)
(A2)

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 U is calculated as follows. The excess (configurational) partition function Z ( V , T ) is the integral (where β = 1 / k B T )

(A3)

Here, Γ indexes points in configuration space and d Γ = d 3 N r V N . In the following, we use the configuration space identity U V = W V valid for changes that scale the coordinates of the microscopic configurations, corresponding to fixed so-called reduced coordinates (compare Appendix A of Paper IV); constant temperature is implicit, as is the dependence of U on Γ and V ,

(A4)
(A5)
(A6)
(A7)
(A8)
(A9)

Thus we have (adding the subscript N V T to the fluctuation expression since this is ensemble dependent)

(A10)

The pressure dependence of V is given by

(A11)

To keep the notation simple, averaging signs are henceforth omitted from simple averages such as V , W , etc. We can write expressions for the variances of U and V in the N p T ensemble using Eq. (49),

(A12)
(A13)

We need also the covariance

(A14)
(A15)

Now we have all we need to construct the H ex V correlation coefficient in the N p T ensemble where H ex = U + p V . The covariance between H ex and V is

(A16)
(A17)
(A18)

where we have used p V = N k B T + W . The variance of H ex is more tedious,

(A19)
(A20)
(A21)

Again using p V = N k B T + W allows some simplication,

(A22)
(A23)

Now we can form the H ex V correlation coefficient,

(A24)
(A25)
(A26)
(A27)

where a = Δ U Δ W N V T + N ( k B T ) 2 and b 2 = K T V k B T ( Δ U ) 2 N V T .

1.
U. R.
Pedersen
,
N. P.
Bailey
,
T. B.
Schrøder
, and
J. C.
Dyre
,
Phys. Rev. Lett.
100
,
015701
(
2008
).
2.
U. R.
Pedersen
,
T.
Christensen
,
T. B.
Schrøder
, and
J. C.
Dyre
,
Phys. Rev. E
77
,
011201
(
2008
).
3.
N. P.
Bailey
,
T.
Christensen
,
B.
Jakobsen
,
K.
Niss
,
N. B.
Olsen
,
U. R.
Pedersen
,
T. B.
Schrøder
, and
J. C.
Dyre
,
J. Phys.: Condens. Matter
20
,
244113
(
2008
).
4.
N. P.
Bailey
,
U. R.
Pedersen
,
N.
Gnan
,
T. B.
Schrøder
, and
J. C.
Dyre
,
J. Chem. Phys.
129
,
184507
(
2008
) (Paper I).
5.
N. P.
Bailey
,
U. R.
Pedersen
,
N.
Gnan
,
T. B.
Schrøder
, and
J. C.
Dyre
,
J. Chem. Phys.
129
,
184508
(
2008
) (Paper II).
6.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon
,
Oxford
,
1987
).
8.
N.
Gnan
,
T. B.
Schrøder
,
U. R.
Pedersen
,
N. P.
Bailey
, and
J. C.
Dyre
,
J. Chem. Phys.
131
,
234504
(
2009
) (Paper IV).
9.
E. D.
Chisolm
and
D. C.
Wallace
,
J. Phys.: Condens. Matter
13
,
R739
(
2001
).
10.
G.
Oster
,
A.
Perelson
, and
A.
Katchalsky
,
Nature (London)
234
,
393
(
1971
).
11.
O.
Klein
,
Medd. Vetenskapsakad. Nobelinst.
5
,
1
(
1919
).
12.
T. H.
Berlin
and
E. W.
Montroll
,
J. Chem. Phys.
20
,
75
(
1952
).
13.
W. G.
Hoover
,
M.
Ross
,
K. W.
Johnson
,
D.
Henderson
,
J. A.
Barker
, and
B. C.
Brown
,
J. Chem. Phys.
52
,
4931
(
1970
).
14.
W. G.
Hoover
,
S. G.
Gray
, and
K. W.
Johnson
,
J. Chem. Phys.
55
,
1128
(
1971
).
15.
Y.
Hiwatari
,
H.
Matsuda
,
T.
Ogawa
,
N.
Ogita
, and
A.
Ueda
,
Prog. Theor. Phys.
52
,
1105
(
1974
).
16.
D.
Ben-Amotz
and
G. J.
Stell
,
J. Chem. Phys.
119
,
10777
(
2003
).
17.
C.
DeMichele
,
F.
Sciortino
, and
A.
Coniglio
,
J. Phys.: Condens. Matter
16
,
L489
(
2004
).
18.
P. E.
Ramirez-Gonzalez
and
M.
Medina-Noyola
,
J. Phys.: Condens. Matter
21
,
075101
(
2009
).
19.
20.
J. D.
Weeks
and
J. Q.
Broughton
,
J. Chem. Phys.
78
,
4197
(
1983
).
21.
J. P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
, 2nd ed. (
Academic
,
New York
,
1986
).
22.
T. B.
Schrøder
,
U. R.
Pedersen
, and
J. C.
Dyre
, e-print arXiv:0803.2199.
23.
T. B.
Schrøder
,
U. R.
Pedersen
,
N. P.
Bailey
,
S.
Toxvaerd
, and
J. C.
Dyre
,
Phys. Rev. E
80
,
041502
(
2009
).
24.
C.
Alba-Simionesco
,
A.
Cailliaux
,
A.
Alegria
, and
G.
Tarjus
,
Europhys. Lett.
68
,
58
(
2004
).
25.
C. M.
Roland
,
S.
Hensel-Bielowka
,
M.
Paluch
, and
R.
Casalini
,
Rep. Prog. Phys.
68
,
1405
(
2005
).
26.
A.
Grzybowski
,
M.
Paluch
, and
K.
Grzybowska
,
J. Phys. Chem. B
113
,
7419
(
2009
).
27.
V.
Molinero
and
E. B.
Moore
,
J. Phys. Chem. B
113
,
4008
(
2009
).
28.
E. A.
Jagla
,
J. Chem. Phys.
111
,
8980
(
1999
).
29.

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 Lewis-Wahnström OTP model (Ref. 30) with the parameters mp=77.106u, σp=0.4963nm, and ϵp=5.726kJmol. A small sphere (mimicking a methyl group) was taken from UA-OPLS having mm=15.035u, σm=0.3910nm, and ϵm=0.66944kJmol. 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).

30.
L. J.
Lewis
and
G.
Wahnström
,
Phys. Rev. E
50
,
3865
(
1994
).
31.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
32.
M.
Goldstein
,
J. Chem. Phys.
51
,
3728
(
1969
).
33.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. A
28
,
2408
(
1983
).
35.
T. B.
Schrøder
,
S.
Sastry
,
J. C.
Dyre
, and
S. C.
Glotzer
,
J. Chem. Phys.
112
,
9834
(
2000
).
36.

In particular, averages of quantities which are the sums of single-particle functions (Ref. 6).

37.
J. L.
Lebowitz
,
J. K.
Perkus
, and
L.
Verlet
,
Phys. Rev.
153
,
250
(
1967
).
38.
M.
Born
and
K.
Huang
,
Dynamical Theory of Crystal Lattices
(
Oxford University Press
,
Oxford U.K.
,
1954
).
39.
N. W.
Ashcroft
and
N. D.
Mermin
,
Solid State Physics
(
Holt, Rinehart and Wiston
,
New York
,
1976
).
40.
D. C.
Wallace
,
Thermodynamics of Crystals
(
Dover
,
New York
,
1972
).
41.
R.
Casalini
,
U.
Mohanty
, and
C. M.
Roland
,
J. Chem. Phys.
125
,
014505
(
2006
).
42.
C. M.
Roland
,
J. L.
Feldman
, and
R.
Casalini
,
J. Non-Cryst. Solids
352
,
4895
(
2006
).
43.
C. M.
Roland
and
R.
Casalini
,
J. Phys.: Condens. Matter
19
,
205118
(
2007
).
44.
H.
Paynter
,
Analysis and Design of Engineering Systems
(
MIT
,
Cambridge, MA
,
1961
).
45.
G. F.
Oster
,
A. S.
Perelson
, and
A.
Katchalsky
,
Q. Rev. Biophys.
6
,
1
(
1973
).
46.
P. V.
Christiansen
,
Dynamik og Diagrammer
(
1978
), IMFUFA Text No. 8, Roskilde.
47.
P. V.
Christiansen
,
Semiotik og Systemegenskaber
(
1979
), IMFUFA Text No. 22, Roskilde.
48.
D. C.
Mikulecky
,
Applications of Network Thermodynamics to Problems in Biomedical Engineering
(
New York University
,
New York
,
1993
).
49.
D. C.
Karnopp
,
D. L.
Margolis
, and
R. C.
Rosenberg
,
System Dynamics: Modeling and Simulation of Mechatronic Systems
(
Wiley
,
New York
,
2006
).
50.
N. L.
Ellegaard
,
T.
Christensen
,
P. V.
Christiansen
,
N. B.
Olsen
,
U. R.
Pedersen
,
T. B.
Schrøder
, and
J. C.
Dyre
,
J. Chem. Phys.
126
,
074502
(
2007
).
51.
T.
Christensen
and
J. C.
Dyre
,
Phys. Rev. E
78
,
021501
(
2008
).
52.
H. J. C.
Berendsen
,
D.
van der Spoel
, and
R.
van Drunen
,
Comput. Phys. Commun.
91
,
43
(
1995
).
53.
E.
Lindahl
,
B.
Hess
, and
D.
van der Spoel
,
J. Mol. Model.
7
,
306
(
2001
).
55.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
56.
B.
Hess
,
H.
Bekker
,
H. J. C.
Berendsen
, and
J. G. E. M.
Fraaije
,
J. Comput. Chem.
18
,
1463
(
1997
).
57.
D.
Coslovich
and
C. M.
Roland
,
J. Chem. Phys.
130
,
014508
(
2009
).