Hidden scale invariance in molecular van der Waals liquids: A simulation study

Results from molecular dynamics simulations of two viscous molecular model liquids -- the Lewis-Wahnstrom model of ortho-terphenyl and an asymmetric dumbbell model -- are reported. We demonstrate that the liquids have a ``hidden'' approximate scale invariance: Equilibrium potential energy fluctuations are accurately described by inverse power law (IPL) potentials, the radial distribution functions are accurately reproduced by the IPL's, and the radial distribution functions obey the IPL predicted scaling properties to a good approximation. IPL scaling of the dynamics also applies -- with the scaling exponent predicted by the equilibrium fluctuations. In contrast, the equation of state does not obey the IPL scaling. We argue that our results are general for van der Waals liquids, but do not apply, e.g., for hydrogen-bonded liquids.


I. INTRODUCTION
A phenomenon is scale invariant if it has no characteristic length or time. Scale invariance emerged as a paradigm in the early 1970s following the tremendous successes of the theory of critical phenomena. Liquids described by scaleinvariant potentials are a theorist's dream by having a number of simple properties ͓1͔. If potentials vary with distance as r −n , the configurational free energy F as a function of density ϵ N / V and temperature T may be written F = Nk B Tf͑ ␥ / T͒ where ␥ = n / 3. In fact, all properties of inverse power-law ͑IPL͒ liquids-in appropriately defined reduced units-are functions of ␥ / T only. Even dynamical properties are simple, e.g., the relaxation time may be written as with the same exponent ␥ that characterizes the thermodynamics, t 0 being the reduced time unit ͑t 0 ϵ −1/3 ͱ m / k B T͒.
Unfortunately the predicted equation of state does not fit data for real fluids, and IPL liquids do not have low-pressure liquid states. Both problems derive from the absence of molecular attractions; these define the spatial scale of one intermolecular distance at low and moderate pressure. Thus real liquids have apparently little in common with scale-invariant liquids. In this paper we demonstrate by example that there is nevertheless a "hidden" approximate scale invariance in a large class of liquids, the van der Waals liquids. Consequently several properties of IPL liquids apply to van der Waals liquids, leading to a number of experimentally testable predictions.

II. STRONGLY CORRELATING LIQUIDS
It was recently shown ͓2-4͔ that several model liquids exhibit strong correlations between the equilibrium fluctuations of the potential energy U and the virial W ͑defining the contribution to pressure coming from the intermolecular interactions via pV = Nk B T + W͒. "Strongly correlating liquids" include ͓4͔ the standard Lennard-Jones liquid, a liquid with exponential short-range repulsion, the Kob-Andersen binary Lennard-Jones liquid, a seven-site united-atom model of toluene, the three-site Lewis-Wahnström model of orthoterphenyl ͑OTP͒, as well as a model consisting of asymmetric "dumbbell" type molecules.
For inverse power-law ͑IPL͒ potentials ϰr −n there is 100% correlation between virial and potential energy fluctuations as function of time: ⌬W͑t͒ = ␥⌬U͑t͒ where ␥ = n / 3 ͓1͔. By reference to the single-component Lennard-Jones system it was argued in Ref. ͓4͔ that strongly correlating liquids are well approximated by inverse power-law potentials as regards their thermal equilibrium fluctuations. The present paper presents concrete proofs that this is the caseeven for molecular model liquids with van der Waals type interactions ͑i.e., excluding hydrogen, ionic, and covalent bonds͒. The main conclusions below are: ͑i͒ At a given density and temperature the energy surface of a strongly correlating molecular liquid may be replaced to a good approximation by that of inverse power-law potentials.
͑ii͒ Strongly correlating liquids "inherit" a number of scaling properties from inverse power-law potentials, but not all such properties.

