Density scaling as a property of strongly correlating viscous liquids

We address a recent conjecture according to which the relaxation time $\tau$ of a viscous liquid obeys density scaling ($\tau=F(\rho^\gamma/T)$ where $\rho$ is density) if the liquid is ``strongly correlating,'' i.e., has almost 100% correlation between equilibrium virial and potential-energy fluctuations [Pedersen {\it et al.}, PRL {\bf 100}, 011201 (2008)]. Computer simulations of two model liquids - an asymmetric dumbbell model and the Lewis-Wahnstr\"om OTP model - confirm the conjecture and demonstrate that the scaling exponent $\gamma$ can be accurately predicted from equilibrium fluctuations.

Understanding the glass transition depends on understanding the preceding highly viscous liquid phase, where the relaxation time upon cooling approaches and exceeds seconds [1]. Increasing the pressure also leads to much slower relaxations. The study of glass-forming liquids under high pressure has recently become popular, and many data are now available on the properties of the alpha and beta processes, etc., under pressure. If the density is ρ and T is temperature, the last few years have shown that for many highly viscous liquids the (alpha) relaxation time follows the scaling expression The state of this rapidly developing field as of 2005 was summarized in the review Ref. [2] by Roland et al. that presented data for more than 50 liquids and polymers. Equation (1) defines what is referred to as thermodynamic or density scaling. Density scaling is a recent discovery that, following pioneering works by Tölle and Dreyfus et al. [3,4], was proposed as a general principle in 2004 in papers by Alba-Simionesco and co-workers and by Casalini and Roland [5,6]. The former authors demonstrated data collapse at varying temperature and density following a more general expression than Eq. (1) used by the latter authors.
Dreyfus and co-workers found γ = 4 for orthoterphenyl (OTP) and argued that this could be linked to the r −12 repulsive term of the Lennard-Jones potential. It turned out, however, that γ = 4 is not a special exponent [2], leaving the question of the microscopic interpretation of γ open. Coslovich and Roland very recently addressed this by computer simulations of binary Lennard-Jones like liquids where the exponent of the repulsive part of the potential took the values 8, 12, 24, 36 [7]. These model systems obey density scaling and the exponent γ is to a good approximation 1/3 of the exponent characterizing the effective power law of the repulsive core of the potential, as expected for an exact inverse power-law potential [7] (see also Ref. [8]).
Recently, simulations of the thermal equilibrium fluctuations of pressure, energy, and volume in different en-sembles revealed that these quantities correlate strongly for a number of model systems [9,10]. For instance, in the NVT ensemble (i.e., constant volume and temperature) the following systems are "strongly correlating" in the sense that they show more than 90% correlation between virial (the non-kinetic part of the pressure) and potential energy: The standard Lennard-Jones liquid, a liquid with exponential short-range repulsion, the Kob-Andersen binary Lennard-Jones liquid, a sevensite united-atom model of toluene, and the model system studied below consisting of asymmetric "dumb-bell" type molecules. The correlations derive from the repulsive core of the intermolecular potential that, interestingly, dominate fluctuations even at zero and slightly negative pressure [10,11]. For the standard Lennard-Jones liquid the repulsive core is approximately described by a repulsive r −18 term [10,12]. The exponent of the approximate power law depends weakly on state point.
In view of the fact that both density scaling and the strong correlations reflects an effective inverse power-law of the repulsive core of the potential, the following conjecture was proposed in Ref. [10]: A viscous liquid is strongly correlating if and only if it obeys density scaling to a good approximation. There is evidence indirectly supporting this: Highly viscous liquids with strong hydrogen bonds are not strongly correlating because of the "competing interactions" [10], and these liquids also do not obey density scaling very well [2,13,14]. In the present publication simulations are presented that supports the conjecture and strengthens it by adding: For strongly correlating viscous liquids density scaling is obeyed with a scaling exponent γ that can be determined from thermal equilibrium virial and potential energy fluctuations. This amounts to the simplest possible assumption; that it is the same part of the potential that dominates equilibrium fluctuations and flow dynamics.
We performed NVT molecular dynamics simulations [15] of 512 asymmetric dumbbell molecules consisting of pairs of Lennard-Jones (LJ) spheres connected by rigid bonds. The dumbbells were parameterized to mimic toluene [19]. Charges of ±q (specified below) were ap- Slope q=0; ρ=1.006 g/cm^3 q=0; ρ=1.109 g/cm^3 q=0; ρ=1.166 g/cm^3 q=0; P = 2.0GPa q=0.5e; ρ=1.006 g/cm^3 q=0.5e; ρ=1.109 g/cm^3 q=0.5e; ρ=1.166 g/cm^3 6.1 4.9 5.9 FIG. 1: Results from equilibrium molecular dynamics simulations of 512 asymmetric dumbbell molecules with, respectively, a strong dipole moment (q = 0.5e, open symbols, three isochores), and zero dipole moment (q = 0e, filled symbols, three isochores and an isobar). (a) Correlation coeffi- plied to the LJ spheres. The model was simulated with two charges: i) q = 0 corresponding to the simulations done in Refs. [9] and [10]. This version of the model is a "strongly correlating viscous liquid". ii) q = 0.5e (e being the elementary charge) resulting in a dipole moment of 7.0D, i.e., almost 20 times stronger than in toluene, and almost 4 times stronger than water. The purpose of using such a large value of q was to break the correlations and thus to have a version of the model that is not strongly correlating. If the virial is denoted by W , and U is the potential energy, ∆W (t) ≡ W (t) − W and ∆U (t) ≡ U (t) − U , where ... indicates thermal average. The correlation coefficient is defined by R ≡ ∆W ∆U / ∆W 2 ∆U 2 .
R is plotted for a several state points in Fig. 1(a). For q = 0 (filled symbols) all investigated state points have R > 0.95. For q = 0.5e (open symbols) the correlation coefficient is significantly smaller; the Coulomb interactions do indeed break the correlations as expected [10].
We define γ ≡ (∆W ) 2 / (∆U ) 2 . If R ≈ 1 it follows that ∆W (t) ≈ γ∆U (t) in their instantaneous fluctuations [10], and consequently we refer to γ as the 'slope'. According to our conjecture, the slope γ is also the scaling exponent in Eq. (1). In Fig. 1(b) we show the slopes for all investigated state points. For q = 0 (filled symbols) there is a small, but significant dependence on density and temperature. Thus if the conjecture is correct, density scaling can not be exact: To apply density scaling a single value of γ is needed, but the γ we get from the fluc-  tuations depends slightly on the state point [10]. In the following we consider for q = 0 two values of the scaling exponent: the slope averaged over all state points with q = 0; γ = 6.1, and the slope averaged over the three data sets with the smallest slopes for q = 0 (ρ = 1.109 g/cm 3 , ρ = 1.166 g/cm 3 and P=2.0GPa); γ = 5.9, i.e., the 'best' compromise if we chose to ignore the ρ = 1.006 g/cm 3 isochore. For q = 0.5e the slopes are less density dependent with a mean value γ = 4.9.
In the following we apply density scaling to the diffusion coefficient estimated from the long-time behavior of the mean-square displacement, ∆r 2 (t) , of the large spheres (the "phenyl group") [22]. The diffusion coefficients for all state points studied are given in Fig. 2.
Following Coslovich and Roland [7] we apply density scaling to the reduced diffusion coefficient, D * ≡ (N/V ) 1/3 (k B T /m) −1/2 D where m is the mass of the molecules. In Fig. 3(a) D * is plotted for q = 0 as a function of 1000ρ γ /T with four different values of γ. Clearly, density scaling works neither with γ = 7.0 or γ = 5.0. Comparing the scaling with γ = 6.1 to the data without scaling (filled symbols in Fig. 2), we find here good data collapse; by far most of the density dependence is captured by the density scaling with γ = 6.1. With γ = 5.9 the data collapse is even better for three of the data sets, whereas one data set deviates slightly from the master curve comprised of these three sets. This is the isochore ρ = 1.006 g/cm 3 , i.e., the one that was ignored when choosing γ = 5.9 ( Fig. 1(b)).
The conclusion drawn from Fig. 3 is two-fold: i) Density scaling is approximate (as discussed above) -for a larger region of state points scaling will be less perfect. ii) For a given range of state points, the scaling exponent can be found by studying equilibrium fluctuations.
In Fig. 4 the reduced diffusion coefficients D * for q = 0.5e are plotted as a function of 1000ρ γ /T with three different values of γ. The value γ = 4.9 chosen from the equilibrium fluctuations ( Fig. 1(b)), is found to be a reasonable scaling exponent. However, as conjectured, the data collapse achieved is inferior to that for the strongly correlating version of the model (q = 0).
To test the generality of our findings, we repeated the analysis for the Lewis-Wahnström OTP model (LW-OTP) [15,23] (see also Refs. [24] and [25]). The results achieved for this model (Fig. 5) are qualitatively similar to the results for the asymmetric dumbbell model with q = 0; (i) LW-OTP is strongly correlating (0.91 < R < 0.92) [26], (ii) the slope is slightly state point dependent, (iii) choosing the average slope as scaling exponent γ gives good data collapse.
In summary, we have presented numerical evidence for the conjecture that density scaling is a property of strongly correlating viscous liquids. For the two strongly correlating models investigated density scaling applies with a scaling exponent that can be accurately predicted from the equilibrium fluctuations. This represents a step forward in the theoretical understanding of density scaling. In particular, the scaling exponent γ should no longer be regarded as an empirical fitting parameter. In computer simulations γ can be estimated directly from the equilibrium fluctuations. Via the fluctuation dissipation theorem the equilibrium fluctuations determine the frequency-dependent linear thermoviscoelastic response functions -this fact provides a possible future route for independent experimental estimation of the scaling exponent γ [9,27].