Thermodynamic scaling of diffusion in supercooled Lennard-Jones liquids
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p Thermodynamic Scaling of Diffusion in Supercooled Lennard-Jones Liquids
D. Coslovich
1, 2 and C. M. Roland Dipartimento di Fisica Teorica, Universit`a di Trieste – Strada Costiera 11, 34100 Trieste, Italy CNR-INFM Democritos National Simulation Center – Via Beirut 2-4, 34014 Trieste, Italy Naval Research Laboratory, Code 6120, Washington, DC 20375-5342, USA (Dated: October 31, 2018)The manner in which the intermolecular potential u ( r ) governs structural relaxation in liquids isa long standing problem in condensed matter physics. Herein we show that diffusion coefficientsfor simulated Lennard-Jones m –6 liquids (8 ≤ m ≤
36) in normal and moderately supercooledstates are a unique function of the variable ρ γ /T , where ρ is density and T is temperature. Thescaling exponent γ is a material specific constant whose magnitude is related to the steepness of therepulsive part of u ( r ), evaluated around the distance of closest approach between particles probedin the supercooled regime. Approximations of u ( r ) in terms of inverse power laws are also discussed. PACS numbers: 61.43.Fs, 61.20.Lc, 64.70.Pf, 61.20.Ja
Establishing a quantitative connection between the re-laxation properties of a liquid and the interactions amongits constituent molecules is the sine qua non for funda-mental understanding and prediction of the dynamicalproperties. The supercooled regime is of particular in-terest, since both intermolecular forces and steric con-straints (excluded volume) exert significant effects on thedynamics. This makes temperature, pressure, and vol-ume essential experimental variables to characterize therelaxation properties. One successful approach to at leastcategorize dynamic properties of supercooled liquids andpolymers is by expressing them as a function of the ratioof mass density, ρ , to temperature, T , with the formerraised to a material specific constant γ , viz. x = ℑ ( ρ γ /T ) (1)where x is the dynamic quantity under consideration,such as the structural relaxation time τ , the viscosity η ,or the diffusion coefficient D , and ℑ is a function. Thisscaling was first applied to a Lennard-Jones (LJ) fluid,with γ = 4 yielding approximate master curves of thereduced “excess” viscosity for different thermodynamicconditions [1]. More recently Eq. (1) has been shown tosuperpose relaxation times measured by neutron scatter-ing [2], light scattering [3], viscosity [4], and dielectricspectroscopy [5, 6, 7, 8, 9] for a broad range of mate-rials, including polymer blends and ionic liquids. Thescaling exponent γ , which varies in the range from 0.13to 8.5 [10], is a measure of the contribution of density (orvolume) to the dynamics, relative to that due to tem-perature. The only breakdown of the scaling is observedfor hydrogen-bonded liquids, in which the concentrationof H-bonds changes with T and P , causing τ to deviatefrom Eq. (1) [4].The function ℑ in Eq. (1) is unknown a priori . Itsform can be derived from entropy models for the glasstransition, leading to an exponential dependence of log τ on ρ γ /T [11, 12]. Another interpretation of the scaling isthat the supercooled dynamics is governed by activated processes with an effective activation energy E ( ρ, T ) [13],in which the ρ -dependence of E ( ρ, T ) can be factoredand expressed in terms of a power law of ρ . The powerlaw scaling arose from the idea that the intermolecularpotential for liquids can be approximated as a repulsiveinverse power law (IPL), with the weaker attractive forcestreated as a spatially-uniform background term [14, 15,16] u ( r ) ∼ r − m + const (2)where r is the intermolecular distance. In the case ofan IPL, in fact, all reduced dynamical quantities [17]can be cast in the form of Eq. (1) with γ = m/
3, i.e.,the thermodynamic scaling is strictly obeyed. For in-stance, this applies to the reduced diffusion coefficient D ∗ ∼ ( ρ / T − / ) D ∼ ( T /m − / ) D . A similar reductionof D by macroscopic variables ( ρ and T ) has also beenemployed in entropy scaling laws of diffusion [18].The IPL approximation emphasizes the dominant roleof the short-range repulsive interactions for local proper-ties such as structural relaxation. Various groups haveexplored through numerical simulations the relationshipof the steepness of the repulsive potential to propertiessuch as the equation of state [19, 20, 21], longitudinalwave transmission [22], vibrational spectrum [23], liq-uid [24] and gaseous [25] transport, the correlation be-tween fluctuations of energy and pressure [26], and thefragility [27, 28]. Recently two simulations have appearedin which Eq. (1) was used to superpose dynamical datafor polymer chains described using an LJ m –6 potentialwith m = 12 and an added term for the intrachain in-teractions. The results appear contradictory: Tsolou etal. [29] obtained a scaling exponent of γ = 2 . γ = 6 when attrac-tion were included in the simulation and γ = 12 when theywere omitted. Thus, the scaling exponent γ is either lessthan [29] or greater than [30] m/ m – n liquids, in which the repulsive exponent m was systemat-ically varied. Our models are binary mixtures composedof N = 500 particles enclosed in a cubic box with peri-odic boundary conditions and interacting with a LJ m – n potential u αβ ( r ) = 4 ǫ αβ [( σ αβ /r ) m − ( σ αβ /r ) n ] (3)where α, β = 1 , n = 6, as in the standard LJ poten-tial, and varied m = 8 , , ,
36. The potential u αβ ( r )was smoothed at r c = 2 . σ αβ using the cutoff scheme ofStoddard and Ford [31]. Reduced LJ units are used, as-suming σ , ǫ and p m σ /ǫ as units of distance,energy and time respectively. The mixture on which wefocus is an additive, equimolar mixture with size ratio λ = σ /σ = 0 .
64, equal masses m = m = 1 . ǫ αβ = 1 .
0. The choice m = 12 cor-responds to the AMLJ-0.64 mixture studied in [32, 33].The samples were quenched isobarically at different pres-sures P = 5 , ,
20 by coupling the system to Berendsenthermostat and barostat during equilibration (see [32]for details), and performing the production runs in theNVE ensemble using the Velocity-Verlet algorithm. Thetimestep δt was varied according to the repulsive expo-nent, ranging from 0.001 ( m = 36) to 0.004 ( m = 8) athigh T , and from 0.003 ( m = 36) to 0.008 ( m = 8) atlow T . The equilibration criteria were similar to the onesused in [32].The effectiveness of the thermodynamic scaling for LJ m − m . For each m , reduced diffu-sion coefficients D ∗ = ( ρ / T − / ) D were gathered alongdifferent isobaric paths ( P = 5 , ,
20) and the materialspecific scaling exponent γ was obtained by maximizingthe overlap between different sets of data, plotted as afunction of ρ γ /T . Repeating the analysis for D , insteadof D ∗ , yields very similar values of γ , but the quality ofthe scaling for D ∗ is slightly superior. Our data spanroughly 5 decades of variation of D , over about two ofwhich the temperature is lower than the so-called onsettemperature T O [34], where non-exponential relaxationtypical of the supercooled regime first becomes appar-ent upon cooling the liquid. Analyzing the variation ofthe scaling exponent in our models, we find that γ in-creases with increasing m , but its actual value is system-atically larger than m/
3. For instance, in the case m = 12we obtain γ = 5 .
0, a value which we also found to pro-vide scaling of D ∗ for other supercooled Lennard-Jones( m = 12) mixtures, such as the AMLJ-0.76 mixture intro-duced in [32] and the mixture of Kob and Andersen [35].The origin of the discrepancy between γ and m/ -5 -4 -3 -2 -1
0 5 10 15 20 25 30 35 40 D * ρ γ /Tm=8 γ =3.5 P=5P=10P=2010 -5 -4 -3 -2 -1
0 5 10 15 20 25 30 D * ρ γ /Tm=12 γ =5.0 P=5P=10P=2010 -5 -4 -3 -2 -1
0 20 40 60 80 100 120 140 D * ρ γ /Tm=24 γ =9.1 P=5P=10P=2010 -5 -4 -3 -2 -1
0 200 400 600 800 1000 1200 D * ρ γ /Tm=36 γ =13.4 P=5P=10P=20 FIG. 1: (color online). Reduced diffusion coefficients D ∗ as afunction of ρ γ /T for different values of the repulsive exponent m at different pressures: P = 5 (squares), P = 10 (circles), and P = 20 (triangles). From top to bottom: m = 36 ( γ = 13 . m = 24 ( γ = 9 . m = 12 ( γ = 5 . m = 8 ( γ = 3 . γ is ± . ± . m = 36). ticle distances, in which u ( r ) ∼ r − m , is not dynamicallyaccessible in normal simulation conditions. The presenceof the fixed attractive term in the potential (Eq. (3))gives rise to an effective IPL which is steeper in the re-gion of r close to the minimum than in the r → m = 24.The lower panel of Fig. 2 shows a fit of the pair poten-tial u ( r ) to an IPL (Eq. (2)) performed in the range g (r) r m=24P=10, T=0.84P=10, T=1.00P=10, T=1.50 0 1 2 3 4 5 6 0.9 1.0 1.1 1.2 1.3 1.4 1.5 g (r) r m=24P=5, T=0.75P=10, T=0.84P=20, T=1.05-4-2 0 2 4 6 8 10 0.9 1.0 1.1 1.2 1.3 1.4 1.5 u (r) r LJ m=24IPL m=27.5 FIG. 2: (color online). Top panel: radial distribution func-tions between large particles g ( r ) at P = 10 for T . T O : T = 1 .
20 (dotted), T = 1 .
00 (dashed), and T = 0 .
84 (solid).Middle panel: g ( r ) at the lowest equilibrated T : T = 0 . P = 5 (dotted), T = 0 .
84 at P = 10 (dashed), and T = 1 . P = 20 (solid). Bottom panel: pair potential u ( r ) (solid)and fitted IPL (dotted) in the range [0 .
95 : 1 . [ r : r ], with r = 0 .
95 and r = 1 .
01. The value m = 27 . m = 24 and is in very good agreement with the value ex-pected from the dynamical scaling (3 γ = 27 . ± . r : r ] corresponds to typical distances of closestapproach between particles probed within our simulationconditions, as it can be seen by inspection of the radialdistribution functions g ( r ) (see upper panels of Fig. 2).Extending the range for the fit up to r = 1 .
06, whichis close to the average position of the first peak in the g ( r ), yields a larger value m = 28 .
8, revealing how γ isdictated by the portion of r around the distance of closestapproach in the supercooled regime.To proceed in a more systematic way, we consideredall α − β pairs (1-1, 1-2 and 2-2) in the potential u αβ ( r )and performed a simultaneous fit to the following IPL u αβ ( r ) = ǫ ( σ αβ /r ) m + k . (4)The range for fitting was defined by two conventionaldistances determined from the radial distribution func- TABLE I: Parameters of IPL approximations for u αβ ( r ). Theeffective exponent m is obtained from fitting to Eq. (4),whereas ǫ , k , and x are the optimal values for Eq. (5). m γ m x ǫ k − . − . − . − . tions g αβ ( r ): the distance of closest approach betweenparticles, r , (i.e., the value of r for which the g αβ ( r )first becomes non-zero) and the position correspond-ing to half of the height of the first peak, r , (i.e., g αβ ( r ) = g αβ ( r m ) / r m is the position of thefirst peak and r < r < r m ). These quantities depend onthe thermodynamic state under consideration, but theirvariation with P and T is mild within our simulation con-ditions [38]. Our interest being the supercooled regime,we simply consider the interval [ r : r ] obtained fromthe low- T behavior of the g αβ ( r ). For each α − β pair weused the corresponding range [ r : r ] for fitting. In gen-eral, the fitted values of m are in good agreement with3 γ (see Table I) for all values of m . Thus, the scalingexponent can be reasonably accounted for in terms of anIPL approximation of the pair potential, provided that asensible choice of the relevant range of distances is made.The above procedure suggests that a model of soft-spheres (SS) with m = 3 γ should provide a good refer-ence system for the LJ m –6 mixtures. To this aim, weapproximate Eq. (3) with v αβ ( r ) = ( ǫ ( σ αβ /r ) m + k r < xσ αβ u αβ ( r ) r ≥ xσ αβ (5)where m , ǫ , and k are expressed in terms of x by requir-ing continuity of 0th, 1th, and 2th derivatives of v αβ ( r )at r = xσ αβ . The value of x is then fixed by requiringthat 3 γ = m ( x ) = ( m /x m +1 − n /x n +1 ) / ( m/x m +1 − n/x n +1 ). The parameters defining the reference SS mod-els for all values of m are reported in Table I. We checkedthat the distance xσ αβ always lies in the range [ r : r ]defined above. Diffusivity data for the LJ 12–6 mixtureare compared in Fig. 3 to those of the corresponding ref-erence SS mixture along two isochores ( ρ = 1 . ρ = 1 . T by the LJ system (at constant P ). The trend of D ( T ) forthe reference system closely follows the one for the fullLJ system. As expected, the SS mixture has a larger dif-fusion coefficient for a given thermodynamic state. Thecontribution to D due to the attractive part of the po-tential could also be explicitly included using a WCA-likesplitting of v αβ ( r ) [36]. For the present purposes, how-ever, it is more useful to note that a simple rescaling of -5 -4 -3 -2 -1
0 1 2 3 4 5 6 D ρ =1.5LJ ρ =1.7SS ρ =1.5SS ρ =1.7 10 -4 -3 -2 -1
0 5 10 15 20 25 30 D * ρ /T FIG. 3: (color online). Arrhenius plot of diffusion coeffi-cient D for the LJ 12–6 mixture and the reference SS mix-ture ( m = 15 . ǫ = 1 .
74) along two isochores: ρ = 1 . ρ = 1 .