III. FLUCTUATIONS AND STRUCTURE
We performed NVT molecular dynamics simulations ͓5-8͔ of two molecular liquids: ͑i͒ 512 asymmetric dumbbell molecules consisting of pairs of different sized Lennard-Jones ͑LJ͒ spheres connected by rigid bonds; the dumbbells were parameterized to mimic toluene ͓9͔; ͑ii͒ the Lewis-Wahnström OTP model ͓10͔ ͑N = 324͒. Both models are strongly correlating for the range of state points investigated here. For the dumbbell model we find for the correlation coefficient 0.95Ͻ R Ͻ 0.97; for the OTP model 0.91Ͻ R Ͻ 0.92. Figure 1͑a͒ illustrates the strong WU correlation for the asymmetric dumbbell model. Since R Ϸ 1 it follows that ⌬W͑t͒Ϸ␥⌬U͑t͒ in their instantaneous fluctuations ͓3͔, consequently ␥ is referred to as the "slope." Figure 1͑b͒ details the origin of the strong correlations, focusing on the strongest interaction in the dumbbell model, the "BB" interaction ͑"B" indicating the large "phenyl" spheres͒. The fluctuations in U BB ͑t͒ are well described by an IPL pair potential IPL,BB ͑r͒ ϰ r −18 . The same conclusion ap-plies for the AA and AB interactions ͑data not shown͒. The reason that Diff,␣␤ ͑r͒ϵ LJ,␣␤ ͑r͒ − IPL,␣␤ ͑r͒ contributes little to the fluctuations at constant volume is that this term is approximately linear in the region of the first peak of the pair correlation function. This implies that, when one pair distance is increased, other pair distance͑s͒ decrease by approximately the same amount, keeping the change in U Diff,␣␤ ͑t͒ ϵ͚ iϾj Diff,␣␤ ͓r ij ͑t͔͒ small ͓4͔. Switching the ensemble from NVT to NpT ͑p = 1.32 GPa͒ reduces the WU correlation from 0.959 to 0.866, illustrating the importance of the constant volume constraint ͑the importance increases with decreasing pressure͒.
The fluctuations tell us that the relevant part of the potential energy surface is well approximated by a corresponding "IPL" system, i.e., the system where the LJ pair potentials are replaced by the IPL potentials discussed in Fig. 1͑b͒. To test how far this correspondence holds, we simulated the IPL system at the same density and temperature. Figure 2͑a͒ compares the pair distribution functions. The agreement is striking. To the best of our knowledge this is the first proof that the structure of a molecular liquid is well reproduced by an IPL liquid. In contrast, the average pressures of the two models are quite different ͓see Fig. 2͑a͔͒, because the IPL system does not include contributions to the pressure from Diff,␣␤ ͑r͒. Consider now changing density and temperature, conserving the parameter ⌫ = ␥ / T. For a pure n =3␥ IPL liquid this means that one is effectively studying the same system, but with scaled length and time units ͑times scaled by t 0 , lengths by l 0 ϵ −1/3 ͒ ͓1͔: two state points with the same ⌫ have 3N-dimensional potential energy surfaces that are identical in reduced units. Switching back to the LJ system amounts to adding the "diff" terms of the potentials, which-to a good approximation-simply shifts the whole energy surface by an amount that only depends on the volume of the state point ͓4͔. This changes neither structure nor dynamics, and consequently these properties are predicted to "inherit" the scaling properties of the IPL potentials. In contrast, the shift of the energy surface does change the absolute values of potential energy and pressure ͑because the shift depends on volume͒, and these properties are not predicted to inherit the IPL scaling.
The IPL scaling of the structure is tested in Fig. 2͑b͒ for the asymmetric dumbbell model, and in Fig. 3 for the LW-OTP model. For both models the IPL scaling works well; there is only a small effect of the fixed bond connecting the LJ spheres. Note that average pressure and potential energy do not follow the IPL scaling.

