Beyond power-law density scaling: Theory, Simulation, and Experiment

Supercooled liquids are characterized by relaxation times that increase dramatically by cooling or compression. Many liquids have been shown to obey power-law density scaling, according to which the relaxation time is a function of density to some power over temperature. We show that power-law density scaling breaks down for larger density variations than usually studied. This is demonstrated by simulations of the Kob-Andersen binary Lennard-Jones mixture and two molecular models, as well as by experimental results for two van der Waals liquids. A more general form of density scaling is derived, which is consistent with results for all the systems studied. An analytical expression for the scaling function for liquids of particles interacting via generalized Lennard-Jones potentials is derived and shown to agree very well with simulations. This effectively reduces the problem of understanding the viscous slowing down from being a quest for a function of two variables to a search for a single-variable function.

The relaxation time of a supercooled liquid increases markedly upon cooling, in some cases by a factor of ten or more when temperature decreases by just one percent. [1][2][3][4][5][6][7][8][9][10][11] This phenomenon lies behind glass formation, which takes place when a liquid upon cooling is no longer able to equilibrate on laboratory time scales due to its extremely long relaxation time. It has long been known that increasing pressure at constant temperature increases the relaxation time in much the same way as does cooling at ambient pressure. Only during the last decade, however, have large amounts of data become available on how the relaxation time varies with temperature and density. Following pioneering works by Tölle, [12] it was demonstrated by Dreyfus et al., [13] Alba-Simionesco et al., [14], as well as Casalini and Roland, [15] that for many liquids and polymers the relaxation time is a function of a single variable. Roland et al. [16] reviewed the field, and demonstrated that for a large number of molecular liquids and polymers the relaxation time is a function of ρ γ /T , where γ is an empirical material-dependent parameter. For recent works on this "power-law density scaling" or "thermodynamic scaling" see, e.g., Refs. 17-20. The consensus is now that van der Waals liquids and most polymers conform to the scaling, whereas hydrogen-bonding liquids disobey it.
A standard model in simulation studies of viscous liquids is the Kob-Andersen binary Lennard-Jones (KABLJ) mixture, [21] which can be cooled to a highly viscous state and only crystallizes for extraordinarily long simulations. [22] The system consists of 80% large Lennard-Jones (LJ) particles (A) interacting strongly with 20% smaller LJ particles (B). [23] The KABLJ mixture was shown by Coslovich and Roland [24] to obey power-law density scaling to a good approximation with γ = 5.1 for the density range ρ ≡ N/V = 1.15 to ρ = 1.35, whereas Pedersen et al. [25] used the slightly different value γ = 5.16 to scale the density range 1.1 to 1.4. Figure 1 demonstrates, however, that power-law density scaling breaks down when considering a larger density Breakdown of power-law density scaling for the reduced structural relaxation timeτα in the KABLJ mixture,τα ≡ ρ 1/3 (kBT /m) 1/2 τα, where τα is the time at which the self-intermediate scattering function (Fs(q, t), q = 7.25(ρ/1.2) 1/3 ) for the A particles has decayed to 1/e. Molecular dynamics (MD) simulations in the NVT ensemble (N = 1000) were performed using RUMD, a MD package optimized for state of the art GPU computing. [26]τα is plotted for three isochores as a function of the density-scaling variable ρ γ /T , where γ is an empirical scaling parameter. Left panel: γ = 4.90 collapses the data for the two lowest densities. Right panel: γ = 4.45 collapses the data for the two highest densities. It is not possible to find a single exponent that collapses all the data; even though the two exponents differ by only 10%, power-law scaling with a single exponent clearly fails. range. Relaxation time data for the isochores ρ = 1.2 and ρ = 1.6 collapse very well using γ = 4.90, whereas the isochores ρ = 1.6 and ρ = 2.0 collapse using γ = 4.45; in both cases the third isochore deviates significantly. In the following we show that power-law density scaling is an approximation to a more general form of scaling, which is derived from the theory of isomorphs. [27,28] We further show that given that the two lowest densities of Fig. 1 obey power-law density scaling with γ = 4.90, the isomorph theory predicts that the two highest densities scales with γ = 4.45, as seen in Fig. 1.
What causes power-law density scaling and its breakdown for large density variations? A justification of density scaling may be given by reference to inverse powerlaw (IPL) potentials (∝ r −n ) where r is the distance between particles. For such unrealistic, purely repulsive systems density scaling is rigorously obeyed with γ = n/3 [29]. Assuming that power-law density scaling reflects some sort of underlying effective power-law potential, the scaling exponent γ can be found from the NVT equilibrium fluctuations in the potential energy U and the virial W = pV − N kT (IPL potentials have W = (n/3)U ) as follows: This was confirmed for the KABLJ mixture by Refs. 24 and 25. Ref. 25 further supported this "hidden scale invariance" explanation by demonstrating that for the investigated density range the dynamics and structure of the KABLJ mixture is accurately reproduced by an IPL mixture with exponent chosen in accordance with Eq. (1). The hidden scale invariance is not just a feature of the KABLJ mixture but of "strongly correlating liquids" in general. [30][31][32][33] These are defined by having strong correlations between equilibrium NVT fluctuations of the potential energy and the virial (correlation coefficients larger than 0.9). Also molecular models can be strongly correlating; examples include the Lewis-Wahnstrom model of ortho-terphenyl and an asymmetric dumbell model. Both models are strongly correlating and obey power-law density scaling with exponents consistent with Eq. (1) for density increases of 8% and 16% respectively. [34] Very recently Gundermann et al. [20] investigated the van der Waals liquid tetramethyltetraphenyl-trisiloxane, and gave the first experimental confirmation of the relation between the power-law scaling exponent and Eq. (1).
For any potential that is not an IPL potential the exponent γ as calculated from Eq. (1) depends on the state point. Power-law density scaling corresponds to disregarding this state-point dependence. It is thus not surprising that power-law density scaling breaks down when large density changes are involved, but interestingly a more fundamental and robust form of scaling can be derived as we now proceed to show.
Strongly correlating liquids have "isomorphs" in their phase diagram, which are curves along which structure and dynamics in reduced units are invariant to a good approximation. [27,28] The invariance of structure implies invariance of the configurational ("excess") entropy, S ex , thus the isomorphs are nothing but the configurational adiabats. Ref. 27 discussed in detail the consequences of a liquid having isomorphs and also showed that a liquid is strongly correlating if and only if it has isomorphs to a good approximation. The precise statistical-mechanical definition of an isomorph [27] is that this is an equivalence class of state points, where two state points are termed equivalent ("isomorphic") if all pairs of physically relevant micro-configurations of the two state points, which trivially scale into one another, have proportional configurational Boltzmann's factors. From this single assumption several predictions can be derived, including isomorph invariance of structure and dynamics in reduced units and that jumps between isomorphic state points take the system instantaneously to equilibrium. [27] Letting R denote a micro-configuration (all particle coordinates) of a reference state point (ρ * , T * ), the condition for state point (ρ, T ) to be isomorphic with (ρ * , T * ), i.e., the proportionality of pairs of Boltzmann's factors, can, by taking the logarithm and rearranging, be expressed as: where K is a constant that only depends on the two state points. Equation (2) is the basis of the so called direct isomorph check : [27] a) draw micro-configurations R from a simulation at (ρ * , T * ), b) evaluate the potential energies of these configurations scaled to density ρ, and plot them in a scatter-plot against the potential energies at ρ * . If a state point (ρ, T ) exists that to a good approximation is isomorphic with (ρ * , T * ), this scatter plot will be close to a straight line and T can be found from the slope. In the following we consider systems where the interaction potential between particles i and j is given by a sum of inverse power laws: This includes the standard 12-6 LJ potential, but also, e.g., potentials with more than two power-law terms. We note that some systems interacting via Eq. (3) will not be strongly correlating and thus not have good isomorphs. In the following properties are derived for those systems that do have good isomorphs. The total potential energy of a given microconfiguration R at density ρ * is a sum over contributions from the power-law terms, U = n U n . When scaling R to the density ρ, keeping the structure invariant in reduced units, each power-law term is scaled bỹ ρ n 3 = (ρ/ρ * ) n/3 , and the potential energy at the new density U ′ = U ρ −1/3 R is: [28] Thus the linear regression slope of the U ′ , U -scatter plot in the direct isomorph check is given by (where all aver-ages refer to the reference state point (ρ * , T * )): Using Einstein's fluctuation formula for the excess isochoric heat capacity and the corresponding formula for the "partial" heat capacities (which can be negative), we get an expression for the new temperature T relative to the reference temperature T * (compare Eq. (2)): Since C ex v = n C ex v,n the number of parameters in the scaling function g(ρ) is one less than the number of terms in the potential (Eq. (3)). In particular, for the standard 12-6 LJ potential, g(ρ) contains just a single parameter: Using that U 12 = W/2 − U for 12-6 LJ systems, [28] g(ρ) can conveniently be expressed in terms of γ * defined as Eq. (1) evaluated at the reference density ρ * : g(ρ) provides a new and convenient method for numerically tracing out an isomorph -previously this could only be done by changing density by a small amount, e.g., 1%, and then adjusting temperature to keep the excess entropy constant, using that γ (Eq. (1)) also can be expressed as: [27] It is a prediction of the isomorph theory that γ depends only on density. [27,35] This means that the same differential equation, Eq. (10), determines the temperature on all isomorphs, implying that g(ρ) is the same for all isomorphs -what changes between different isomorphs is T * . Thus g(ρ)/T is an isomorph invariant (compare Eq. (7)), which can be used as a scaling variable for the reduced relaxation timeτ that is also an isomorph invariant [27] τ = f (g(ρ)/T ). This form of scaling was first proposed by Alba-Simionesco et al. [14] Here a theoretical derivation has been provided, as well as an explicit expression for g(ρ) for systems interacting via generalized LJ potentials (Eq. (3)). We note further that since the solid-liquid coexistence lines of strongly correlating liquids are predicted to be isomorphs [27,28], this immediately explains the recent observation of Khrapak and coworkers that the solid-liquid coexistence line of an LJ liquid is given by C 1 ρ 4 − C 2 ρ 2 /T = Const. [36]  From the scalingτ = f (g(ρ)/T ) follows that pairs of isochores obey power-law density scaling, since a powerlaw Aρ γ that is equal to g(ρ) at the two densities involved always exists. This is indeed what is seen in Fig. 1. Choosing ρ * = 1.6 as reference density, the scaling in the left panel of Fig. 1 corresponds to g(1.2/1.6) = (1.2/1.6) 4.90 which via Eq. (9) leads to γ * = 4.59. Using this value we find g(2.0/1.6) = 2.70 = (2.0/1.6) 4.45 . Thus from one scaling exponent in Fig. 1, the other is uniquely predicted. Moreover, the value γ * = 4.59 is consistent with what is found by evaluating Eq. (1) at the reference isochore ρ * = 1.6 in the temperature range T = 1.7 to 5, which leads to values of γ * decreasing from 4.6 to 4.5. In the following figures reporting results for the KABLJ mixture, we use γ * = 4.59 as estimated from the left panel of Fig. 1, i.e., no further fitting or adjustment of parameters was applied.
As mentioned, the scaling function g(ρ) was derived assuming that good isomorphs exist. In Fig. 2 we test  1)) to the prediction of the isomorph theory, γ = d ln g/d ln ρ (black curve). Isomorphs denoted by a reduced relaxation time are each generated from a single reference point using Eqs. (7) and (8), with '±' quantifying the resulting variation of the reduced relaxation time on the isomorph. Isochores are plotted with the same symbols as in the main figure. this for the KABLJ mixture using the most sensitive isomorph invariant -the dynamics of viscous states. State points with the same g(ρ)/T , predicted to be on the same isomorph, are seen to have very similar dynamics even though density varies from 1.2 to 2.0. Figure 3 tests the proposed scaling for the KABLJ mixture using the data of Fig. 1. Clearly the new form of scaling works well. Combining Eq. (10) with the definition of g(ρ) (Eq. (7)) shows that γ is the logarithmic derivative of g(ρ), γ = d ln g/d ln ρ. The inset of Fig. 3 demonstrates that this prediction agrees well with simulations, even when going to large densities where the purely repulsive r −12 limit is approached.
We now turn briefly to model molecular systems. In this case it is still a prediction of the isomorph theory that an expression of the form g(ρ)/T is the right scaling variable [35], but we do not have an explicit expression for g(ρ). Figure 4 demonstrates how power-law density scaling breaks down for the Lewis-Wahnstrom model of ortho-terphenyl and an asymmetric dumbell model when considering larger density changes than previously studied. [34] Like in Fig. 1, power-law scaling works when considering pairs of isochores, consistent with the right scaling variable being of the form g(ρ)/T . The insets of Fig. 4 test the isomorph prediction that γ to a good approximation is a function of density only, the assumption used to derive the new scaling. The prediction agrees well with simulations: γ is found to be much more dependent on density than on temperature. For more results on isomorphs in these model molecular liquids see Ref. 37.
Finally, we present experimental results for the two van and (∂ ln γ/∂ ln ρ)T plotted against ln ρ (stars). Consistent with the isomorph theory, γ is found to be much more dependent on density than on temperature.
der Waals liquids dibutyl phthalate (DBP) and decahydroisoquinoline (DHIQ), using larger density increases than usually studied in scaling experiments (20% and 18% respectively). Using earlier reported dielectric data for DBP [38] and DHIQ [39] and values of the Tait equation of state parameters for DBP [40] and DHIQ [41], we find that the isochronal dependences log 10 T vs log 10 ρ determined at given structural relaxation times in reduced units, where v m and m are molecular volume and mass, are nonlinear ( Fig. 5(b) and (e)). This implies breakdown of power-law density scaling. The isochrones can be superposed after vertical shifting, however (Fig. 5(c) and (f)), which implies that a scaling variable exists of the form g(ρ)/T . The isochrones can be described by a phenomenological form of the scal- ing function log 10 g(ρ) = A 1 log 10 ρ + A 2 log 2 10 ρ, chosen simply to take into account the curvature. In the case of DBP, our results are consistent with those reported by Niss et al. [43] We conclude that for DBP and DHIQ power-law density scaling breaks down at large density variations in the way predicted by the isomorph theory. Interestingly, the density dependence of γ is stronger than for the model systems, and in the opposite direction; for the experimental systems γ increases with density.
What are the perspectives of our findings? Based on the theory of isomorphs in dense liquids we have now a form of density scaling that is more fundamental and more robust than power-law density scaling, and which is consistent with both simulations and experiments. This "isomorph scaling" -originally proposed by Alba-Simionesco et al. [14] -is predicted to apply for all strongly correlating liquids, i.e., van der Waals and metallic liquids, but not, e.g., for hydrogen-bonding liquids. Our results should not be used to abandon powerlaw density scaling -it is a useful approximation to isomorph scaling when the scaling function g(ρ) is unkown and only moderate density changes are considered. Under these conditions the isomorph theory predicts that power-law density scaling works with an exponent determined by Eq. (1). Isomorph scaling provides a deeper understanding of why -and when -power-law density scaling works.
Isomorph scaling has important consequenses regarding the most fundamental open question in the field of viscous liquids and the glass transition: what controls the dramatic viscous slowing down? In general this question has to be considered as a search for a physically justified function of two variables, temperature and density (or temperature and pressure). Our results imply that this problem is now effectively reduced to a search for a function of a single variable, at least for the class of strongly correlating liquids. This is particularly striking for LJ type systems like the KABLJ mixture, where we have a prediction for the scaling function, that agrees very well with simulations for much larger density variations than usually considered. The fact that the LJ scaling function contains just a single parameter -i.e., no more parameters than power-law density scaling -confirms that isomorph scaling is more fundamental and not merely a higher-order approximation compared to power-law den-sity scaling. Isomorph scaling must be taken into account in theories of the viscous slowing down: since the relaxation time in reduced units obeys isomorph scaling, any quantity proposed to control the relaxation time must also obey isomorph scaling.