Density-scaling exponents and virial potential-energy correlation coefficients for the (2n,n) Lennard-Jones system Ida M. Friisberg,∗ Lorenzo Costigliola,† and Jeppe C. Dyre‡ “Glass and Time”, IMFUFA, Department of Science and Environment, Roskilde University, P.O. Box 260, DK-4000 Roskilde, Denmark (Dated: March 6, 2017) This paper investigates the relation between the density-scaling exponent γ and the virial potential-energy correlation coefficient R at several thermodynamic state points in three dimen- sions for the generalized (2n,n) Lennard-Jones (LJ) system for n = 4,9,12,18, as well as for the standard n = 6 LJ system in two, three, and four dimensions. The state points studied include many low-density states at which the virial potential-energy correlations are not strong. For these 7 statepointswefindtheroughlylinearrelationγ ∼=3nR/dinddimensions. Thisresultisdiscussed 1 in light of the approximate “extended inverse power law” description of generalized LJ potentials 0 [N. P. Bailey et al., J. Chem. Phys. 129, 184508 (2008)]. In the plot of γ versus R there is in all 2 cases a transition around R ≈ 0.9, above which γ starts to decrease as R approaches unity. This r is consistent with the fact that γ → 2n/d for R → 1, a limit that is approached at high densities a and/or high temperatures at which the repulsive r−2n term dominates the physics. M 3 ] t f o s . t a m - d n o c [ 2 v 0 7 4 5 0 . 1 0 7 1 : v i X r a 2 I. INTRODUCTION Inthepastdecadeaclassofsystemshasbeenidentifiedthatiswelldescribedbytheisomorphtheory1–8,sometimes referredtoasR(Roskilde)simplesystems9–18. ThisclassisbelievedtoincludemostvanderWaalsbondedandmetallic liquids and solids, as well as most weakly dipolar or ionic systems1–8. The class does not include systems with strong directional bonds like hydrogen or covalently bonded systems8. To determine whether the isomorph theory holds for a given system one calculates the strength of the virial potential-energy correlations at the state points in question in NVT simulations, i.e., in the canonical ensemble1,2. Isomorphsarecurvesofinvariantstructureanddynamicsinthethermodynamicphasediagram,seeSec. II.Asystem hasisomorphsifandonlyifithasstrongvirialpotential-energycorrelations4. Virialpotential-energycorrelationsare quantified by the state-point dependent Pearson correlation coefficient R defined as (cid:104)∆W∆U(cid:105) R(ρ,T)= . (1) (cid:113) (cid:104)(∆W)2(cid:105)(cid:104)(∆U)2(cid:105) InEq.(1)ρisthe(number)density(ρ=N/V withN particlesinvolumeV),T isthetemperature,W andU arethe virialandpotentialenergy,thebracketsdenoteNVT averages,and∆istheinstantaneousdeviationfromequilibrium mean values. The isomorph theory is exact only for systems with a potential that is a constant plus an Euler-homogeneous function, in which case the correlation coefficient R in Eq. (1) is unity for all ρ and T4. An example is the inverse- power-law (IPL) pair-potential systems. For non-Euler-homogeneous potentials, which includes all realistic systems, thequantityR islessthanunity. Inthatcasetheapplicabilityoftheisomorphtheoryisrestrictedtoacertainregion of the phase diagram, usually the condensed-phase region encompassing the solid and “ordinary” liquid phase. A system is defined to be R simple whenever R > 0.91, because this ensures that isomorphs exist to a good approximation4. TheR=0.9thresholdissomewhatarbitrary,however. Thepresentpaperinvestigateswhathappens when R falls below this threshold. From the isomorph theory there is little help, because it generally breaks down when correlations are no longer strong, i.e., when the correlation coefficient goes significantly below 0.92,4,7. In 2014, however, an interesting paper appeared by Prasad and Chakravarty13 establishing that Rosenfeld’s excess entropy scalingaswellasdensity-scaling19, whichmaybothbederivedfromtheisomorphtheory4, neednotbreakdowneven if the virial potential-energy correlations are weak. Reference 13 studied the transition to simple liquid behavior in computer simulations of modified water models. We have taken inspiration from this work to systematically study the region of the phase diagram where virial potential-energy correlations are not strong for a class of models that – in contrast to water models – do have sizable regions of strong correlations in the thermodynamic phase diagram. This is done by studying Lennard-Jones type systems at lower densities than have previously been in focus. This paper presents a study of the LJ(2n,n) class of potentials defined as single-component systems with LJ-type radially symmetric pair potentials with a repulsive term proportional to r−2n and an attractive term proportional to r−n (r being the interparticle distance). Results are presented for the cases n = 4,9,12,18 in three dimensions, as well as for the standard LJ potential (n = 6) in two, three, and four spatial dimensions (Sec. III). Before doing this the isomorph theory is briefly reviewed (Sec. II). Section IV summarizes the paper’s findings. II. THEORETICAL/COMPUTATIONAL A. Aspects of the isomorph theory R simple systems – previously termed “strongly correlating” which, however, led to confusion with strongly cor- related quantum systems – were defined in 2009 by reference to strong virial potential-energy correlations1. This section presents the alternative definition of the same class of systems given in 201414 that refers to their “hidden scaleinvariance”7. Theoriginalformulationoftheisomorphtheory4 isrecoveredviaafirst-orderTaylorexpansion14. Consider an N-particle system in d spatial dimensions. The system is defined to be R simple if for any two configurations corresponding to the same density that obey U(R ) < U(R ), this ordering is preserved after a a b uniform scaling of the two configurations14 (R is the dN dimensional vector describing a configuration of N particles in d spatial dimensions). Formally, this condition is written14 U(R )<U(R )⇒U(λR )<U(λR ) . (2) a b a b 3 In computer simulations periodic boundary conditions are applied, and it is understood that the box size scales with R. Equation (2) need not be obeyed for all configurations; as long as it applies for most of the physically relevant configurations,thepredictionsofisomorphtheoryareobeyedtoagoodapproximation. InEq.(2)thescalingfactorλ can be any positive real number, and the equation holds for scaling ”both ways”. Because of this an R simple system obeys14 U(R )=U(R )⇒U(λR )=U(λR ) . (3) a b a b Thus if two configurations have the same potential energy, their scaled potential energies are also identical. R simple systems are characterized by strong correlations between virial and potential-energy fluctuations in the NVT ensemble. This can be derived from Eq. (3) in the following way. For two configurations, R and R , at the a b same density with the same potential energy, Eq. (3) states that U(λR ) = U(λR ). Taking the derivative of this a b identity with respect to λ one obtains R ·∇U(λR )=R ·∇U(λR ). (4) a a b b Using the definition of virial W(R)≡R·∇U(R)/d20, one gets for λ=1 W(R )=W(R ). (5) a b Thus U(R) determines W(R), implying perfect correlation. This applies for systems that satisfy Eq. (2) or, equiva- lentlyEq.(3),forallconfigurations;itholdsforrealisticRsimplesystemstoagoodapproximation. Asaconsequence, for such systems the correlation coefficient R of Eq. (1) is close to unity, but not exactly unity. As mentioned, the threshold defining R simple systems has usually been taken to be R=0.91. Two configurations of density ρ and ρ , which scale uniformly into one another, i.e., R = λR for some λ, are 1 2 2 1 related through the following equality (where d is the dimension): ρ1/dR =ρ1/dR ≡R˜ . (6) 1 1 2 2 The last equality defines the “reduced” (dimensionless) configuration vector R˜. Henceforth reduced quantities are marked by a tilde. The reduced version of a physical quantity is obtained by using the “macroscopic” unit system in √ whichthelengthunitisρ−1/d,theenergyunitisk T,andthetimeunitisρ−1/d/ mk T (mistheparticlemass)21. B B The entropy of any physical system can be written as an ideal-gas term plus an excess term S . The ideal-gas ex term is the entropy of an ideal gas at the density and temperature in question, the excess term is the contribution to entropy from the interactions between the particles. We define the microscopic excess entropy function, S (R), to ex be the thermodynamic excess entropy14,22 of the potential-energy surface the configuration belongs to, i.e., S (R)≡S (ρ,U(R)) . (7) ex ex The function on the right-hand side, S (ρ,U), is the thermodynamic excess entropy of the state point with density ex ρ and average potential energy U(R). In other words, the microscopic excess entropy of a configuration R is defined as the excess entropy of a thermodynamic equilibrium system of same density with average energy precisely equal to U(R). Inverting Eq. (7) we get U(R)=U(ρ,S (R)) (8) ex in which the function on the right-hand side, U(ρ,S ), is the thermodynamic average potential energy of the state ex point with density ρ and excess entropy S . ex Equations (7) and (8) apply for any system22, but they are particularly significant for R simple systems for which they imply that the configurational adiabats are curves of invariant structure and dynamics. To prove this, we first note that, as shown in Ref. 14, Eq. (2) implies that the excess-entropy function is invariant under uniform scaling, i.e., it only depends on the reduced coordinate vector R˜: 4 S (R)=S (R˜). (9) ex ex The relation for the potential-energy function Eq. (8) consequently becomes U(R)=U(ρ,S (R˜)). (10) ex Isomorphs are defined as the configurational adiabats, i.e., curves along which the excess entropy is constant, in the region of phase diagram where the system is R simple4. To demonstrate invariance of structure and dynamics along the isomorphs we show that Newton’s second law in reduced units is invariant along an isomorph. In complete generality, this law is in reduced units (note that the particle mass is absorbed into the reduced time) d2R˜ F (R) = ≡F˜(R). (11) dt˜2 ρ1/dk T B For R simple systems Eq. (10) implies (cid:18) (cid:19) ∂U F =−∇U =− ∇S (R˜). (12) ∂S ex ex ρ Using ∇=ρ1/d∇˜ and (∂U/∂S ) =T, the above expression becomes ex ρ F =−ρ1/dT∇˜S (R˜) (13) ex or ∇˜S (R˜) F˜ =− ex . (14) k B This expression reveals that for R simple systems the reduced force vector F˜ is a function of the reduced coordinate vector, implying that Eq. (11) is invariant along an isomorph. Thus the dynamics is isomorph invariant in reduced units, which implies that the reduced-unit structure is also isomorph invariant14. B. The density-scaling exponent When strong correlations between virial and potential energy are present, a constant of proportionality between the instantaneous fluctuations of these two quantities’ deviation from their equilibrium values, γ, can be introduced via ∆W (t)(cid:39)γ∆U(t) . (15) The correlation coefficient associated with this linear regression is that given in Eq. (1). The exact definition of γ is the following4: (cid:104)∆W∆U(cid:105) γ = . (16) (cid:104)(∆U)2(cid:105) By applying a standard fluctuation relation and the volume-temperature Maxwell relation one can show4 that (cid:18) (cid:19) ∂lnT γ = . (17) ∂lnρ Sex 5 The number γ is termed the density-scaling exponent4,19. It determines the configurational adiabats, which as mentionedareisomorphsintheRsimpleregionofthephasediagram. Ifvariationsofγ areinsignificant,thesecurves areviaEq.(17)givenbyργ/T =Const. Inthegeneralcase,Eq.(17)canbeusedtotraceoutisomorphsstep-by-step by repeatedly changing density by typically an amount of order 1%, calculating the temperature change via Eq. (17), recalculating γ from Eq. (16) at the new state point, etc. In early publications1,2 we identified the constant of proportionality γ of Eq. (15) by the following symmetric fluctuation expression, which was later renamed γ to distinguish it from γ4: 2 (cid:115) (cid:104)(∆W)2(cid:105) γ = . (18) 2 (cid:104)(∆U)2(cid:105) Whenever the correlation coefficient R is close to unity, one has γ ∼=γ since the following applies 2 γ =γ R. (19) 2 Now that the main ingredients of the isomorph theory have been introduced, we proceed to present the simulation results. III. RESULTS AND DISCUSSION A. Generalized Lennard-Jones pair potentials in three dimensions If r is the interparticle distance, the generalized LJ pair potential is defined as follows ε (cid:104) (cid:16)σ(cid:17)m (cid:16)σ(cid:17)n(cid:105) vLJ (r)= n −m . (20) m,n m−n r r Here m and n are positive integers, and in order to ensure thermodynamic stability it is assumed that m > n. The constantsσ andεdefinethepotential’slengthandenergyscales, respectively, andthenormalizationusedinEq.(20) ensures that the minimum pair potential energy is −ε which is obtained at r = σ. Note that the normalization is differentfromthatusuallyemployedforthestandard12-6LJpairpotentialparametrizedas4ε(cid:2)(r/σ)−12−(r/σ)−6(cid:3). The aim of our study is to investigate whether any relation between the correlation coefficient and the density- scalingexponentcanbedetermined. Thisinvolvessimulatingseveralstatepointsforwhichthevirialpotential-energy correlations are not strong. With this goal in mind, a particular case of the generalized LJ potential was simulated, the case where m=2n: (cid:20)(cid:16)σ(cid:17)2n (cid:16)σ(cid:17)n(cid:21) vLJ (r)=ε −2 . (21) 2n,n r r Thefollowingfourinstancesofthispairpotential,henceforthdenotedbyLJ(2n,n),weresimulatedinthreedimensions: n = 4,n = 9,n = 12 and n = 18 (the n = 6 case corresponding to the standard LJ potential is considered later). Figure 1 shows that these four pair potentials are quite different as functions of the pairwise distance r. The simulations were performed in the NVT ensemble, i.e., for a constant number of particles (N = 864 for n = 4,9,12 and N = 4000 for n = 18) at constant temperature T and constant volume V. The time step was 0.001 in LJ units and the simulations were performed with periodic boundary conditions and a standard shifted potential cut-off at 2.5σ. Using an FCC crystal as starting configuration, each system was equilibrated for 2·106 time steps before the collection of data began. After equilibration, the simulation ran for 50·106 time steps during which data were collected. For each of the four systems, six densities were considered, ρ = 0.25,0.50,0.75,1.00,1.25,1.50 (LJ units); the temperature was varied within the range 0.25 < T < 5.00 (LJ units). At each state point the correlation coefficient and the density-scaling exponent were calculated from Eq. (1) and Eq. (16), respectively. What to expect for the behavior of R when density is lowered? To answer this question we refer to the isomorph theory, which works well for the LJ and related pair potentials in the ordinary liquid phase. According to isomorph theory, if a reference state point (ρ ,T ) and another state point (ρ,T) are on the same isomorph, the ratio between 0 0 T and T defines what may be termed the isomorph shape function h(ρ,S ) via T/T = h(ρ,S )23. Thus each 0 ex 0 ex isomorph is mapped out in the phase diagram by 6 1.5 (a) LJ(8,4) LJ(18,9) 1.0 LJ(24,12) LJ(36,18) 0.5 ) r (n 0.0 LJ2n, v -0.5 -1.0 -1.5 0.5 1.0 1.5 2.0 2.5 3.0 r FIG. 1. The simulated four generalized Lennard-Jones (LJ) pair potentials of Eq. (21) plotted as a function of the pairwise distance r (distance and energy are given in LJ units): n = 4 (black), n = 9 (red), n = 12 (green), and n = 18 (blue). The potentials are quite different. h(ρ,S ) ex = Const. (22) T The original isomorph theory predicted the function h(ρ,S ) to be independent of S , whereas the more correct ex ex theoryfrom2014predictsthathmayvaryslightlyfromisomorphtoisomorph14. Inbothcases,however,theanalytical formofthedensitydependenceisthesame. FortheLJ(2n,n)pairpotentialitcanbeshownthattheisomorphshape function takes the following form in d dimensions14,23–25 (cid:18)d (cid:19)(cid:18) ρ (cid:19)2n/d (cid:18)d (cid:19)(cid:18) ρ (cid:19)n/d h(ρ,S )= γ −1 − γ −2 . (23) ex n 0 ρ n 0 ρ 0 0 Here (ρ ,T ) is the isomorph’s reference state point, at which γ is the density-scaling exponent – the S (isomorph) 0 0 0 ex dependence of the function h is contained in γ . Combining Eqs. (22) and (23) one can map out an isomorph from 0 information obtained by computer simulations at the reference state point, i.e., with no need to use Eq. (17) in a tedious step-by-step process. Note that Eq. (22) and the definition of the density-scaling exponent Eq. (17) implies26 dlnh γ = . (24) dlnρ Thuswhenthedensity-scalingexponentisknownatasinglereferencestatepoint,thetheoryviaEq.(23)andEq.(24) predicts how γ varies with density along the isomorph in question. InEq.(23)therearetwoterms,onepositiveandonenegative. Athighdensitiesthepositivetermdominates. Upon lowering the density a point is reached where h(ρ,S ) changes sign. Below this density the theory implies negative ex isomorph temperatures, compare Eq. (22), which shows that the isomorph theory must break down here. Since the theory works well whenever there are strong virial potential-energy correlations, one concludes that the correlation coefficientmustdecreaseuponloweringdensityalongaconfigurationaladiabat,i.e.,forgeneralizedLJsystemsinany dimension the isomorph theory predicts its own breakdown at sufficiently low densities. To summarize, moving along a given isomorph it is possible to define a density below which the isomorph theory does not hold. The condition that the shape function is non-negative sets the following limit for the isomorph theory to work: (cid:18)dγ −2n(cid:19)d/n ρ > 0 ρ . (25) dγ −n 0 0 7 LJ(8,4) LJ(18,9) 5.0 F F F F F F 5.0 F F F F F C F F F F F F F F F F F C 4.0 F F F F F F 4.0 F F F F F C F F F F F F F F F F F C 3.0 F F F F F F 3.0 F F F F F C T F F F F F F T F F F F F C 2.0 F F F F F F 2.0 F F F F F C L+G L+G L+G L L L F F F F C C F F F F C C 1.0 L+G L+G L+G L+G L L 1.0 F F F F C C F F F F C C L+G L+G L+G L+G L+G C L+G L+G L+G L+G C C C+G C+G C+G C+G C+G C C+G C+G C+G C+G C C 0.0 0.0 0.25 0.50 0.75 1.00 1.25 1.50 0.25 0.50 0.75 1.00 1.25 1.50 ρ ρ LJ(24,12) LJ(36,18) 5.0 F F F F F C 5.0 F F F F C C F F F F F C F F F F C C 4.0 F F F F C C 4.0 F F F F C C F F F F C C F F F F C C 3.0 F F F F C C 3.0 F F F F C C T F F F F C C T F F F F C C 2.0 F F F F C C 2.0 F F F F C C F F F F C C F F F F C C 1.0 F F F F C C 1.0 F F F F C C F F F F C C F F F F C C L+G L+G L+G L C C F F F F C C C+G C+G C+G C+G C C C+G C+G C+G C+G C+G C 0.0 0.0 0.25 0.50 0.75 1.00 1.25 1.50 0.25 0.50 0.75 1.00 1.25 1.50 ρ ρ FIG.2. StatepointssimulatedforthepairpotentialsLJ(8,4),LJ(18,9),LJ(24,12),LJ(36,18)definedbyEq.(21). Ateachstate point the virial potential-energy correlation coefficient (Eq. (1)) and the density-scaling exponent (Eq. (16)) were evaluated. Most state points are at supercritical temperatures (at which the liquid and gas phases merge) and marked by an F (F means fluid,Lmeansliquid,Cmeanscrystalline,L+Gmeansastatepointofcoexistingliquidandgasphases,etc). Thephaseswere identified by visual inspection of selected configurations. WehaveseenthatRdecreasesasdensityisloweredalongaconfigurationaladiabat(isomorph). Moreover,according to Eq. (16) the scaling exponent must approach zero if R → 0. Since R = 0 implies γ = 0, one expects that γ also decreasesuponloweringthedensity. Ifdensityisincreased,ontheotherhand,thepositivetermofEq.(23)dominates and the density-scaling exponent eventually approaches 2n/d due to the dominance of the repulsive r−2n term23–25. At the same time R tends to unity. After considering these two limits, we turn to the simulation results. Figure 2 shows the state points studied, marking in each case whether it is a solid or liquid state point, etc. Most state points are supercritical in which case we marked them by an F indicating “fluid”. Figure 3 plots γ versus R alongtheisochoresfortheLJ(2n,n)potentialswithn=4,9,12,18. Ourfindingsareconsistentwiththetheoretically predictedresultthatasR→1, γ approachesthehigh-densitylimitingvalues8/3, 6, 8, and12, respectively, although wenevercomereallyclosetothesevaluesthatapplywhenevertheattractivetermofthepotentialmaybecompletely ignored. The density-scaling exponent changes behavior and starts to decrease with increasing R in the region where R is around 0.9 which, as mentioned earlier, is the threshold usually used to define R simple systems. In the opposite limit, R→0, we find γ →0, as expected. Our data show that the change of behavior in γ versus R takes place at the boundary of the region in which the systemisRsimple. AmoreextensivestudyisneededtodeterminewhetherthisholdsalsoforotherRsimplesystems. If confirmed, the relation between the change in the behavior of the density-scaling exponent and the value of the correlation coefficient could be useful in practice, because R is not accessible in experiment while γ is27. A change in the density or temperature dependence of γ might therefore be used for identifying in which regions of the phase diagramaliquidisexpectedtoobeyisomorphtheorypredictionslikeisochronalsuperposition,excessentropyscaling, etc7. Figure4summarizesthedataofFig.3. Thereisaroughlylinearbehaviorbetweenγ andR,butitdoesnotextend above R≈0.9 where γ as mentioned starts to decrease towards the value 2n/3 predicted at high density (where the r−2n repulsive term dominates and the virial potential-energy correlations become virtually perfect). For R > 0.9 the system is R simple and the density variation of γ is well described by Eqs. (23) and (24). For this reason we 8 4.0 (a) ρ = 0.25 (b) ρ = 0.25 ρ = 0.25 ρ = 0.25 ρ = 0.25 10.0 ρ = 0.25 3.0 3.5 ρ = 0.50 10.0 ρ = 0.50 γ ρ = 0.50 γ ρ = 0.50 ρ = 0.50 7.5 7.5 ρ = 0.50 2.0 ρ = 0.75 ρ = 0.75 3.0 ρ = 0.75 0.8 0.9 1.0 ρ = 0.75 γ 0.8 0.9 1.0 ρ = 0.75 γ 5.0 R ρ = 0.75 R ρ = 1.00 ρ = 1.00 1.0 ρ = 1.00 ρ = 1.00 ρ = 1.00 ρ = 1.00 ρ = 1.25 2.5 ρ = 1.25 0.0 ρ = 1.25 ρ = 1.25 LJ(8, 4) LJ(18,9) ρ = 1.25 ρ = 1.50 ρ = 1.50 0.0 -1.0 ρ = 1.50 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 R R (c) ρ = 0.25 (d) ρ = 0.25 15 12 ρ = 0.25 20 ρ = 0.25 ρ = 0.25 20 ρ = 0.50 11 ρ = 0.50 ρ = 0.50 γ10 ρ = 0.50 15 γ ρ = 0.75 ρ = 0.50 15 ρ = 0.75 10 9 ρ = 0.75 ρ = 1.00 0.8 0.9 1.0 ρ = 0.75 0.8 0.9 1.0 ρ = 1.00 γ R ρ = 0.75 γ10 R ρ = 1.25 ρ = 1.00 ρ = 1.25 5 ρ = 1.00 ρ = 1.50 ρ = 1.25 5 ρ = 1.25 ρ = 1.50 LJ(24,12) LJ(36,18) 0 0 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 R R FIG. 3. The density-scaling exponent, γ, versus the correlation coefficient, R, for generalized LJ potentials: (a) LJ(8,4); (b) LJ(18,9); (c) LJ(24,12); (d) LJ(36,18). Full squares mark crystalline states, full circles are liquid, fluid, or gas states – open squares are coexisting crystal-gas states; open circles are coexisting liquid-gas states. The simulated state points are given in Fig. 2. The inset in each figure focuses on the crossover region. Along each isochore the temperature was varied in the range 0.25 < T < 5.00. In (b), (c), and (d) it is possible to distinguish two different behaviors of the variation of γ with R, which reflects liquid contra crystal phases (the latter being the steepest). The state points with negative or close to zero correlation coefficient are all in the gas-liquid or gas-solid coexistence regions. 1.5 20 (a) (b) LJ(8,4) LJ(8,4) 1.0 15 LJ(18,9) LJ(18,9) LJ(24,12) LJ(24,12) LJ(36,18) LJ(36,18) n γ10 / γ0.5 γ /n = R 5 0.0 0 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 R R FIG.4. (a)AlldataofFig.3inacommonplotinwhichopentrianglesarecrystallinestatesandfulltrianglesareliquid,fluid, gas,orcoexistencestates. (b)Thesamedatawithγ/nplottedversusR. Theblackdashedlineconnectsthepoints(0,0)with (1,1), demonstrating the universal, roughly linear dependence of Eq. (26). Significant deviations fromthis are seen mainlyfor R>0.9andR<0. TheformercasereferstotheRsimpleregionofthephasediagram,thelattertostatesofcoexistingliquid and gas phases (compare to Fig. 2). 9 henceforth focus on the region in which R<0.9. Figure 4(b) plots the data differently and establishes for most of the R < 0.9 data the following approximate proportionality: γ/n(cid:39)R. (26) Howtounderstandthis? InRef. 2itwasshownthatanLJ-typepairpotentialmaybeapproximatedbythe“extended inverse power law (eIPL)” pair potential with an effective state-point-dependent exponent p that is not simply the exponent of the repulsive term of the LJ-type potential: v(r) ∼= Ar−p+B+Cr. (27) This approximation usually works very well within the entire first coordination shell2. At typical condensed-matter low and moderate pressure state points the term B+Cr contributes little to the forces on a given particle, because thesumofthesetermsoverallnearestneighborstendstobealmostconstant. Thisreflectsthefactthatiftheparticle is moved slightly, some nearest-neighbor distances increase and some decrease, but their sum stays almost constant2. This implies that the total force contribution from the B+Cr term is small and the physics is dominated by the r−p IPL term. For an IPL pair potential ∝r−p the density-scaling exponent is given in d dimensions by2 p γ = . (28) d As mentioned, the effective exponent p of Eq. (27) varies with state point. In large parts of the phase diagram of the LJ(m,n) pair potential the effective exponent p is considerably larger than m2. This is because even within the repulsive part of the pair potential, i.e., at distances below the potential minimum, there is a sizable contribution from the attractive r−n term making the repulsions significantly steeper than expected from the r−m repulsive term alone. At densities and temperatures dominated by the pair potential minimum, i.e., at typical low or moderate pressurecondensed-matterstatepoints, itmaybeshownfromacurvatureargumentreferringtothepotential-energy minimum that the effective exponent of the generalized LJ pair potential Eq. (20) is given2 by p = m+n. (29) Atsufficientlyhighdensityand/ortemperaturetherepulsiver−m termdominatesthephysics, however, andhereone finds that p→m. These conditions arise as R→1. For the standard LJ pair potential at typical state points the density-scaling exponent is between 5 and 6, corre- sponding via Eq. (28) to p between 15 and 181,2. Based on the above, we expect that the LJ(2n,n) pair potential in a large part of its phase diagram is equivalent to an eIPL pair potential with p (cid:39) 2n+n = 3n. Via Eq. (28), this corresponds in three dimensions to γ (cid:39)n. (30) Insummary,R→1atstatepointsofhighdensityand/ortemperaturewherether−2n termdominatesthephysics. This behavior is hinted at in Fig. 4(b) above R(cid:39)0.9 where γ/n starts to decrease and systematically deviates from the black dashed line. The limit γ (cid:39) 2n/3 was not reached in our simulations, however, which shows that only at very high density or temperature the attractive term of the generalized LJ pair potential may be ignored. A linear extrapolation of most of the data to R=1 gives the limit value suggested by the eIPL approximation in conjunction with Eq. (29). We take this as an indication that the eIPL approximation works well at the state points simulated, whether or not these are characterized by strong virial potential-energy correlations. Havingidentifiedthelimitsγ →0forR→0andγ →nforR→1(aslongasR<0.9,i.e.,formostdata),thenext question one may ask is: why do the data follow the approximate linear relation Eq. (26)? A possible explanation refersagaintotheeIPLapproximation. Thevirialiscalculatedinddimensionsfromthepairpotentialasasumover terms of the form rv(cid:48)(r)/d1. For the eIPL pair potential Eq. (27) this gives Apr−p/d plus an approximate constant (compare the above argument). Thus for fluctuations away from the equilibrium value, if the linear term of Eq. (27) is uncorrelated to the IPL term, the eIPL approximation suggests that one may have effectively 10 ∆W (cid:39)(p/d)∆U (31) with p ∼= 3n for d = 3. This must be, admittedly, a very rough approximation when the W U correlations are not strong, but if one nevertheless makes it, Eq. (1) and Eq. (16) imply (cid:112) (cid:112) γ/R(cid:39) (cid:104)(∆W)2(cid:105)/ (cid:104)(∆U)2(cid:105)(cid:39)p/d(cid:39)n. (32) This is Eq. (26) that summarizes our findings for most state points in three dimensions. A concise way of expressing thisisthattheexponentγ ofEq.(18)isroughlyconstantthroughoutthenon-R-simpleregionofthephasediagram, 2 compare Eq. (19). B. The Lennard-Jones pair potential in two, three, and four dimensions We proceed to report results for the standard LJ system in two, three, and four dimensions. Details on how the simulations were performed can be found in Refs. 25 and 28. Figure 5(a) shows the density-scaling exponent versus the correlation coefficient along three configurational adiabats in the liquid phase and two in the crystalline phase for the two-dimensional LJ system; (b) shows a similar plot along two configurational adiabats in the liquid phase of the three-dimensional LJ system. 9.5 (a) 2d (b) 3d 9.0 5.6 8.5 5.4 8.0 5.2 γ7.5 γ 5.0 7.0 T0=0.515 ρ0=0.600 T=1.312 ρ=0.632 T=0.515 ρ=0.710 0 0 0 0 4.8 T=1.312 ρ=0.948 6.5 T0=0.515 ρ0=0.750 0 0 T=0.515 ρ=0.900 0 0 6.0 T0=0.515 ρ0=1.065 4.6 0.6 0.7 0.8 0.9 1.0 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.00 R R FIG.5. Thedensity-scalingexponent,γ,versusthecorrelationcoefficient,R,alongconfigurationaladiabatsintwoandthree dimensionswithlegendsgivingthereferencestatepoints. (a)Resultsforthreeconfigurationaladiabatsintheliquidphaseand two in the crystalline phase generated using Eq. (17) for the 2d LJ system. Triangles represent liquid states, crosses represent solid states. The density-scaling exponent increases with increasing R for values of R below the 0.9 threshold and then starts to decrease. For the crystalline phase, R is always larger than 0.9 and γ decreases with increasing R (towards 12/2=6). (b) Two configurational adiabats for the liquid phase of the 3d LJ system, constructed as in the 2d case. The decrease in γ starts a bit after the 0.9 threshold. Figure 6 shows the density-scaling exponent versus the correlation coefficient along the configurational adiabats defined by the critical points for the LJ system in two, three, and four dimensions. Data for the critical points are taken from Refs. 29–31. For the configurational adiabats in the liquid phase the relation between R crossing the thresholdvalue0.9andthestartofthedecreasingbehaviorofγ isobserved. Forthecrystallinephase,thecorrelation coefficient is above 0.9 for all the state points studied and no crossover is observed. In d dimensions the eIPL-justified relation Eq. (32) implies via p=2n+n=3n that γd ∼= 3nR. (33) This is tested in Fig. 7(a) for the LJ simulations in 2−4 dimensions. We see that Eq. (33) overall works well. Finally, Fig. 7(b) tests all data presented in this paper versus Eq. (33). The data in the region where the system is not R simple conforms roughly to Eq. (33), unless R<0 which are state points of liquid-gas coexistence.