IV. DYNAMICS
We now turn to the long-time dynamics. To the degree that the IPL approximation holds over a range of state points, the long-time dynamics is a function of ⌫ = ␥ / T where ␥ is one third of the IPL exponent n, i.e., the "density" scaling of FIG. 1. ͑Color online͒ ͑a͒ ⌬W͑t͒ / N vs ⌬U͑t͒ / N for the asymmetric dumbbell model ͓9͔. The strong correlations are quantified by the correlation coefficient, R ϵ͗⌬W⌬U͘ / ͱ ͗͑⌬W͒ 2 ͗͑͘⌬U͒ 2 ͘ = 0.959. The dashed line has the slope ␥ ϵ ͱ ͗͑⌬W͒ 2 ͘ / ͗͑⌬U͒ 2 ͘ = 5.90. ͑b͒ For configurations taken from the simulation reported in Fig. 1͑a͒ we evaluated U IPL,BB ͑t͒ϵ͚ iϾj IPL,BB ͑r ij ͑t͒͒ with IPL,BB ͑r͒ϵC BB ⑀ BB ͑ BB / r͒ 18 . The constant C BB = 1.489 ͑C AA = 1.075, C AB = 1.237͒ was chosen to make the variance of U IPL,BB ͑t͒ equal to the variance of U BB ͑t͒ϵ͚ iϾj LJ,BB ͓r ij ͑t͔͒. U IPL,BB ͑t͒ strongly correlates with U BB ͑t͒ with a correlation coefficient of 0.977. In comparison U Diff,BB ͑t͒ϵU BB ͑t͒ − U IPL,BB ͑t͒ has small variance, 0.045 times the variance of U BB ͑t͒ and is only weakly correlated with the latter. Redoing the analysis with n =3 ϫ 5.90 ͓the IPL exponent suggested by the fluctuations, see Fig.  1͑a͔͒ gives almost identical results.
FIG. 2. ͑Color online͒ ͑a͒ Radial distribution functions for the asymmetric dumbbell model ͓9͔ ͑full lines͒, and the corresponding r −18 IPL model ͑dashed lines, see text͒. ͑b͒ Scaled g BB ͑r͒ for three state points; the two state points with the same ⌫ = 5.9 / T ͑black full line and red long dashed line͒ have almost identical ͑scaled͒ structure. When increasing at constant ⌫ there is a slight shift of the first peak to smaller distances, consistent with the existence of the fixed bond. Similar scaling is found for the AA and AB interactions ͑data not shown͒.
Eq. ͑1͒ is predicted to apply. This type of scaling has in recent years become a well established empirical fact for a large number of viscous liquids ͓11,12͔ ͑usually ignoring the insignificant state-point dependence of t 0 ͒. In particular, density scaling applies for van der Waals liquids, but not for liquids with hydrogen bonds or other directional bonds ͓13-16͔. Other forms of scaling have been suggested, but it has been argued that the scaling reflects an underlying IPL ͓3,12,17-20͔, and thus that Eq. ͑1͒ is the correct scaling.
A potential complication for the IPL explanation of density scaling, is that the IPL exponent suggested by the fluctuations ͑n =3␥͒ varies somewhat with state point ͓4,21͔. This is illustrated for the two models in Fig. 4. Note that, while the suggested exponent is weakly state-point dependent for both models, one finds quite different exponents for the two models. For the OTP model we test density scaling using the slope found by averaging over all state points: ␥ = 7.9. For the dumbbell model we consider two values of the scaling exponent: the slope averaged over all state points, ␥ = 6.1, and the slope averaged over the three data sets with the smallest slopes, ␥ = 5.9, i.e., the "best" compromise if the = 1.006 g / cm 3 isochore is left out of consideration. We ap-    ply 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 and D is the diffusion coefficient estimated from the long-time mean-square displacement. Figure 5 applies density scaling to the Lewis-Wahnström OTP model. The scaling works neither with ␥ = 9.0 nor with ␥ = 7.0, whereas it works well with ␥ = 7.9 which is the exponent we calculated from the fluctuations ͑Fig. 4͒. In Fig. 6 density scaling is applied to the asymmetric dumbbell model. Density scaling works neither with ␥ = 7.0 nor with ␥ = 5.0. Comparing the scaling with ␥ = 6.1 to the data without scaling ͑␥ = 0, open symbols͒, we find quite good data collapse-by far most of the density dependence is captured by using ␥ = 6.1. With ␥ = 5.9, however, the data collapse is even better for three of the data sets, whereas one data set deviates from the master curve comprised of these three sets. The latter is the isochore = 1.006 g / cm 3 , i.e., precisely the one that was ignored when choosing ␥ = 5.9 ͑Fig. 4͒.
The above results demonstrate that density scaling is a consequence of an underlying IPL scale invariance which is revealed by studying thermal equilibrium fluctuations at a single state point. Obviously the scaling with a single value of ␥ can only be perfect for true IPL liquids. For "real" liquids a larger region of state points makes the scaling less perfect for two reasons ͑compare Fig. 6, ␥ = 6.1 and 5.9͒: the apparent IPL exponent ͑n =3␥͒ changes more, and the relative effect of intramolecular interactions that do not obey the IPL scaling ͑here the fixed bonds͒ is larger. In this connection it is important to note, however, that even small density changes are experimentally relevant. Thus changing the density of, e.g., phenylphthalein-dimethylether by just 2% can change the relaxation time by more than two orders of magnitude ͓22͔.
As a consequence of the underlying IPL scale invariance, state points with same ␥ / T are expected to have not only the same scaled structure and reduced diffusion coefficient as demonstrated above, every aspect of the dynamics should be the same in reduced units. Figure 7 shows the reduced meansquare displacement and rotational correlation coefficient of the LW-OTP model, at the three state points of Fig. 3. The agreement between the two state points with same ␥ / T ͑full black line and red crosses͒ is striking. In contrast, the third state point ͑blue dashed line͒ has the same temperature as the first state point ͑full black line͒ and the same pressure as the second ͑red crosses͒, but very different dynamics. Clearly, neither pressure nor temperature controls the dynamics-to a very good approximation the controlling parameter is ⌫ = ␥ / T.
We argued above that the IPL scaling properties demonstrated are a feature only of strongly correlating liquids. To illustrate this we modified the asymmetric dumbbell model to make it less correlating. This was done by adding charges of q = Ϯ 0.5e ͑e being the elementary charge͒ to the LJ spheres, resulting in a dipole moment of 7.0 D, i.e., almost four times stronger than water. The correlation coefficients and slopes are shown Fig. 8 ͑open symbols͒, compared to the  Fig. 3. ͑a͒ Meansquare displacement for the noncentral LJ spheres. ͑b͒ Rotational correlation function for the vector perpendicular to the plane of the three LJ spheres, using the first Legendre polynomial. q = 0 version of the model ͑filled symbols͒. The addition of the large dipole moments decreases the correlation, as expected. As before we use the averaged slope as the scaling exponent suggested by the fluctuations, giving ␥ = 4.9 for q = Ϯ 0.5e. The density scaling is applied in Fig. 9. Density scaling with ␥ = 4.9 works better than either ␥ = 6.0 or ␥ = 4.0, but the data collapse is clearly inferior to what is seen in Figs. 5 and 6 ͑note the smaller dynamic range in Fig. 9͒. This confirms that density scaling is indeed a property of strongly correlating liquids, consistent with the experimental finding that, e.g., hydrogen-bonded liquids do not obey density scaling.

