An extreme case of density scaling: The Weeks-Chandler-Andersen system at low temperatures
AAn extreme case of density scaling: The Weeks-Chandler-Andersen system at lowtemperatures
Eman Attia, Jeppe C. Dyre, Ulf R. Pedersen “Glass and Time”, IMFUFA, Dept. of Science and Environment,Roskilde University, P. O. Box 260, DK-4000 Roskilde, Denmark (Dated: January 27, 2021)This paper studies numerically the Weeks-Chandler-Andersen (WCA) system, which is shown toobey hidden scale invariance with a density-scaling exponent that varies from below 5 to above 500.This unprecedented variation makes it advantageous to use the fourth-order Runge-Kutta algorithmfor tracing out isomorphs. Good isomorph invariance of the structure and dynamics is observed overmore than three orders of magnitude temperature variation. For all state points studied, the virialpotential-energy correlation coefficient and the density-scaling exponent are controlled mainly by thetemperature. A mean-field theory is presented based on the assumption of statistically independentpair interactions, which rationalizes this finding. a r X i v : . [ c ond - m a t . s o f t ] J a n I. INTRODUCTION
Density scaling (thermodynamic scaling) is an important experimental discovery of the last 20 years of liquid-stateresearch, which simplifies high-pressure data for hundreds of systems [1–4]. The crucial insight is that, in order tocharacterize a thermodynamic state point, the relevant variable supplementing the temperature T is not the pressure p , but the (number) density ρ ≡ N/V (considering N particles in volume V ) [1, 2, 5, 6]. According to density scalingas it is typically used in the analysis of experimental data, if γ is the density-scaling exponent, plotting data for thedynamics as a function of ρ γ /T results in a collapse [1–4]. In other words, the dynamics depends on the two variablesof the thermodynamic phase diagram only via the single number ρ γ /T . This provides a significant rationalization ofdata, as well as an important hint for theory development. It should be noted, though, that density scaling is notuniversally applicable; for instance, it works better for van der Waals liquids than for hydrogen-bonded liquids [2, 4].Some time after these developments were initiated, a framework for density scaling was presented in terms ofthe isomorph theory [7, 8], which also provides a link to Rosenfeld’s 1977 excess-entropy-scaling principle [9, 10].According to isomorph theory, any system with strong correlations between the constant-volume virial and potential-energy equilibrium fluctuations has curves of invariant structure and dynamics in the thermodynamic phase diagram.These so-called isomorphs [7] are defined as the curves of constant excess entropy, i.e., as the system’s configurationaladiabats [11]. If the potential energy is denoted by U and the virial by W , the U W
Pearson correlation coefficient R is given by R = (cid:104) ∆ U ∆ W (cid:105) (cid:112) (cid:104) (∆ U ) (cid:105)(cid:104) (∆ W ) (cid:105) . (1)Here ∆ denotes the deviation from the thermal average and the sharp brackets denote canonical ( N V T ) averages.The pragmatic criterion used for defining “strong” correlations is
R > . R < . ρ γ /T =Const. witha constant γ [31]. This immediately explains density scaling. Isomorph theory, however, does not require that γ is constant throughout the phase diagram, and γ does indeed vary in most simulations [16, 32–34]. The generalisomorph-theory definition of the density-scaling exponent γ at a given state point in the thermodynamic phasediagram is [7, 10] γ ≡ (cid:18) ∂ ln T∂ ln ρ (cid:19) S ex = (cid:104) ∆ U ∆ W (cid:105)(cid:104) (∆ U ) (cid:105) . (2)We here quote also the statistical-mechanical formula used for calculating γ from canonical ensemble ( N V T ) fluctu-ations of the potential energy and the virial.The question whether the density-scaling exponent in experiments is strictly constant or may vary throughout thephase diagram has recently come into focus [35, 36]. Theoretically, isomorphs are in many cases described by thefollowing equation [1] h ( ρ ) T = Const. (3)in which h ( ρ ) is a function of density that can be calculated in various way. For the Lennard-Jones (LJ) system, onehas h ( ρ ) ∝ ( γ / − ρ/ρ ) − ( γ / − ρ/ρ ) in which γ is the density-scaling exponent at a selected “referencestate point” of density ρ [32, 37]. For isomorphs given by Eq. (3), Eq. (2) implies γ = d ln h ( ρ ) d ln ρ . (4)Clearly, unless h ( ρ ) is a power-law function, the density-scaling exponent changes if density if changed. This appliesfor many of the systems that have been studied by computer simulation [8, 32, 37]. More generally, however, γ alsodepends on the temperature [33]. This applies, for instance, for the LJ system at very high temperatures: for T → ∞ at a fixed density this system becomes dominated by the repulsive r − term of the pair potential, implying that γ must approach 12 / γ is the fact that density often does notvary very much. When extreme pressures are applied, the density-scaling exponent stops being constant as Casaliniand Ransom have recently documented for several systems [36]. These authors noted, moreover, that γ converges to4 at high pressures, an interesting observation which is not discussed further here.Although there is agreement that the density-scaling exponent generally is not constant, the variation of γ is oftenquite modest and as mentioned insignificant in many experiments. The present paper gives an example in which γ varies dramatically. We present a study of the noted Weeks-Chandler-Andersen (WCA) system [38], which is shownto have strong virial potential-energy correlations and thus be R-simple at typical liquid-state densities. We find that γ varies by more than two decades in the investigated part of the phase diagram. To the best of our knowledge, this isa much larger variation than has previously reported for any system in simulations or experiments. For all the statepoints studied below, which include state points on isomorphs, isotherms, and isochores, we find that γ primarilydepends on the temperature. A mean-field theory is presented that explains this observation and which accountsquantitatively for the low-temperature behavior of the system.After providing a few technical details in Sec. II, the paper starts by presenting the state points studied numericallyin the thermodynamic phase diagram (Sec. III). The paper’s focus is on three isomorphs, numbered 1-3. Eachisomorph is associated with one isotherm and one isochore (constant volume). The purpose is to put into perspectivethe variation of structure and dynamics along the isomorphs by comparing to what happens along isotherms andisochores with similar density/temperature variations. In Sec. III we also give data for the virial potential-energycorrelation coefficient R and the density-scaling exponent γ , demonstrating that all state points studied have strongcorrelations ( R > .
9) and that γ varies from about 5 to above 500. An analytical mean-field theory is developedin Sec. IV predicting that R and γ both depend primarily on the temperature. Section V presents simulations ofthe structure and dynamics along isotherms, isochores, and isomorphs. Despite the extreme γ variation, meaningthat an approximate inverse-power-law description fails completely, we find pretty good isomorph invariance of thereduced-unit structure and, in fact, excellent isomorph invariance of the reduced-unit dynamics. Section Sec. VI givesa few concluding remarks. The Appendix details the implementation of the fourth-order Runge-Kutta method fortracing out isomorphs and compares its predictions to those of the previously used simple Euler method. II. MODEL AND SIMULATION DETAILS
Liquid model systems are often defined in terms of a pair potential v ( r ). If r ij = | r i − r j | is the distance betweenparticles i and j , the potential energy U as a function of all particle coordinates R ≡ ( r , r , .., r N ) is given by U ( R ) = (cid:88) i Figure 1(a) shows the thermodynamic phase diagram of the WCA system. The solid yellow/orange lines are thefreezing and melting lines [42, 43]. Three solid lines mark the isomorphs in focus below (1-3), while the red dashedline is a fourth isomorph (0), which is generated from a state point on the melting line running slightly into theliquid-solid coexistence region. Each isomorph was traced out starting from a reference state point of density 0.84.Isomorphs are often identified by integrating Eq. (2) using a simple first-order Euler scheme for a density changeof order one percent [7, 16, 20]. The extreme variation of γ for the WCA system reported below, however, meansthat Euler integration can only be used reliably for much smaller density changes. Consequently, a more accurateintegration scheme is called for. We used standard fourth-order Runge-Kutta integration (RK4) as detailed in theAppendix, in which it is demonstrated that RK4 is significantly more efficient that Euler integration for tracing outisomorphs with a given accuracy.In order to investigate the degree of isomorph invariance of the reduced-unit structure and dynamics (Sec. IV), foreach isomorph we also performed simulations along an isotherm and an isochore through the state point marked bya black star (limiting the simulations to state points in the equilibrium liquid phase). The state points simulated onthe isotherms and isochores are marked as small purple crosses in Fig. 1(a). Figure 1(b) shows the isomorphs and themelting line in a diagram with a logarithmic temperature axis. Note that the melting line is an approximate isomorph[44].A configurational adiabat is only an isomorph for state points with strong virial potential-energy correlations, i.e.,if R > . R is given by Eq. (1). That this is the case is demonstrated in Fig. 2,which shows R for all state points simulated. (a) shows R as a function of density, while (b) shows R as a functionof temperature. We see that R increases with density and temperature, approaching unity as ρ → ∞ or T → ∞ . Density, [ ] C o rr e l a t i o n C o e ff i c i e n t , R (a) R = 0.921 Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 Temperature, T [ / k B ] C o rr e l a t i o n C o e ff i c i e n t , R (b) R = 0.921 Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T=12.1 FIG. 2. The virial potential-energy Pearson correlation coefficient R for all simulated state points on the isomorphs, isotherms,and isochores (compare Fig. 1). There are strong correlations ( R > . 9) everywhere. The horizontal dashed lines mark thelow-temperature, low-density limit of the mean-field prediction (Eq. (23)). (a) R as a function of the density. (b) R as afunction of the temperature. Clearly, the correlations are mainly controlled by the temperature. This reflects the fact that the ( r/σ ) − term of the pair potential dominates the interactions in these limits and thatany inverse power-law pair potential has R = 1. An interesting observation of Fig. 2 is that strong correlations aremaintained even at low densities and temperatures. Comparing (a) to (b) reveals that R is primarily controlled by thetemperature. This may be understood from a simple analytical mean-field theory, which assumes that the interactionsat low temperatures and low densities are dominated by single-pair interactions (Sec. IV). Density, [ ] g a mm a , ( a ) Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 Pressure, p [ / ] g a mm a , ( b ) Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 Temperature, T [ / k B ] g a mm a , ( c ) Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 Density, [ ] g a mm a , ( d ) Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 Pressure, p [ / ] g a mm a , ( e ) Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 Temperature, T [ / k B ] g a mm a , ( f ) T Isomorph 0Isomorph 1Isomorph 2Isomorph 3Isochore: =0.84Isochore: =1.0Isochore: =1.2Isotherm: T =0.6Isotherm: T =2.7Isotherm: T =12.1 FIG. 3. The density-scaling exponent γ at all state points on the isomorphs, isotherms, and isochores studied (compare Fig. 1).(a) γ as a function of the density. (b) γ as a function of the pressure. (c) γ as a function of the temperature. (d) γ as afunction of the density in a log-log plot, showing data for all the state points simulated. (e) γ as a function of the pressure ina log-log plot, showing data for all the state points simulated. (f) γ as a function of the temperature in a log-log plot, showingdata for the all state points simulated. We see that γ is primarily a function of the temperature. The dashed line in (f) marksthe mean-field prediction in the limit T → 0, which here fits data well (Eq. (22)). Figure 3 gives data for the density-scaling exponent γ at all the state points simulated plotted in different ways usingthe same color coding as in Fig. 2. We note that γ increases monotonically as the density, pressure, or temperature islowered, eventually reaching values above 500. Figure 3(a) shows γ as a function of the density ρ . Clearly, knowledgeof ρ is not enough to determine γ . This means that Eq. (3) does not apply for the WCA system. It has been suggestedthat γ is determined exclusively by the pressure [45]. This works better than the density for collapsing data, thoughstill with considerable scatter ((b)). Figure 3(c) plots γ as a function of the temperature. We here observe a quitegood collapse (see also (f) with a logarithmic temperature axis). IV. MEAN-FIELD THEORY FOR R AND γ AT LOW DENSITIES This section presents a pair-based theory for estimating the virial potential-energy correlation coefficient R and thedensity-scaling exponent γ . The theory, which is inspired by Refs. 18, 33, and 46, assumes that the individual pairinteractions are statistically independent. This is only a good approximation at low densities, but as becomes clearfrom the below treatment what “low” means in this context depends a lot on the temperature.In MD units the truncated WCA potential of Eq. (6) is v ( r ) = 4 r − − r − + 1 for r < r c = 2 / = 1 . . . . (7)and zero otherwise. This potential is our focus henceforth, but it should be emphasized that the theory developedapplies for any purely repulsive truncated potential.The virial of the configuration R is given by W ( R ) = (cid:80) Ni>j w ( r ij ) in which the pair virial is defined as w ( r ) ≡− ( r/ v (cid:48) ( r ) [40]. If pair interactions are statistically independent, in order to calculate R and γ it is enough toconsider the potential-energy and virial fluctuations of a single pair. The remaining particles are treated in a mean-field approximation (Fig. 4). Let one particle be fixed at the origin. The particles it interacts with are found in asphere of radius r s , while the states with r > r s are excluded by the mean-field of the other particles. We expect thisapproximation to be good at low densities where at any given point in time a particle at most interacts with a singleother particle.The partition function of the canonical ensemble is given by Z ( ρ, T ) = (cid:90) r s p ( r ) dr (8)in which p ( r ) = 4 πr exp( − v ( r ) /k B T ) (9)is the un-normalized probability distribution. Since v ( r ) = 0 for r > r c , we can split the integral of the partitionfunction into two parts, writing Z ( ρ, T ) = Z ( T ) + V free ( ρ ) (10)in which Z ( T ) = (cid:90) r c p ( r ) dr . (11)Here Z ( T ) is independent of the density and V free ( ρ ) = (cid:82) r s r c πr dr is the volume in which the particles do not interact(the white region in Fig. 4). The partial volume available to each particle is 1 /ρ . The volume excluded by any givenparticle is the volume of a sphere of diameter r c , i.e. πr c / 6. The difference between these two volume gives the partof phase space where the particles do not interact, V free ( ρ ) = 1 /ρ − πr c / . (12)At high densities near random close packing (rcp) or in a dense face centered cubic (fcc) crystal V free ( ρ ) will besmaller due to many-body steric interactions. A high-density approximation for V free ( ρ ) could be made by includingknowledge of the fcc or rcp packing fraction, but here we are interested in the low density limit.Expectation values are now computed as (cid:104) A (cid:105) = (cid:90) r c A ( r ) p ( r ) dr/Z ( ρ, T ) (13) r r c r s V free FIG. 4. Sketch of a pair of particles within the mean-field approximation. The light blue/grey area represents the volume inwhich the particles interact. The white area represents the free volume V free . The brown area represents the volume excludedby other particles treated as a mean field. The dashed lines represent the volumes excluded by each particle, i.e., spheres ofdiameter r c . where A ( r ) is either v ( r ), w ( r ), [ v ( r )] , [ w ( r )] or w ( r ) v ( r ). In this way one calculates γ ( ρ, T ) = (cid:104) wv (cid:105) − (cid:104) w (cid:105)(cid:104) v (cid:105)(cid:104) v (cid:105) − (cid:104) v (cid:105) (14)and R ( ρ, T ) = (cid:104) wv (cid:105) − (cid:104) w (cid:105)(cid:104) v (cid:105) (cid:112) ( (cid:104) w (cid:105) − (cid:104) w (cid:105) )( (cid:104) v (cid:105) − (cid:104) v (cid:105) ) . (15)Figure 5 compares the predictions of the mean-field theory to data for isomorphs and isochores. There is goodoverall agreement. Systematic deviations are visible in (b) and (d), however, that focus on densities higher than thosefor which the mean-field theory may reasonably be expected to apply.We proceed to discuss the low-density limit of the mean-field theory. In the limit ρ → /Z → 0. Theterms that involve a single expectation value ( (cid:104) v (cid:105) , (cid:104) w (cid:105) , and (cid:104) wu (cid:105) ) scale as 1 /Z while the terms that involve amultiplication of expectation values ( (cid:104) v (cid:105) , (cid:104) w (cid:105) , and (cid:104) v (cid:105)(cid:104) w (cid:105) ) scale as [1 /Z ] . At low densities Z is large and one canneglect terms that involve multiplications of expectation values, leading [18, 33] to γ ( T ) = (cid:104) wv (cid:105) / (cid:104) v (cid:105) for ρ → R ( T ) = (cid:104) wv (cid:105) / (cid:112) (cid:104) w (cid:105)(cid:104) v (cid:105) for ρ → . (17)Notice that these averages do not depend on Z since both the numerators and denominators scale as 1 /Z . This iswhy γ and R at low densities both become independent of ρ , i.e., functions only of T .At low temperatures the probability distribution p ( r ) concentrate near r c since v ( r ) is monotonic and decreaseswith increasing r . Thus, one can make an expansion around x = r c − r = 0, writing the pair potential as v ( x ) = k x + k x / k x / ..... (18)The pair virial then becomes [15] w ( x ) = k ( r c − x ) + ( r c − x ) k x/ k r c x / O ( x ) . (19) Temperature, T [ / k B ] C o rr e l a t i o n C o e ff i c i e n t , R (a) Isomorph 1Isomorph 2Isomorph 3Isomorph 1 (Theory)Isomorph 2 (Theory)Isomorph 3 (Theory) Temperature, T [ / k B ] C o rr e l a t i o n C o e ff i c i e n t , R (b) Isochore: =0.84Isochore: =1.0Isochore: =1.2Isochore 1 (Theory)Isochore 2 (Theory)Isochore 3 (Theory) Temperature, T [ / k B ] g a mm a , (c) Isomorph 1Isomorph 2Isomorph 3Isomorph 1 (Theory)Isomorph 2 (Theory)Isomorph 3 (Theory) Temperature, T [ / k B ] g a mm a , (d) Isochore: =0.84Isochore: =1.0Isochore: =1.2Isochore 1 (Theory)Isochore 2 (Theory)Isochore 3 (Theory) FIG. 5. Comparing predictions of the mean-field theory for R and γ to the simulation results. (a) and (c) show results as afunction of temperature along each of the four isomorphs. (b) and (d) show results along the three isochores, allowing for afocus on the higher-temperature regime where the mean-field theory is expected to break down gradually. We conclude thatthe theory works well at low temperatures. For the WCA potential one has k = 0 and k = 36 √ 4, and expectation values can now be written as (cid:104) A (cid:105) = (cid:90) ∞ A ( x ) p ( x ) dx/Z ( T → 0) (20)in which p ( x ) = 4 π ( r c − x ) exp( − k x / (2 k B T )) ∼ = 4 πr c exp( − k x / (2 k B T )) . (21)Since p ( x ) is concentrated near x = 0, the upper limit of the integral Eq. (20) is extended to infinity. This implies thatthe integral of Z ( T ) and the expectations values of interest are Gaussian integrals that can be evaluated analytically.In the low-temperature limit, γ and R become γ = 4 r c √ k √ πT = 163 √ πT (22)and R = (cid:114) π = 0 . . . . (23)Figure 6 shows the mean-field prediction for γ and R at T = 0 . 01 plotted as a function of the density. We see thatthe theory works well, even though one is still far from the T → V. VARIATION OF STRUCTURE AND DYNAMICS ALONG ISOTHERMS, ISOCHORES, ANDISOMORPHS This section investigates to which degree the reduced-unit structure and dynamics are invariant along three iso-morphs. This is not obvious in view of the considerable γ variation, which means that the situation is very far from Density, [ ] g a mm a , (a)T =0.01 = 16/3 T from mean field theory from simulations from simulations- Crystalized Density, [ ] C o rr e l a t i o n C o e ff i c i e n t , R (b)T =0.01 R = 8/3 R from mean field theory R from simulations R from simulations- Crystalized FIG. 6. Low-temperature ( T = 0 . 01) density dependence of (a) γ and (b) R . Isotherm: T = 0.6Isotherm: T = 0.6Isotherm: T = 0.6 Isotherm: T = 2.7Isotherm: T = 2.7Isotherm: T = 2.7 Isotherm: T = 12.1Isotherm: T = 12.1Isotherm: T = 12.1 R a d i a l d i s t r i b u t i o n f un c t i o n , g ( r ) Isochore: = 0.84Isochore: = 0.84Isochore: = 0.84 Isochore: = 1.0Isochore: = 1.0Isochore: = 1.0 Isochore: = 1.2Isochore: = 1.2Isochore: = 1.2 Isomorph 1Isomorph 1Isomorph 1 Pair distance, r [ ] Isomorph 2Isomorph 2Isomorph 2 Isomorph 3Isomorph 3Isomorph 3 FIG. 7. Radial distribution functions (RDF) for the three isotherms, isochores, and isomorphs. Although the height of the firstpeak is not isomorph invariant, in comparison to isotherms and isochores we see a significantly better invariance for the RDFalong the isomorphs. This is the case even though the density variation for the isotherm and the temperature variation for theisochore are both smaller than those of the isomorph (compare Fig. 1). that described by an approximate Euler-homogeneous potential-energy function. Isomorph invariance is rarely exact,so in order to have something to compare to, we present also results for the variation of the reduced-unit structureand dynamics along isotherms and isochores (Fig. 1). As a measure of the structure, we look at the reduced radialdistribution function (RDF) g (˜ r ). As a measure of the dynamics, we look at the reduced mean-square displacement(MSD) as a function of reduced time, as well as on the reduced diffusion coefficient identified from the long-timebehavior of the MSD.Starting with structure, Fig. 7 shows reduced-unit RDF data along the three isotherms, isochores, and isomorphs01-3 (Fig. 1). The isotherms span roughly the same density range and the isochores roughly the same temperaturerange as the corresponding isomorphs, subject to the requirement that we only give data for the equilibrium liquidphase, i.e., above the freezing line. The first observation is that the RDF shows some variation along the isomorph(lowest row of Fig. 7). In comparison to the isotherms and isochores, however, there is clearly approximate isomorphinvariance of the RDF, in particular beyond the first minimum. The first peak height is not isomorph invariant; forall three isomorphs it is significantly higher at low temperatures than at higher temperatures. This is an effect ofthe steepness of the effective pair potential, with large γ resulting in a higher first peak. The argument is as follows.Consider the IPL pair potential ∝ r − n system. It has γ = n/ n is, the moreharshly repulsive are the forces. From the Boltzmann probability factor of finding two particles at distance r , whichis ∝ exp( − v ( r ) /k B T ), it is clear that particle near encounters become less likely as n → ∞ , thus suppressing theRDF at small distances. If there is isomorph invariance of the number of particles within the first coordination shell,as n (or γ = n/ 3) increases some of the RDF must therefore move from small r to somewhat larger r . This results ina higher first peak. Very recently, this argument has been strengthened and formalized by the observation that theso-called bridge function is isomorph invariant to a very good approximation [48].A similar increase with increasing γ of the height of the first RDF peak has been reported for the EXP system (Fig.5 in Ref. 33). In that case it was a much less dramatic effect, though, because the EXP system’s γ variation at theinvestigated state points is less than a factor of 3 compared to more than a factor of 100 for the present WCA statepoints. Interestingly, for both systems the data indicate that γ → ∞ as T → Isotherm: T = 0.6Isotherm: T = 0.6Isotherm: T = 0.6 Isotherm: T = 2.7Isotherm: T = 2.7Isotherm: T = 2.7 Isotherm: T = 12.1Isotherm: T = 12.1Isotherm: T = 12.1 M e a n s q u a r e d d i s p l a c e m e n t Isochore: = 0.84Isochore: = 0.84Isochore: = 0.84 Isochore: = 1.0Isochore: = 1.0Isochore: = 1.0 Isochore: = 1.2Isochore: = 1.2Isochore: = 1.2 Isomorph 1Isomorph 1Isomorph 1 Time, t [ m / T ] Isomorph 2Isomorph 2Isomorph 2 Isomorph 3Isomorph 3Isomorph 3 FIG. 8. The mean-square displacement (MSD) plotted against time in reduced units for the three isotherms, isochores, andisomorphs. The dynamics is isomorph invariant to a good approximation. From end point to end point of the isomorphs,compare Fig. 1, the variation in ˜ D (determined from the MSD long-time variation) is, respectively, 39%, 23%, and 14%. Thecorresponding numbers are 1000%, 880%, and 549% along the isochores, and 214%, 893%, and 305% along the isotherms. Proceeding to investigate the dynamics, Fig. 8 shows data for the reduced-unit MSD as a function of the reducedtime along the three isotherms, isochores, and isomorphs. There is only good invariance along the isomorphs. In MDunits, the MSDs are also not invariant along the isotherms or isochores (data not shown). This means that the lack1of invariance for the isotherms and isochores is not a consequence of the use of reduced units. With Fig. 7 in mind,we conclude that the varying first-peak heights of the RDFs along the isomorphs do not influence the dynamics verymuch. This is consistent with expectations from liquid-state quasiuniversality, according to which many systems havestructure and dynamics similar to that of the hard-sphere or EXP generic liquid systems [18, 49]. D [ / m ] Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1Isomorph 1 Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2Isomorph 2 Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3Isomorph 3 D = D / m / k B T Temperature, T [ / k B ] FIG. 9. Diffusion coefficients along the three isomorphs 1-3 in MD units (upper row) and in reduced units (lower row), plottedas functions of the logarithm of the temperature. The diffusion coefficients vary a lot along the isomorphs when given in MDunits, but are fairly constant in reduced units. This demonstrates the importance of using reduced units when checking forisomorph invariance. The reduced diffusion coefficient ˜ D is extracted from the data in Fig. 8 by making use of the fact that the long-timereduced MSD is 6 ˜ D ˜ t . We show in Fig. 9 how both D and ˜ D vary along the three isomorphs. The upper figuresdemonstrate a large variation in D along each isomorph. The lower figures show ˜ D , which is rigorously invariant ifthe system had perfect isomorphs. This is not the case, but the variation is below 40% for all three isomorphs insituations where temperature varies by roughly a factor of 10000. We conclude that the reduced diffusion coefficientis isomorph invariant to a good approximation.Though there is some numerical uncertainty, it appears from Fig. 9 that ˜ D stabilizes at the lowest temperature,and for each isomorph one can tentatively identify a low-temperature limit. Figure 10 plots the low-temperaturelimits of the three isomorphs versus data obtained at the lowest density simulated (supplemented by data for onemore isomorph). We see that higher low-temperature densities correlate with lower ˜ D . An obvious question is: whichdensity corresponds to ˜ D = 0? At very low temperatures the WCA system behaves more and more like a system ofhard spheres (HS) because γ becomes very large. The disordered HS system has a maximum density, the randomclosed-packed (rcp) density of roughly 63% packing fraction. In Fig. 10, the open black star at ˜ D = 0 marks thedensity corresponding to RCP. Within the uncertainty, the data are consistent with a convergence to this point. VI. SUMMARY This paper has shown that the WCA systems presents a striking example of the density-scaling exponent being farfrom approximately constant throughout the thermodynamic phase diagram [35, 36]. We have also shown that theWCA system is R-simple and has the expected isomorph invariances, although the first peak of the structure factordoes vary significantly over the roughly four decades of temperature variation covered.We studied three isomorphs of the WCA system and showed that the density-scaling exponent along them vary bymore than a factor of 100. This extreme variation means that, in contrast to, e.g., the LJ system, the WCA systemcannot be described even approximately as an effective IPL system [15]. In the LJ case, the pair potential may be2 D = D m / k B T D e n s i t y , [ ] Isomorph 0Isomorph 1Isomorph 2Isomorph 3 FIG. 10. Low-temperature limiting reduced diffusion coefficient for the three isomorphs 1-3, supplemented by data for onemore isomorph, plotted versus the density of the lowest-temperature state point simulated on the isomorph in question. Thepoints are fitted by a cubic spline function (curve), which by construction goes through the calculated random close-packingdensity ( ρ = 0 . D → ρ = 0 . 864 is calculated at as follows. With r c = 2 / = 1 . V sphere = πr c / . ρV sphere , one arrives at ρ = 0 . approximated by the so-called extended IPL pair potential, which is a sum of an IPL term ∼ r − , a constant, anda term proportional to r [15]. The latter two terms contribute little to the fluctuations of virial and potential energy[15], which explains the strong correlations for the LJ system and why one here observes a γ that is fairly close to6 (not to 4 as one might guess from the repulsive r − term of the potential). The WCA situation is very different.Because that system has no liquid-gas phase transition and thus no liquid-gas coexistence region, isomorphs may bestudied over several orders of magnitude of temperature and, in particular, followed to very low temperatures withoutloosing the strong-correlation property. In this process, γ increases in an unprecedented fashion. Nevertheless, thereduced-unit structure and dynamics are fairly invariant along the isomorphs. The result that R and γ are bothprimarily functions of the temperature can be explained by a simple mean-field theory based on the assumption ofstatistically independent pair interactions. ACKNOWLEDGMENTS This work was supported by the VILLUM Foundation’s Matter grant (16515). APPENDIX: USING THE RUNGE-KUTTA METHOD FOR EFFICIENTLY TRACING OUTCONFIGURATIONAL ADIABATS The density-scaling exponent γ is the slope of the lines of constant S ex in the (ln T, ln ρ ) plane (Eq. (2)). Thus bynumerical integration one can from Eq. (2) compute the lines of constant S ex , the configurational adiabats (which areisomorphs for any R-simple system). The density-scaling exponents required for the integration can be determinedfrom the thermal equilibrium virial potential-energy fluctuations in an N V T simulation (Eq. (2)). In the followingwe denote the theoretical slope by f , i.e., the slope without the unavoidable statistical noise of any MD simulation.Let ( x, y ) be (ln ρ, ln T ) (occasionally it is better to choose ( x, y ) = (ln T, ln ρ ) instead). In this notation, let dydx = f ( x, y ) (24)be the first-order differential equation one needs to integrate. Several methods have been developed to integrate thisnumerically [50]. The simplest one is Euler’s method: Imagine that one has estimated the slope at some point ( x i , y i )by computing γ = f ( x i , y i ) from the virial potential-energy fluctuations by means of Eq. (2). The point ( x i +1 , y i +1 )is then calculated from x i +1 = x i + hy i +1 = y i + hf ( x i , y i ) + O ( h ) . (25)3Here, h is the size of the numerical integration step along x . The “local” truncation error on the estimated y i +1 scalesas h . The statistical error on the numerical calculation of the slope f scales as 1 / √ τ where τ is the simulation time.Thus, the local statistical error on y i +1 scales as h/ √ τ (rounding errors from the finite machine precision are notrelevant for the h ’s we investigate here). The scaling of the total local error is thus proportional to h + ch/ √ τ inwhich c is a constant. However, we are interested in the global truncation error, i.e., the accumulated error for someintegration length ∆ x . Let N = ∆ x/h be the number of steps needed to complete the integration. If the wall-clockcomputation time is kept constant, the total simulation time available is t = N ( τ + τ eq ) where τ eq is the time it takesfor the system to come into equilibrium when temperature and density are changed. Thus τ = t/N − τ eq , and with h = ∆ x/N the local statistical error on y is ch/ √ τ = c ∆ x/ (cid:112) N t − N τ eq . The global error from truncation scales as N since it is systematic, while the statistical error scales as √ N due to its randomness. Thus, the total global erroris proportional to (∆ x ) /N + c ∆ x/ (cid:112) t − N τ eq . The first term, related to truncation error, is lowered by making N large, while the second term favors small N ’s and diverges as N → t/τ eq . Thus, since c is in general unknown, theoptimal choice of N for a given t and ∆ x is non-trivial to determine. In the following, we give a recipe for the optimalparameter choice. First, however, we reduce the truncation error by adopting a higher-order integration method,using the often favored fourth-order Runge-Kutta (RK4) method: For a given point ( x i , y i ), if one defines k = hf ( x i , y i ) k = hf ( x i + h/ , y i + k / k = hf ( x i + h/ , y i + k / 2) (26) k = hf ( x i + h, y i + k ) , the next point ( x i +1 , y i +1 ) is computed as x i +1 = x i + hy i +1 = y i + k / k / k / k / O ( h ) . (27)While the simple Euler method has a local truncation error scaling as O ( h ), the truncation error of RK4 scales as O ( h ). This allows for significantly larger steps along x and thus smaller N . From the same type of arguments as givenabove for the Euler method, the global error of the RK4 method scales approximately as (∆ x ) /N + c ∆ x/ (cid:112) t − N τ eq (here c is new unknown constant).To illustrate the methods’ accuracy, we use them for integrating from the initial state point ( ρ, T ) = (0 . , . . γ variation from 6.825at the intial density to 4.539 at ρ = 1 . 25. The difference between the final temperature of the down integration andthe initial temperature, denoted by ∆ T , provides a measure of the maximum temperature error. Ideally ∆ T = 0, butin practice ∆ T deviates from zero. Since the RK4 involves four simulations per step, we compare its accuracy where h is four times larger than for the Euler method, giving approximately the same wall-clock time for the computation.With this constraint, the RK4 is still about two orders of magnitude more accurate: we find ∆ T = 0 . 186 for the Euleralgorithm and ∆ T ∼ = 0 . 002 for RK4. Figure 12 shows estimates of the maximum error ∆ T for several values of h . Tofocus on the truncation error, we have here performed long-time simulations with τ ∼ = 650. Nonetheless, this analysisdemonstrates that a significantly smaller N (larger h ) is allowed for the RK4.Since the RK4 algorithm allows for large steps in h , it can be advantageous to perform interpolations to identifyadditional state-points on the configurational adiabat. The solid lines in Fig. 11 show such interpolations using acubic Hermite spline: Define x φ as a point between the two adjacent points x i and x i +1 , i.e. let x i ≤ x φ < x i +1 where x φ = x i + φ [ x i +1 − x i ] and 0 ≤ φ ≤ 1. The interpolated y φ is given by the smooth third-degree polynomial: y φ = Ax φ + Bx φ + Cx φ + D where y φ = y i + [ y i +1 − y i ][ aφ + bφ + cφ ]. For simplicity we introduce the notation y (cid:48) φ = [ y φ − y i ] / [ y i +1 − y i ] and write the polynomial as y (cid:48) φ = aφ + bφ + cφ . The coefficients yielding smooth firstderivative are a = f (cid:48) i + f (cid:48) i +1 − b = 3 − f (cid:48) i − f (cid:48) i +1 and c = f (cid:48) i where f (cid:48) i = f i x i +1 − x i y i +1 − y i and f (cid:48) i +1 = f i +1 x i +1 − x i y i +1 − y i are “reduced” slopes at the start- and end point, respectively. The f (cid:48) slopes are given by known γ ’s along theconfigurational adiabat, thus no extra simulations are needed to evaluate the interpolation.We investigated the local error by comparing a full h step to two half steps of size h/ 2. The black dots in Fig. 11(b)shows the result of such two half-steps. The truncation error for the half-step approach is then raised to the sixth order[50], one order higher than RK4 (the price is that we have to perform twice as many simulations for each integrationstep). The triangles on Fig. 13 show the resulting T i +1 starting from the reference state-point ( ρ, T ) = (0 . , . h = 0 . τ ’s. For comparison, the dashed line results from long-time simulationsusing the half-step algorithm. In effect, the distance from triangles to the dashed line estimates the total local error.4 . Density, [ ] T e m p e r a t u r e , T [ / k B ] (a) ForwardBackwardStart point 0.9 1.0 1.1 1.2 Density, [ ] T e m p e r a t u r e , T [ / k B ] (b) ForwardBackwardHalf step, h Start point FIG. 11. Configurational adiabat of the WCA system traced out in the thermodynamic phase diagram as functions of the steplength h = ∆ x/N in which N is the number of steps involved in the integration. (a) The Euler method; (b) the RK4 method.The Euler integration uses a log-density step of size h = 0 . e . − (cid:39) h = 0 . 4, corresponding to density variation of e . − (cid:39) T is a measure of the maximum error of the predicted temperature. We find ∆ T ∼ = 0 . 186 forthe Euler algorithm and ∆ T ∼ = 0 . 002 for the RK4 algorithm. The solid lines are interpolations using a cubic Hermite spline.. Step Size, h T e m p e r a t u r e E rr o r , T [ / k B ] (a) 0.84 1.30 T hT h EulerRK4 10 Step Size, h T e m p e r a t u r e E rr o r , T [ / k B ] (b) 0.58 0.84 T h T h EulerRK4 FIG. 12. (a) The temperature difference ∆ T of the forward-backward integration in Fig. 11, for different steps sizes h . Theblue dots show results for Euler and the orange dots for RK4 integration. The temperature difference is a measure of themaximum error of the integration interval 0 . ≤ ρ ≤ . 30. The RK4 is significantly more accurate than the Euler algorithm.This allows for larger h steps. The dashed lines indicate the expected scaling of the global error from truncation – deviationsstem from statistical errors on the estimated slopes (slopes are evaluated using simulations lengths of τ = 655). The arrowconnects two calculations with approximately the same computational cost (see Fig.11). (b) Same analysis for the integrationinterval 0 . ≤ ρ ≤ . For short simulation times (small τ ’s) the statistical error dominates as seen by the scatter. The truncation errordominates when using long simulation times, as seen by the triangles’ systematic deviation from the dashed line. Forefficient calculation we suggest choosing h and τ such that the statistical and truncation errors are of the same orderof magnitude. The red × on Fig. 13 indicate the simulation time τ used for figures in the main paper.The above analysis to arrive at the optimal computation time τ is tedious and involves computationally expensivelong-time simulations. We proceed to suggest an efficient optimization recipe that utilizes the fact that the localstatistical error of the slopes can be estimated by dividing a given simulation into blocks. If the simulation time foreach block is sufficiently long, the blocks can be considered as statistically independent. The 67% confidence standarderror is then given by SE( γ ) = (cid:112) VAR( γ ) /N B where VAR( γ ) is the variance of the γ ’s using N B blocks [51]. If theblocks are independent, VAR( γ ) scales as N B and SE( γ ) will be independent of the number of blocks. If we divide thesimulation into a few blocks, VAR( γ ) may give a bad estimate of the underlying distribution’s theoretical variance.On the other hand, if one divides the simulation into many blocks, the simulation time for each block ( τ /N B ) may bebrief and the blocks are not independent. In effect, the above formula for SE( γ ) gives an overestimate. The optimal N B is determined by tests of several different N B , as shown on Fig. 14 (the red × corresponds to a good choice of N B = 128). The statistical error on y i +1 can now be estimated as SE( y i +1 ) = SE( γ ) h/ 2. Here, 2 = √ √ Simulation Time, T i + , [ / k B ] i +1 = 1.25 Temperature with h = 0.4Half Step Approach FIG. 13. The difference in temperature between using a full step of h = 0 . h = 0 . ρ = 0 . 84 upto ρ = 1 . 25, plotted against simulation time per slope evaluation. The desired h can change and the simulation time changesaccordingly. The error bar at τ = 51 . × indicates thesimulation time used in the paper’s simulations. -2 -1 Block time, τ/N B Number of blocks, N B E s t i m a t e d S E ( γ ) τ = 51.2γ = 6.82 ± 0.05 FIG. 14. Estimate of the statistical error on γ from the blocking method. The analysis indicates that N B = 128 is an goodchoice for the number of blocks. This gives SE ( γ ) = 0 . 03 on the estimated γ = 6 . tional adiabat:1. Make an N V T simulation at a reference state point of temperature T and density ρ . The simulation time τ should be sufficiently long that the equilibration time τ eq can be determined using any standard method (e.g.,as the time when the mean-squared displacement has reached the diffusive limit). Use the block method todetermine SE( γ ), using only the equilibrated part of the trajectory.2. Choose an h step. Make a full RK4 step and estimate the local statistical error using SE( y i +1 ) = h SE( γ ) / √ τ if the statistical error is of the same magnitude as the total error, or(b) decrease h if the total error is larger than the statistical error.Small errors suggest that the simulation time, τ , could be decreased, or that h can be increased to make thecalculation more efficient. Moreover, h may safely be increased or τ decreased if the statistical and total errorsare of similar magnitude.3. Compute adiabatic state-points using the RK4 algorithm using the parameters determined in the above steps.Based on these, a continuous curve can be produced by interpolation of using a cubic spline.4. Estimate the maximum error by integrating backwards. This error estimate quantifies the accuracy of thecomputed adiabat.As a consistency check of this recipe, Fig. 15 shows the excess entropy from the equation of state (EOS) of thesingle-component LJ system by in Ref. 52. The agreement with the configurational adiabat of this EOS is excellent.6 FIG. 15. The excess entropy values plotted against the densities of the state points on the configurational adiabat traced outfor the single-component LJ system starting from the triple point ( ρ = 0 . T = 0 . h = 0 . 04. The valuesare zoomed in to see the deviation from the average value, the black dotted line. [1] C. Alba-Simionesco, A. Cailliaux, A. Alegria, and G. Tarjus, “Scaling out the density dependence of the alpha relaxationin glass-forming polymers,” Europhys. Lett. , 58–64 (2004).[2] C. M. Roland, S. Hensel-Bielowka, M. Paluch, and R. Casalini, “Supercooled dynamics of glass-forming liquids andpolymers under hydrostatic pressure,” Rep. Prog. Phys. , 1405–1478 (2005).[3] E. R. Lopez, A. S Pensado, J. Fernandez, and K. R. Harris, “On the Density Scaling of pVT Data and Transport Propertiesfor Molecular and Ionic Liquids,” J. Chem. Phys. , 214502 (2012).[4] K. Adrjanowicz, M. Paluch, and J. Pionteck, “Isochronal superposition and density scaling of the intermolecular dynamicsin glass-forming liquids with varying hydrogen bonding propensity,” RSC Adv. , 49370 (2016).[5] D. Gundermann, U. R. Pedersen, T. Hecksher, N. P. Bailey, B. Jakobsen, T. Christensen, N. B. Olsen, T. B. Schrøder,D. Fragiadakis, R. Casalini, C. M. Roland, J. C. Dyre, and K. Niss, “Predicting the density–scaling exponent of aglass–forming liquid from Prigogine–Defay ratio measurements,” Nat. Phys. , 816–821 (2011).[6] D. Kivelson, G. Tarjus, X. Zhao, and S. A. Kivelson, “Fitting of Viscosity: Distinguishing the Temperature DependencesPredicted by Various Models of Supercooled Liquids,” Phys. Rev. E , 751–758 (1996).[7] N. Gnan, T. B. Schrøder, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, “Pressure-energy correlations in liquids. IV.“Isomorphs” in liquid phase diagrams,” J. Chem. Phys. , 234504 (2009).[8] J. C. Dyre, “Hidden scale envariance in condensed matter,” J. Phys. Chem. B , 10007–10024 (2014).[9] Y. Rosenfeld, “Relation between the transport coefficients and the internal entropy of simple systems,” Phys. Rev. A ,2545–2549 (1977).[10] J. C. Dyre, “Perspective: Excess-entropy scaling,” J. Chem. Phys. , 210901 (2018).[11] T. B. Schrøder and J. C. Dyre, “Simplicity of condensed matter at its core: Generic definition of a Roskilde-simple system,”J. Chem. Phys. , 204502 (2014).[12] U. R. Pedersen, N. P. Bailey, T. B. Schrøder, and J. C. Dyre, “Strong pressure-energy correlations in van der Waalsliquids,” Phys. Rev. Lett. , 015701 (2008).[13] N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrøder, and J. C. Dyre, “Pressure-energy correlations in liquids. I. Resultsfrom computer simulations,” J. Chem. Phys. , 184507 (2008).[14] D. M. Heyes and A. C. Branka, “Physical properties of soft repulsive particle fluids,” Phys. Chem. Chem. Phys. , 5570–5575 (2007).[15] N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrøder, and J. C. Dyre, “Pressure-energy correlations in liquids. II.Analysis and consequences,” J. Chem. Phys. , 184508 (2008).[16] T. B. Schrøder, N. Gnan, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, “Pressure-energy correlations in liquids. V.Isomorphs in generalized Lennard-Jones systems,” J. Chem. Phys. , 164505 (2011).[17] T.-J. Yoon, M. Y. Ha, E. A. Lazar, W. B. Lee, and Y.-W. Lee, “Topological extension of the isomorph theory based onthe shannon entropy,” Phys. Rev. E , 012118 (2019).[18] A. K. Bacher, T. B. Schrøder, and J. C. Dyre, “The EXP pair-potential system. I. Fluid phase isotherms, isochores, andquasiuniversality,” J. Chem. Phys. , 114501 (2019).[19] A. K. Bacher, U. R. Pedersen, T. B. Schrøder, and J. C. Dyre, “The EXP pair-potential system. IV. Isotherms, isochores,and isomorphs in the two crystalline phases,” J. Chem. Phys. , 094505 (2020).[20] T. S. Ingebrigtsen, T. B. Schrøder, and J. C. Dyre, “Isomorphs in model molecular liquids,” J. Phys. Chem. B ,1018–1034 (2012).[21] D. Fragiadakis and C. M. Roland, “Intermolecular distance and density scaling of dynamics in molecular liquids,” J. Chem.Phys. , 204501 (2019).[22] K. Koperwas, A. Grzybowski, and M. Paluch, “Virial-potential energy correlation and its relation to the density scalingfor quasi-real model systems,” (2020), arXiv:2004.04499 [cond-mat.soft].[23] T. S. Ingebrigtsen and H. Tanaka, “Effect of size polydispersity on the nature of Lennard-Jones liquids,” J. Phys. Chem.B , 11052–11062 (2015).[24] D. E. Albrechtsen, A. E. Olsen, U. R Pedersen, T. B. Schrøder, and J. C. Dyre, “Isomorph invariance of the structureand dynamics of classical crystals,” Phys. Rev. B , 094106 (2014).[25] T. S. Ingebrigtsen, J. R Errington, T. M. Truskett, and J. C. Dyre, “Predicting how nanoconfinement changes therelaxation time of a supercooled liquid,” Phys. Rev. Lett. , 235901 (2013).[26] A. A. Veldhorst, J. C. Dyre, and T. B. Schrøder, “Scaling of the dynamics of flexible Lennard-Jones chains,” J. Chem.Phys. , 054904 (2014).[27] F. Hummel, G. Kresse, J. C. Dyre, and U. R. Pedersen, “Hidden scale invariance of metals,” Phys. Rev. B , 174116(2015).[28] L. Friedeheim, J. C. Dyre, and N. P. Bailey, “Hidden scale invariance at high pressures in gold and five other face-centered-cubic metal crystals,” Phys. Rev. E , 022142 (2019).[29] A. A. Veldhorst, T. B Schrøder, and J. C. Dyre, “Invariants in the Yukawa system’s thermodynamic phase diagram,”Phys. Plasmas , 073705 (2015).[30] P. Tolias and F. L. Castello, “Isomorph-based empirically modified hypernetted-chain approach for strongly coupled Yukawaone-component plasmas,” Phys. Plasmas , 043703 (2019).[31] T. B. Schrøder, U. R. Pedersen, N. P. Bailey, S. Toxvaerd, and J. C. Dyre, “Hidden Scale Invariance in Molecular van derWaals Liquids: A Simulation Study,” Phys. Rev. E , 041502 (2009). [32] T. S. Ingebrigtsen, L. Bøhling, T. B. Schrøder, and J. C. Dyre, “Thermodynamics of condensed matter with strongpressure-energy correlations,” J. Chem. Phys. , 061102 (2012).[33] A. K. Bacher, T. B. Schrøder, and J. C. Dyre, “The EXP pair-potential system. II. Fluid phase isomorphs,” J. Chem.Phys. , 114502 (2018).[34] D. M. Heyes, D. Dini, L. Costigliola, and J. C. Dyre, “Transport coefficients of the Lennard-Jones fluid close to thefreezing line,” J. Chem. Phys. , 204502 (2019).[35] A. Sanz, T. Hecksher, H. W. Hansen, J. C. Dyre, K. Niss, and U. R. Pedersen, “Experimental evidence for a state-point-dependent density-scaling exponent of liquid dynamics,” Phys. Rev. Lett. , 055501 (2019).[36] R. Casalini and T. C. Ransom, “On the experimental determination of the repulsive component of the potential from highpressure measurements: What is special about twelve?” J. Chem. Phys. , 194504 (2019).[37] L. Bøhling, T. S. Ingebrigtsen, A. Grzybowski, M. Paluch, J. C. Dyre, and T. B. Schrøder, “Scaling of viscous dynamicsin simple liquids: Theory, simulation and experiment,” New J. Phys. , 113035 (2012).[38] J. D. Weeks, D. Chandler, and H. C. Andersen, “Role of repulsive forces in determining the equilibrium structure ofsimple liquids,” J. Chem. Phys. , 5237–5247 (1971).[39] J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids: With Applications to Soft Matter , 4th ed. (Academic, NewYork, 2013).[40] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications (Oxford), 1987).[41] N. P. Bailey, T. S. Ingebrigtsen, J. S. Hansen, A. A. Veldhorst, L. Bøhling, C. A. Lemarchand, A. E. Olsen, A. K. Bacher,L. Costigliola, U. R. Pedersen, H. Larsen, J. C. Dyre, and T. B. Schrøder, “RUMD: A general purpose molecular dynamicspackage optimized to utilize GPU hardware down to a few thousand particles,” Scipost Phys. , 038 (2017).[42] A. de Kuijper, J. A. Schouten, and J. P. J. Michels, “The melting line of the Weeks–Chandler–Anderson Lennard-Jonesreference system,” J. Chem. Phys. , 3515–3519 (1990).[43] A. Ahmed and R. J. Sadus, “Phase diagram of the Weeks-Chandler-Andersen potential from very low to high temperaturesand pressures,” Phys. Rev. E , 061101 (2009).[44] U. R. Pedersen, L. Costigliola, N. P. Bailey, T. B Schrøder, and J. C. Dyre, “Thermodynamics of freezing and melting,”Nat. Commun. , 12386 (2016).[45] R. Casalini and T. C. Ransom, “On the pressure dependence of the thermodynamical scaling exponent γ ,” Soft Matter , 4625–4631 (2020).[46] T. Maimbourg, J. C. Dyre, and L. Costigliola, “Density scaling of generalized Lennard-Jones fluids in different dimensions,”SciPost Phys. , 90 (2020).[47] D. M. Heyes, D. Dini, and A. C. Branka, “Scaling of Lennard-Jones liquid elastic moduli, viscoelasticity and otherproperties along fluid-solid coexistence,” Phys. Status Solidi (b) , 1514–1525 (2015).[48] F. L. Castello, P. Tolias, and J. C. Dyre, “Testing the isomorph invariance of the bridge functions of Yukawa one-componentplasmas,” J. Chem. Phys. , 034501 (2021).[49] J. C. Dyre, “Simple liquids’ quasiuniversality and the hard-sphere paradigm,” J. Phys. Condens. Matter , 323001 (2016).[50] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes: the art of scientific computing ,3rd ed. (Cambridge University Press, 2007).[51] H. Flyvbjerg and H. G. Petersen, “Error estimates on averages of correlated data,” J. Chem. Phys. , 461–466 (1989).[52] J. Kolafa and I. Nezbeda, “The Lennard-Jones fluid: An accurate analytic and theoretically-based equation of state,” FluidPhase Equilib.100