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 and pressure have contributions both from particle momenta and positions,
For a liquid with pair interactions, if is the pair potential and is the distance between particles and , we have
Strong , 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 , defined1,4 by
Here and henceforth, unless otherwise specified, angular brackets ⟨ ⟩ denote thermal ensemble averages and denotes deviation from the average value of the quantity in question. We call liquids with strongly correlating. Another characteristic quantity is the “slope” , which we here define1,4 as the ratio of standard deviations,
In the limit of perfect correlation , becomes equal to the slope of average as a function of average at fixed volume, .
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 are perfectly correlating, the results of this analysis can be summarized as follows: (1) almost all of the fluctuations in and come from interparticle separations in the region of the first peak of the radial distribution function ; (2) in this region, the LJ potential is approximated very well by the sum of an IPL with exponent and a linear term ; and (3) when volume is fixed, the parts of and that come from the linear term are almost constant. Our initial and simpler explanation of strong 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 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 ensemble than in the 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 , showing that the slope [Eq. (6)] is larger than 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
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 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 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 identical particles in volume interacting by a pair potential of the form ; 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 is conveniently written6,21 as the ideal gas term plus the nontrivial “excess” free energy, . The first term is the free energy of an ideal gas at the same volume and temperature, , where is the particle number density and is the thermal de Broglie wavelength. The excess free energy is given6,21 by
Whenever , this expression leads to a free energy with a well-defined extensive thermodynamic limit .11,12
Equation (8) implies that a number of derived quantities are also functions of . As important examples, recall the following standard identities: the excess entropy: , the potential energy: , the virial , the excess isothermal bulk modulus: , the excess isochoric specific heat per unit volume: , the excess pressure coefficient: . These quantities are all functions of the single variable ; more accurately, one has [where , , , and ],
The functions and all depend on , but, for simplicity of notation, we have not indicated this explicitly. Dividing across by the dimensional factors on the right hand side (for example, in the case of potential energy and virial), one arrives at dimensionless forms of the excess entropy, potential energy, etc., which are functions of only.
Turning now to the dynamics, consider the standard molecular dynamics (MD) case where the equations of motion are Newton’s equations. Suppose is a solution to Newton’s equations. Then, it is straightforward to show that is also a solution if . In particular, if refers to equilibrium ( or ) dynamics at a state point with density and temperature , then refers to equilibrium dynamics at density and temperature (temperature scales as the mean-square velocity and velocities get a factor ). Using the above relation between and , this implies that
This means that two states with different densities and temperatures but the same have dynamics that scale into one another by simple scalings of space and time. In particular, if for any quantity one defines the relaxation time via , it follows that any two states with same have the same “reduced” (dimensionless) relaxation time , if this quantity is defined by where the characteristic time is defined by . Similarly, if one defines the reduced diffusion constant by , where , then 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 are identical provided that lengths are scaled by .
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 and , rewriting it (as can always be done) as a sum of an IPL potential plus the difference, denoted as “diff,”
For any microscopic configuration , 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 to a good approximation is only a function of volume: —at least for states that carry Boltzmann weights of any significance. The approximate identity 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 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 . For the pure IPL system , the two state points have the same dynamics and structure as argued in the previous section. Letting 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 (using 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 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 single-pair effects are insufficient to explain the strong 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 -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 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, indicates an optional constant. To start with, how should the exponent and the coefficients and be chosen? An obvious choice, followed in the first part of Paper II, is to require that the two potentials, and , should agree as much as possible around a particular value of , denoted . Given , if we require the functions and their first two derivatives to match at , this determines all three parameters , , and . The exponent is given (Paper II) by
Here, the notation refers to one kind of -dependent effective IPL exponent, based on the ratio of the second and first derivatives. For this simply returns . Otherwise, it gives a local matching of the to . This leaves effectively one parameter to vary, namely , which must be less than the location of the minimum of the LJ potential , where diverges. The parameter may be chosen to optimize the match of the fluctuations in the total energy and virial. For an simulation at and near-zero pressure, the best fit was obtained with (while the exponent implied by the slope [Eq. (6)], , 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 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 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 -dependent effective exponent is the higher order defined by
This also returns for , but since it does not involve the first derivative, it returns even if a linear function of 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 , , and . 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 and . For simplicity, we take the exponent directly from the observed fluctuations: , where is defined in Eq. (6). To fix the coefficient , 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 and the power-law potential energy , similarly the corresponding virials and . The difference quantities and are defined as
A perfect match of the fluctuations would mean that and have zero variance. Therefore, we choose to minimize the sum of the relative diff variances,
For the near-zero pressure state point used in Fig. 1 of Paper I, the exponent determined from is and the optimal value of is . Before examining the difference potential, what do we get if we simulate with the matched IPL potential? Figure 2 shows the radial distributions obtained for the above state point and another with a higher density and temperature, for three potentials: LJ, the repulsive term of the LJ potential, and the optimal IPL potential with . We used the same IPL potential at both state points (i.e., we did not adjust and to match the second state point). The first thing to note is that the potential gives a structure much closer to that of the LJ potential than does the repulsive 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 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 is nearly linear. Thus, a good approximation to the real potential is the eIPL potential of Eq. (33) for less than a cutoff and zero otherwise. Neglecting the small value of , the cutoff is given by . For the fit shown in Fig. 3(b), .
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 , which means that argument that the linear term contributes little to the fluctuations depends on all bond lengths being less than . 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 . 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 as the lengths of bonds fluctuate back and forth across , so the sum of bond lengths, which are less than , 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 does go to zero after the peak—the correlation coefficient becomes very high, over 99.5%, as (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 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 ). The 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 simulation at the same state point, i.e., with pressure chosen as the average pressure of the corresponding 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 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 , while the linear term gets multiplied by 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 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 case, on the other hand, the diff contributions to the fluctuations of are more than double the IPL ones; this is not the case for the diff contributions to , though they are still a larger fraction of the total than in the 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 is smaller than the sum of those of and , in the case, indicates that the latter two are negatively correlated. This is presumably due to bond lengths around which alternately are counted as part of and as part of as they fluctuate. This effect is less noticeable in the case; there we see clearly that fluctuations in account for most of those in .
|Ensemble .||Quantity .||LJ .||IPL .||diff .||lin .||rest .|
|Ensemble .||Quantity .||LJ .||IPL .||diff .||lin .||rest .|
Based on the above, we can now answer the question: Is it possible to predict whether or not a liquid is strongly correlating by inspection of its potential (i.e., without simulating virial and potential energy fluctuations)? For liquids with particles interacting by pair potentials, the answer is in the affirmative: the liquid is strongly correlating if the potential around the first peak of the structure factor (the typical interparticle distance) may be fitted well by an eIPL and if the potential involves only one characteristic length. For more general potentials, the situation is more complex. Thus, it is possible to construct many-body potentials with angular dependencies which scale like IPL pair potentials; these have 100% 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 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).
IV. OUT-OF-EQUILIBRIUM DYNAMICS IN MOLECULAR MODELS
According to the eIPL explanation detailed in Sec. III, strong 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 fluctuations.
A. Temperature down-jump simulations of three molecular model liquids
Figure 6(a) shows the results for a temperature down-jump at constant volume, starting and ending in equilibrium ( 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 . The green ellipse is the scatter plot of simultaneous values of and in equilibrium at . The strong correlation is revealed by the elongation of the ellipse (; ). When the liquid is similarly equilibrated at , the red blob appears. To test for correlation in an out-of-equilibrium situation, we changed temperature abruptly from the equilibrium situation to . The blue points show how virial and potential energy evolve following the temperature down-jump. Clearly, strong 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 (; ). The colors are as in Fig. 6(a): green gives a high-temperature equilibrium state , red gives a low-temperature equilibrium state , and the blue points show the aging toward equilibrium after changing temperature from . 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 and after the temperature jump for the data in Fig. 6(a) for the asymmetric dumbbell liquid. and 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 (, ); in the second case, correlations are somewhat stronger (, ), 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, and 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 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 and . 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 , so strong 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.
C. Glasses and inherent states
The above out-of-equilibrium simulations show that the property of strong 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 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.
V. ENSEMBLE DEPENDENCE OF THE CORRELATION COEFFICIENT
Most of the simulations of Paper I were done in the 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 and ensembles, then the and 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 held fixed and one with its conjugate intensive variable held fixed (defined so their product is dimensionless). The other parameters defining the ensembles are the same. The covariance of observables and in the two ensembles are then related as6,37
To compare the and ensembles, we take and as the energy and the inverse temperature , respectively, keeping the volume fixed. Noting that and , where is the extensive isochoric specific heat, , the covariance of and is given by
Here, we have introduced the excess parts of the isochoric specific heat and pressure coefficient, and , 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 -correlation coefficient (in the following subscripts or indicate the or ensemble, respectively) is given by
We wish to express the right side in terms of the coefficient . The definition of implies that
Inserting this into Eq. (44) and making use of the fluctuation relations and (see, e.g., Appendix B of Paper I) gives
After squaring and cancelling factors of and , we get an expression relating to ,
To get a feel for the relation, divide by . This yields
First one notes that when , the denominator on the right side becomes equal to the numerator and . That is, the property of perfect correlation is independent of (fixed volume) ensemble. When , the denominator becomes greater than the numerator, and so . That is, the correlation coefficient is smaller in the ensemble than in the one. How can we understand this? Consider the set of points sampled by the system during an trajectory; this is an elongated blob in the 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 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 .||.||.||.||.||[Eq. (47)] .||(measured) .|
|System .||.||.||.||.||[Eq. (47)] .||(measured) .|
In the ensemble, where volume is allowed to fluctuate, we must consider different variables. The natural variables to correlate are the excess enthalpy and the volume . We use again Eq. (37), but now take as and as , 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 correlation coefficient is
Notice that is strictly less than unity—even for perfectly correlating liquids (that is, with perfect correlations in the and ensembles). For the LJ simulation of Fig. 1(b), (the correlation coefficient for the same state point is 0.94). Unlike the situation when comparing the and ensembles, there does not seem to be a simple relation between the two correlation coefficients. It seems likely, though, that the correlation in the ensemble is generally smaller than the correlation in the ensemble; thus, the property of strong correlation is less evident in the 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 is the frequency of the normal mode. By assuming that is the same for all modes and denoting the common value by , one can derive the Mie–Grüneisen equation of state38 as
Here, is the pressure, with is the “static” energy of the crystal per atom (the energy of the force-free configuration about which vibrational motion occurs), and is the vibrational energy. In general depends on volume, but this dependence is typically small enough that it can be neglected. From Eq. (54) it follows that, if is the total, thermally averaged internal energy, one has , i.e.,
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 ). If is the coefficient of thermal expansion, is the isothermal bulk modulus, and is the isochoric specific heat per unit volume, Eq. (55) implies via standard thermodynamic identities
There have been suggestions of how to connect the so-called density scaling exponent24,25—the one controlling the relaxation time via the variable —with the Grüneisen parameter, notably by Roland and co-workers.41–43 In Ref. 41, equality of 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 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 ), by introducing an excess or corrected Grüneisen parameter defined via
Here, is the difference in the isochoric specific heats per unit volume between the liquid and the glass. Because is smaller than , one has . 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 is constant, Roland and Casalini derive
Hence, 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
Expressing the temperature derivatives in terms of fluctuations (Appendix B of Paper I) gives
In the limit of strong correlation, one can replace with . Writing the resulting expression in terms of the excess (configurational) specific heat gives
In the harmonic approximation, good for many simple liquids close to their melting point,9 (while it is generally larger in the supercooled state). Thus, the term 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 one-half in the harmonic approximation, otherwise larger.
B. Energy-bond formulation of the strongly correlation property
This section derives a general thermodynamic condition of the property of strong 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 and a “flow” variable , where is the free energy transferred into the system per unit time. The “displacement” is the time-integrated flow, i.e., . 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 energy bonds and external control of the effort variables. If is the equilibrium flow autocorrelation function, the average flow at time is given by
If the arbitrary additive constants of the displacements are chosen such that on average, the time-integrated 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 (or , 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 , defined by for a periodic situation with infinitesimal perturbations around equilibrium, , 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 .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 (, ), the mechanical energy bond has pressure variation as the effort variable and the negative volume variation as the displacement variable (, ). 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)
The system’s relaxing properties, i.e., its noninstantaneous responses, are described50 by a single variable as follow:(68)
In these equations the ’s are the compliances referring to the short-time, nonrelaxing response (the high-frequency response). Note that by the FD theorem.
- 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 that .
The dynamic Prigogine–Defay ratio50 is unity at all frequencies (where double prime denotes the imaginary part),
Proof that . 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 and are dimensionless. Define , , and . Introducing the variables and , it follows that . This means that we are in the situation described by Eq. (68) with a common relaxing variable to the two energy bonds, , and symmetric short-time compliances, .
Proof that . In terms of functional derivatives with respect to the efforts at an earlier time since for small variations in the effort variables is linear in these, via the FD theorem time-reversal invariance implies that . Thus, Eq. (68) implies that . From this, Eqs. (69) and (70) now follow via the FD theorem and Eq. (68).
Proof that . According to the FD theorem, the compliance matrix imaginary parts are given by . In conjunction with Eqs. (69) and (70), this implies that at all frequencies.
Proof that . 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 )
How does this all relate to strong 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.
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 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 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).
APPENDIX: CALCULATING THE CORRELATION COEFFICIENT IN THE ENSEMBLE
Here, we provide the details of the calculation of the correlation coefficient between volume and excess enthalpy in the ensemble. We apply Eq. (49) with . First, we need the pressure derivatives at constant temperature of and . Taking first, we have (noting that for simple averages like , it is not necessary to specify the ensemble because of equivalence of ensembles)
Note that in the derivative is without averaging signs since there it is a parameter of the relevant ensemble . The volume derivative of is calculated as follows. The excess (configurational) partition function is the integral (where )
Here, indexes points in configuration space and . In the following, we use the configuration space identity 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 on and ,
Thus we have (adding the subscript to the fluctuation expression since this is ensemble dependent)
The pressure dependence of is given by
To keep the notation simple, averaging signs are henceforth omitted from simple averages such as , , etc. We can write expressions for the variances of and in the ensemble using Eq. (49),
We need also the covariance
Now we have all we need to construct the correlation coefficient in the ensemble where . The covariance between and is
where we have used . The variance of is more tedious,
Again using allows some simplication,
Now we can form the correlation coefficient,
where and .
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 , , and . A small sphere (mimicking a methyl group) was taken from UA-OPLS having , , and . The bonds were kept rigid with a bond length of . The volume was , giving an average pressure of approximately . The temperature was held constant at using the Nosé–Hoover thermostat. 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).