V. CONCLUSIONS
We have provided direct evidence for an underlying IPL scale invariance in model molecular van der Waals liquids, leading to three experimentally testable predictions for these liquids: ͑i͒ the density scaling exponent ␥ can be identified from the equilibrium fluctuations via the fluctuationdissipation theorem determining the frequency-dependent linear thermoviscoelastic response functions ͓2,23͔; ͑ii͒ different state points with same ␥ / T-same values of D ‫ء‬ and ␣ -are "isomorphic" ͓24͔, i.e., have same scaled structure. ͑iii͒ Isomorphic state points must have the same relaxation spectra. This is consistent with recent experimental findings for molecular glass formers and amorphous polymers ͑excluding hydrogen bonds͒: varying pressure and temperature, state points with the same relaxation time have the same shape of the dielectric ␣-relaxation function ͓25,26͔. Finally, the hidden scale invariance explanation of density scaling explains why one observes density scaling of the dynamics, without the same scaling applying to the pressure ͓11,27͔.
The hidden scale invariance is a consequence of the strong virial/potential energy correlations ͓3,4͔ characterizing van der Waals liquids, as well as most or all metallic liquids. Thus the results from computer simulations of two model liquids presented here point to a new physics of these large classes of liquids, which is simpler than previously believed.

ACKNOWLEDGMENTS
Useful discussions with Kristine Niss are gratefully acknowledged. This work was supported by the Danish National Research Foundation's ͑DNRF͒ center for viscous liquid dynamics "Glass and Time." FIG. 8. ͑Color online͒ Results from equilibrium molecular dynamics simulations of 512 asymmetric dumbbell molecules ͓9͔ 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 coefficients, R ϵ͗⌬W⌬U͘ / ͱ ͗⌬W 2 ͗͘⌬U 2 ͘.