Density scaling of generalized Lennard-Jones fluids in different dimensions

Liquids displaying strong virial-potential energy correlations conform to an approximate density scaling of their structural and dynamical observables. These scaling properties do not extend to the entire phase diagram, in general. The validity of the scaling can be quantified by a correlation coefficient. In this work a simple scheme to predict the correlation coefficient and the density-scaling exponent is presented. Although this scheme is exact only in the dilute gas regime or in high dimension d, the comparison with results from molecular dynamics simulations in d = 1 to 4 shows that it reproduces well the behavior of generalized Lennard-Jones systems in a large portion of the fluid phase.


Introduction
The past 20 years have developed an increasing interest in the so-called density scaling approach. Starting from the experiments of Tölle, Dreyfus and Alba-Simionesco [1][2][3], it has been found that many liquids at different state points (ρ, T ) in their phase diagram exhibit a constant relaxation time along lines with fixed ratio ρ γ /T . T and ρ are the temperature and density, respectively, and γ is the so-called density-scaling exponent. This behavior has been often interpreted as the result of the repulsive part of the interaction potential giving the dominant contribution of the dynamics, allowing for mapping to an inverse-power-law (IPL) potential: v IPL (r) = ε(σ/r) λ [4][5][6][7][8][9][10]. Knowledge of the density-scaling exponent at a given state point allows one to transfer information about the dynamics of the system to other state points, thus drastically reducing the amount of experiments needed to measure properties of a system in its phase diagram. It allows as well predictions of the dynamics in regimes that are difficult to probe experimentally. For an IPL system the density-scaling exponent is related to the IPL exponent λ through the relation γ = λ/d, where d is the spatial dimension. It is therefore constant throughout the entire phase diagram. As such, the IPL potential is paradigmatic as the density scaling holds exactly. Then one may use this scaling approach for other systems. This has been successful in many situations [11][12][13] but is undermined by several weak points: 1. The very idea of mapping of a real material's interaction to an IPL pair potential is not satisfying as real materials display gas-liquid and gas-crystal coexistence not present in the phase diagram of the IPL system; 2. There are several evidences [13][14][15][16] that density scaling does not work for all systems. For real systems it simply cannot apply in the entire phase diagram and the density scaling exponent is state-point dependent [17][18][19]; 3. It has been argued by Bøhling et al. [20] that the standard assumption that the repulsive and attractive parts of the potential play separate roles is difficult to justify.
A formulation of density scaling that eliminates the above-mentioned problems and contradictions is provided by the isomorph theory [11,21]. In this framework, the density-scaling exponent γ is generally state-point dependent. The invariance of static and dynamic properties along the lines of constant excess entropy (the so-called isomorphs), which defines the density-scaling exponent, is related to a scale invariance of the potential-energy hypersurface [22]. According to isomorph theory γ can be obtained from equilibrium simulations at the state point in question [11]; for a real material it can also be determined by measuring several quantities at the state point [17,23].
The density-scaling exponent is not the only relevant state-point dependent quantity in the isomorph theory: it is paired with the virial potential-energy correlation coefficient R ∈ [−1, 1], the value of which indicates whether or not density scaling is satisfied in the proximity of the state point in question. It can be shown that perfect invariance of structure and dynamics along isomorphs is ensured if the correlation coefficient is equal to unity, R = 1 [11]. This identity holds only for Euler homogeneous potentials [11], as IPL for instance, while for a real system R = 1 never applies. Therefore, a threshold value of 0.9 has been chosen; whenever R > 0.9, density scaling is expected to work well. Both quantities are defined through equilibrium correlations of fluctuations: in which U denotes the potential energy, W is the virial 1 (i.e. , the excess part of P V with respect to the ideal gas, where P is the pressure and V the volume), and ∆ is the instantaneous deviation of a given quantity from its thermodynamic equilibrium value. The quantities γ and R are well defined in the entire phase diagram, not only in the region where density scaling works. While γ can be determined in experiments [6,10,17,19,23], this is not the case for R. A crucial prediction for density scaling would therefore be to relate the condition of R crossing some threshold value with properties of the state-point dependence of γ, or even to compute in a simple manner the state-point dependence of R. A first step in this direction was taken by Friisberg et al. [25], who showed that for generalized LJ potentials with exponents (2n, n) -defined below, in Eq. (3) -, the value R = 0.9 roughly identifies the state points with highest density-scaling exponent γ. A convenient way of showing this is to plot γ(ρ, T ) versus R(ρ, T ). It was also shown that, in this system, γ is a unique function of R: γ(ρ, T ) = F (R(ρ, T )), to a good approximation. Such a relation allows one to better understand the density scaling throughout the phase diagram. We emphasize that such a plot of γ versus R is related to the system not being a simple IPL system since in the IPL case all state points map onto the single point (R, γ) = (1, λ/d).
In the present work, generalized Lennard-Jones systems are studied theoretically and numerically, extending the results of Ref. [25]. The density-scaling quantities mentioned above are derived from excess thermodynamic observables, which are naturally expressed in the language of the virial expansion [24]. We thus devise a simple low-density approximation for the density-scaling exponent and the virial-potential energy correlation coefficient, and compare it to computer simulations, which show very good agreement. Then, as this approximation becomes also exact in another limit of infinite dimension, we connect these results to the recent finding that isomorph invariance is exactly achieved for a large class of potentials in d → ∞, beyond the Euler-homogeneous ones like IPLs, in the liquid and glass phases [26]. The generalized Lennard-Jones potentials, albeit very important in practice, do not belong to this class; but computer simulations in Ref. [27] provided evidence that density-scaling becomes more robust increasing the number of dimensions, including close to the liquid-gas coexistence. We thus clarify and extend the results of Refs. [26,27]. We observe a convergence to the large-d limit, starting already from dimension d = 2.
The structure of the paper is the following: in the next section we introduce the system studied and summarize past observations relevant to the present article. Then in the third section we present the analytic approximation which is confronted in the fourth section to  Figure 1: Different generalized LJ(m , n ) pair potentials. (a) m = 2n for several choices of the value of n = 4, 9, 12, 18, 100. Upon increasing n the potential becomes sharper and sharper, approaching a sort of sticky-sphere limit (similar in spirit to the discontinuous version of Baxter [28]). (b) m = Xn for X = 2, 3, 4, 5, 50 and fixed n = 6. For large X the attractive tails of the potential tend to a limiting shape, while the repulsive part approaches a hard wall at r/σ = 1. molecular dynamics simulations in two, three, and four dimensions, as well as in a wide range of temperature and density. A concluding discussion ends the paper.

Previous results on generalized Lennard-Jones potentials
Following Friisberg et al. [25], the class of generalized Lennard-Jones potentials will be studied in this work. It is defined as follows: and will be denoted as a LJ(m , n ). This class of potentials is also known as Mie potentials [29], providing a generalization of the standard LJ potential that is widely used in the liquid-state literature. The definition Eq. (3) implies that the minimum of the LJ(m , n ) potential is located at r = σ where the potential's value is − . This normalization is different from the standard definition of the 12-6 LJ potential; it ensures that the minimum does not shift when the exponents vary. The conclusions in the following do not depend on this particular choice. The units of length and energy, respectively σ and ) are set to 1. Incidentally, the standard 12-6 LJ in three dimensions corresponds to a simple rescaling of these units in the LJ (12,6) potential. Temperature is measured in energy units (k B = 1). In this work the exponents m and n are changed independently, in contrast to Ref. [25]. The LJ(m , n ) potential is displayed in Fig. 1 for several choices of the integers (m , n ). In the following we consider only integer values, but the conclusions do not depend on this limitation. For convenience we define the ratio The value X = 1 is excluded because the LJ(n , n ) potential is not defined. Note that at low density for fixed temperature (or high temperature for fixed density) the repulsive IPL term dominates, implying R → 1, γ → m /d. Although in general the correlation coefficient is below unity, it was shown in Refs. [11,14,30] that the standard 12-6 LJ potential is strongly correlating (i.e. R > 0.9) in the region above the freezing line in the (ρ, T ) phase diagram. Furthermore, in Ref. [25] the LJ(2n , n ) potential was studied at many state points in dimensions two to four, displaying either gas, liquid/fluid, crystal phases or coexistence between these. Away from strongly-correlating regions, it was found that the simple relation holds to a good approximation when varying dimensions d or exponent n . As pointed out in the introduction, such a relation is most welcome in order to assess quickly isomorph invariance, given by the R value, from experimental measure of γ. An analytically manageable formula for R and γ would be useful for better understanding why and how these two quantities behave and correlate. The simple dimensional dependence, Eq. (5), similar to the prediction of an IPL system with exponent 3n [30], is intriguing. It was shown in Ref. [26] that in the dense liquid regime, a sizeable class of potentials are effectively IPLs in large dimension and as such display isomorph invariance. They are therefore of little use for understanding the variations of (R, γ) observed in real systems. This is, however, not the case for the LJ(m , n ) potentials [26], which retain their 'non-IPL' behavior even in high dimensions. For this limit to be well-defined, one needs to scale the exponents with the dimension and thus consider LJ(md, nd) pair potentials 2 [26,[34][35][36]. In this case, one expects the large-dimensional limit to exhibit scaling properties closer to the finite-dimensional systems than the potentials studied in Ref. [26], while still being considerably simpler. Molecular dynamics simulations in Ref. [27] showed that the standard 12-6 LJ fluid above (in density or temperature) the liquid-vapor critical point manifests increasingly better isomorph invariance as d goes from 2 to 4, i.e. increasing R. This observation tends to broaden the conclusions of Ref. [26], which hold in the liquid and glassy regimes, to the fluid region close to the liquid-vapor critical point, which cannot be described by the arguments developed in Ref. [26]. Indeed they rely on a truncation of the virial expansion to the second virial coefficient, valid for dense liquid regimes in high d [37][38][39], whereas the third virial coefficient cannot be neglected close to the liquid-vapor critical point [24]. It turns out that, in the infinite-dimensional limit, the liquid-vapor critical point regime and the dense (possibly supercooled) liquid regimes occur at density ranges exponentially separated in d [40], displaying different behaviors. We point out that in these simulations (Ref. [27]), the exponents were not scaled with d, which is mandatory to compare with high-d theory. We take this into account in the present work, and consequently investigate the effect of dimensionality in light of a systematic comparison to the mean-field prediction for LJ(md, nd) interactions.

Analytic expressions derived from the virial expansion
In order to get an analytic handle on the quantities γ and R, we study here their lowest order in the virial expansion of liquid-state theory [24]. The reason is threefold: First, the quantities studied here are excess thermodynamic observables, for which this framework has been designed. One expects that this expansion is well suited to the study of gases and dilute liquids as the density is small, while it could perform poorly for dense liquids in low dimensions where other perturbation strategies or re-summations work better [24]. Second, it is also a right choice for the mean-field limit of large dimension [36][37][38][39]41], allowing us to observe the impact of the parameter d. Third, the observed density dependence of the (R, γ) diagram is found to be very mild, irrespective of the thermodynamic phase [25], and one may thus hope that the virial approximation provides a good starting point.

The lowest order in the virial expansion in any dimension
In this section we give the first-order virial expansion of certain thermodynamic fluctuation averages. The number density ρ = N/V is regarded to be a small parameter. The lowest order of any observable is then its ideal gas value; we are interested in the first non-trivial small-density correction. The equilibrium fluctuations computed below are the ones entering in the definition of the virial-potential energy correlation coefficient R and the density-scaling exponent γ, i.e. the fluctuations ∆U ∆W , (∆U ) 2 and (∆W ) 2 (see Eq. (1)). We shall here just give the main steps as a detailed derivation is found in [26,Appendix]. The first two fluctuations are rewritten using the definition of canonical equilibrium averages as follows, where β = 1/k B T is the inverse temperature, k B is the Boltzmann constant and F is the Helmholtz free energy. The average virial is related to the pressure by the virial equation of state [24]. The expansion of the latter and the free energy follows from a standard computation [24,38,39], with v(r) the radial pair potential: Combining Eqs. (6),(7) we get ∆W ∆U and (∆U ) 2 . The derivation of (∆W ) 2 from the lowest-order virial expansion of the two-point distribution function is slightly more involved [26,Appendix], but can be guessed from the relationship between W and U , Eq. (2). The final result, for all three fluctuations, is where we performed an integration by part in the last line of Eq. (8). We can now give the following first-order virial expressions from the definition Eq. (1): These expressions apply for any isotropic pair-potential liquid.
A check of the validity of Eq. (9) can be obtained by considering the IPL potential v IPL (r) = ε(σ/r) λ for which we already know the result. Since rv IPL (r) = −λv IPL (r), IPL systems have perfect virial-potential energy correlations (R = 1). Using spherical coordinates and standard properties of the Euler Gamma function one arrives at which is the general result for IPL systems in d dimensions [11], as mentioned in the introduction. It must be recovered from Eq. (10) because this general results holds in particular at low densities where the virial approximation is exact.
In the case of the LJ(md, nd) potential in d dimensions, i.e. Eq. (3) with m = md and n = nd, the quantities R and γ can be simplified further from Eq. (9). This result is obtained using similar considerations to the IPL case, e.g.
in which X = m/n, V d (r) is the volume of an hypersphere of radius r in d-dimensional space, and Ω d is the solid angle in such a space given by Note that the condition n > 1 discussed in footnote 2 implies that the y 1−1/n divergence at small y (long distances) is integrable (as well as y −1/n ), whereas our restriction X > 1 implies integrability at large y (short distances). Similar considerations for the other integrals involved lead to the following expressions for R and γ in which y = z = (σ/r) nd . These equations establish the lowest order in the small-density expansion. We note that: • R and γ at this lowest order are independent of the density (albeit keeping in mind that its validity is guaranteed only for low enough density). Only the temperature and the details of the potential enter as parameters. This is a direct consequence of the fact that they are expressed as ratios of extensive quantities. Indeed, any extensive observable O scales in the thermodynamic limit as O ∝ N ∝ ρ and thus O = O(ρ) in the virial expansion.
• R and γ are also independent of the space dimensionality. This is due to the fact that the potential is built with IPLs whose exponents are proportional to d. The numerical computation of these quantities in any dimension is therefore straightforward.
• The explicit dependence upon n is rather mild, unlike the X dependence.
These facts imply that a comparison to simulation data in any dimension is easily achieved. This is yet another reason to scale LJ(md,nd) with dimension d in this way, in addition to ensuring the existence of a well-defined thermodynamic limit. Besides, we expect any thermodynamic observable constructed as a ratio between extensive quantities involving only the pair potential v md,nd (r) to exhibit the properties mentioned in the above first two points. As the dependence upon n is quite mild, one can simplify the analytical expressions by considering the large n limit at fixed X. The three different integrals become The corresponding values of R and γ are Both quantities are of order 1 with respect to n, which is why we divided γ by a quantity proportional to n; we chose m + n for better comparison with Ref. [25], as the linear relation found there reads γ = (m + n)R.
Finally, in the case of the LJ(2nd, nd) potential (X = 2), the integrals become Gaussian and one arrives at the very simple expressions R is a strictly increasing function of temperature, so that γ(T ) can be parametrized instead by R and the function γ(R) is well defined. Its curve for X = 2 and any n from Eq. (13) can as well be drawn by a parametric temperature plot. It has a non-trivial shape that is well described by the above n → ∞ expression: indeed already for n = 2, the curve γ m+n (R) has the same shape and the difference with the n → ∞ curve is always below 10%. In Sec. 4 we will confirm this statement comparing these analytical results to the simulation data.

The high-dimensional limit
The considerations of the last section are independent of the dimension d but hold only in the dilute regime of the liquid phase. The virial truncation is, however, exact in the whole (possibly supercooled) liquid phase for mean-field situations, such as when the spatial dimension goes to infinity [36][37][38][39]41], or in the infinite-range Mari-Kurchan model considering random shifts of the particles in any dimension [42]. We thus expect smaller deviations from the exact d → ∞ reference situation as dimension increases, but one does not know a priori how large the dimension must be to ensure satisfactory convergence to the mean-field prediction. This is the aim of a following section, Sec. 4.2.
We now mention what can be expected from Eq. (13), focusing on the limit d → ∞. Let us first restrict to the liquid phase. For large d the physics of the LJ(md, nd) potential (3) is as follows [26]. Both thermodynamics and dynamics are dominated by a fluctuating region of order 1/d around a length scale r * defined by requiring βv md,nd (r * ) to be of order unity with respect to d. The length scale r * plays the role of an effective particle diameter. At temperatures exponentially large in d, r * < σ and the attractive IPL term is exponentially suppressed. The system is then dominated by the repulsive IPL term, and it thus has exact isomorphs in this regime: R = 1, γ = m (compare Eq. (10)). Conversely, for exponentially small temperatures, the system will be dominated by the attractive IPL term, which has a non-stable thermodynamic limit. The non-trivial regime, closer to the finite-dimensional system, is that of O(1) temperature in which r * = σ. In this case both terms in the potential compete and the potential has the following behavior • for r < σ, the interaction is effectively hard core: βv md,nd (r) → ∞ • in the 1/d region around σ, i.e. r = σ(1 + r/d) with r a reduced distance of order unity, the potential is a sum of two competing repulsive and attractive exponentials • for r > σ, there is effectively no interaction and βv md,nd (r) is exponentially vanishing.
As in finite dimension, both attractive and repulsive terms contribute in this regime. Consequently no exact isomorphs are found, i.e., no rescaling of the density and the temperature can make the free energy and dynamics invariant [26]. One has R < 1 except in the infinitetemperature limit where R → 1, which is the case in which only the repulsive IPL dominates.
The precise values and evolution of R will be provided in the next section. From e.g. Eq. (16), it is clear that for low temperature (R, γ) → (0, 0), and R increases monotonically to 1 at high temperature. Yet in the next section, there exist negative values of R and γ, or large values of γ at small R, as outcome of simulations in d 4. As detailed in the next section, they are due to state points in a liquid-vapor or solid-vapor coexistence phase. Such values are absent from Eqs. (9). Indeed, even in large dimensions the validity of the latter equations is restricted to the uniform liquid phase, and Eqs. (9) do not necessarily extend to the crystal or liquid-vapor coexistence regimes 3 . The reason is that in deriving the above equations we implicitly assumed that the thermodynamic density profile is uniform. On the one hand, while crystalline phases and the liquid-crystal transition have been studied for d = 4, 5, 6 in Refs. [43,44], precise descriptions of large-d crystalline phases are not yet possible as the equilibrium crystalline configurations in large d are not known. This is due to the daunting geometry of spherical packings in high d [45] and because numerically nucleating the crystal phase through compression of the liquid has proven extremely difficult when the dimension is larger than three [44,46]. On the other hand, following a Landau approach for large d, the liquid-vapor coexistence regime may be studied by including one additional order in density in the virial expansion [24]. But then the density window considered is exponentially separated from the one of the liquid/fluid phase [38,39], which makes comparison to the finite-d system harder. Indeed, the regime close to the critical point has been studied by Mon and Percus [40] for a square-well potential, where one can analytically extract a criticalpoint density whose exponential dependence is given through the effective packing fraction 4 ϕ c ≈ (4/ √ 3) −d 2 −d , much lower than the (dense) liquid phase scaling ϕ = O(d/2 d ) [36,41]. The critical-point temperature is moreover not of order unity, as T c = O(1/d). We expect the same dimensional scaling of the critical density to hold for a LJ(md, nd) potential.

Molecular Dynamics simulations in dimension one to four
In this section the expressions previously derived from the virial expansion are compared with Molecular Dynamics simulation data for different LJ potentials in 2d, 3d, and 4d. These results include the standard 12-6 LJ potential in three dimension. The variation with the main parameters -the density, the potential exponents, and the dimension -will be analyzed. We find qualitative (if not quantitative) agreement with the analytic expressions of R and γ, Eqs. (13), in a large density regime.
Two different Molecular Dynamics codes were put to work. The 2d and 3d simulations employed the GPU-based Roskilde University Molecular Dynamics (RUMD) code [47], while the 1d and 4d simulations used an ad hoc CPU-based code (more details can be found in Ref. [48]). All simulations were performed in the N V T ensemble in which temperature was controlled using a standard Nosé-Hoover thermostat [49]. The time step ∆t is kept constant in reduced units (∆t) when exploring the phase diagram, i.e. gets rescaled at each state point via the only potential-independent time unit of the system: ∆t = ρ −1/d (M/k B T ) 1/2 ∆t, with M the particle mass. In particular this scaling ensures that the ballistic part of the mean-squared displacement collapses onto a single master curve for all state points. ∆t values are given in Tab. 1 together with the system size N and the number of simulation steps. The potentials studied in this work were cut-off at distance r cut = 2.5σ. The choice of the parameters in Tab. 1 is mainly dictated by running time limitations. Our ad hoc code is less optimized than RUMD and simulations cannot be run as extensively. We checked that the duration of the simulations in 4d is not an issue to calculate accurately thermodynamic quantities. In 4d the system size is chosen to be the smallest necessary to reliably study the highest-density state points (ρ = 1.5). In 3d simulations, two different system sizes are listed to ensure the absence of system size dependence. Similarly in 1d the results are reported for N = 40000 particles, but other simulations were ran for smaller system size (N = 500).
At each state point (ρ, T ) we compute the correlation coefficient R and the scaling exponent γ (Eq. (1)). Density lies within the interval ρ ∈ [0, 1.5] and temperature T ∈ [0.5, 5] , which is enough to retrieve all range of values of R and γ. A complete list of the state points studied in this work can be found at the data repository of the Glass and Time group (see Acknowledgements).  Table 1: Simulation parameters. The system size N , the reduced time step ∆t, and the number of time steps simulated at each state point for all four dimensions considered numerically in this work. For the LJ(300, 6) potential, we resorted to a smaller reduced time step ∆t = 10 −4 because of the potential steepness at small distances.
Some of the state points exhibits phase coexistence; these are associated with low values of R [12]. The density-scaling exponent γ is not uniquely defined when different phases coexist. We find large γ and small R in the crystal-vapor coexistence, whereas in the liquid-vapor coexistence both R and γ are small (close to zero or negative). As mentioned in Sec. 3.2, state points in the coexistence region cannot be interpreted with the analytical equations obtained in this work because they were derived based on the implicit assumption of a single, isotropic phase. We therefore left out these state points in the following analyses except in Fig. 2 (c), where they are left as an example.

Density dependence of the γ(R) relationship
In Fig. 2 we study the (R, γ) diagram of the LJ(4d, 2d) (n = X = 2) for d = 1 -4. As will become clear from Fig. 3 below, the shape for this particular potential is representative of a wide range of values of n and X. In fact, it corresponds to the standard LJ potential for d = 3 with rescaled units of length (thus density) and temperature (see Sec. 2). Apart from a rescaling of the isochore value, the curves γ(R) are thus the same for both potentials , the potential is v 12,6 (r), which can be mapped to the standard LJ potential by scaling density and temperature. In the figure data for both the standard 12-6 LJ (labelled as "std") and the generalized LJ(12,6) are presented. The unit length for these two potentials is different, resulting in different definitions for the density. The lowest temperature on each isochore is either along the liquid-gas or along the liquid-solid coexistence curves in the case of the standard 12-6 LJ while for the generalized LJ(12,6) the same densities and temperatures as in d = 1, 2 are considered. In the case of the standard LJ potential the phase diagram is therefore explored in a more comprehensive manner: the highest temperature on each isochore is roughly 200 times the starting one and thus well above T = 5. At such a high temperature the density scaling exponent γ is close to 4 and those state points are consistent with the IPL limit (R, γ) → (1, 4). The similarity between the shape of the data in this figure and in (a),(b),(d) shows that it is not needed to consider such a wide region of the phase diagram to get a satisfying overview of the behavior of the γ(R) curve. (d) d = 4, the potential is v 16,8 (r). in d = 3 when exploring all the many different state points plotted below. As discussed in Ref. [25], these curves are similar for all isochores and dimensions. At low densities the curves are very close to the prediction of Eq. (13) for n = 2 and X = 2, as expected. As the density is increased, the curve stays still close to the analytic prediction, although the values of the density-scaling exponent become slightly higher and the state points are more confined into the region around R = 1. This is because we do not take into account liquid-vapour coexistence state points, which have a low temperature. We are then left with higher temperatures, for which the repulsive IPL term is more dominant, resulting in a stronger correlation R. Note that the attractive term cannot be ignored (being γ different from 4). Also, for d = 1, 2, 3 ordered phases are found at high densities for which the correlation coefficient R is almost unity. For a better quantitative agreement at high densities, the discrepancy can be reduced through a simple rescaling procedure, with the knowledge of the maximum density-scaling exponent, denoted by γ peak . This is shown in Fig. 6 where the data of Fig. 2 are shown for d = 2 and d = 3 rescaled by γ peak . These two dimensions are displayed as they are the most relevant and because only in these cases it is possible to clearly identify a maximum value in the (R, γ) curve for some isochores. Near-quantitative agreement for all state points is obtained between the analytic virial expression and the numerical values of γ(R)/γ peak .
The dependence of the previous diagrams on the exponents n and X, at fixed dimension d = 3, is explored in Fig. 3. In Fig. 3(a) the potential LJ(2nd, nd) is studied varying n from 4/3 to 6 in 3d, showing little influence of the value of n on the relation between γ and R, as expected from Eq. (13). In Figs. 3(b)(c)(d) the influence of the ratio X = m/n is studied: the results for the LJ(Xnd, nd) potential with (n, d) = (2, 3) are shown for many values of the ratio X. For moderate values, almost no difference is seen, and one must reach values of X of order 50 to see a qualitative difference in the shape. At very high X ≈ 10 3 , the curve converges to a parabolic shape γ = mR 2 as argued in Appendix B, as the non-parabolic part of the curve gets shifted to the R = 1 region as X increases and gradually disappears. There seems to be no distinction in the (R, γ) diagram between the large-X LJ(2Xd, 2d) and a continuous 'sticky-sphere' potential LJ(Xnd, nd) with n and X → ∞, i.e. for large X the tail shape does not matter much in the liquid regime (see also Fig. 1).

Dimensional dependence of the γ(R) relationship
We next examine how fast the curves obtained above at finite dimension converge towards the analytical prediction of Eqs. (9) - (16). In order to investigate the deviations from the infinite-dimensional (R, γ) diagram, we quantify the deviation by subtracting the d → ∞ curve from the finite-d one, γ d (R) − γ ∞ (R) in which γ d (R) is isochore dependent and γ ∞ (R) is obtained from Eq. (13) with n = X = 2. When increasing d, we must compare different isochores in different dimensions, and we shall now detail several strategies employed so as to choose the isochores in a physically meaningful way.
We first plot in Fig. 4 the deviations from d = 1 to 4 at fixed density. This choice is natural as the interval of values of the virial-potential energy correlation coefficient plotted in this way is roughly the same for all the dimensions considered. For all densities the deviations decrease as d increases from 1 to 4. This comparison is somehow not satisfactory because fixing the density means comparing different physical situations. The volume occupied by a particle, ∼ V d (σ/2), where the interaction range of the potential σ is interpreted as a particle diameter, decreases monotonically with d and therefore the system gets more diluted as d it is shown that when X becomes large, a more natural scaling for γ is γ/m. Note that in (b) the R = 1 endpoint at infinite temperature varies slightly with X, as γ/(m + n) = X/(X + 1); this is not the case when γ is scaled as suggested by the large X limit. At large X, the shape of the curve also changes. The peak is gradually shifting to higher R and disappears for X = 50. This is shown both by MD simulations (c) and from the virial expressions of Eqs. (13) for n = 2 (d). The curves decrease monotonically towards a parabola γ = mR 2 (Appendix B).
increases. A plausible interpretation of the reduced deviations with increasing d in Fig. 4 is that the system in d + 1 dimensions appears to be less dense than its counterpart in dimension d. As a result it could be that the virial truncation performs better not because the meanfield d → ∞ approximation improves, but primarily as an indirect 'low-density' effect. To The horizontal dashed line at the origin of the vertical axis corresponds to the exact d → ∞ curve. Each plot from (a) to (e) is given for a fixed isochore with density varying from ρ = 0.25 to 1.25 in steps of 0.25. The higher the density, the closer state points are from the R = 1 boundary. Simulation data before subtraction are displayed in Fig. 2. compensate for this, the density must increase when comparing the d to d + 1 data, as has been observed for the liquid-crystal transition [27,43] or the liquid-glass transition [44,46].  Figure 6: The data sets of Fig. 2 for which it is possible to identify a maximum value for γ (termed γ peak ) are here shown again with the same color coding, namely (a) d = 2 ( Fig. 2(b)) and (b) d = 3 (Fig. 2(c)). For each isochore ρ the density-scaling exponent γ is scaled by γ peak . This rescaling partially accounts for the density corrections not taken into account in a first-order virial expansion.
Note that in the large-dimensional limit, the dense liquid region emerges for densities scaling exponentially in d [35,36,38,39,41].
There is no a priori simple way to compute densities of related physical regimes in different dimensions. In order to get reasonable values, we attempted two different strategies, focusing on fairly dense liquid regimes, since their large-dimensional limit is well understood [26,36,39,41]. The first strategy computes ratios between the dynamical glass transition densities ρ HS d (d) for hard spheres. This is the only potential for which such a transition has been determined in dimensions ranging from 2 to 12 [46,50]; temperature does not influence this scaling. The second strategy amounts to estimating as a function of the dimension the density at which the first-order virial expansion breaks down, providing meaningful density ratios for the regime we are after. This characteristic density, ρ Z (d, T ), is estimated by comparing the first-and second-order virial coefficients, yielding ρ Z (d, T ) ∼ |B 2 (T )/B 3 (T )| in which B 2,3 (T ) are the second and third virial coefficients, respectively [24]. The density ratios between different dimensions are fairly independent of temperature far from liquid-gas coexistence (Appendix A). The dimension-dependent density ratios provided by both methods are listed in Tab. 2.   [46,50]. One-dimensional hard spheres do not exhibit such a slowing down, so for d = 1 we considered the dense regime to occur at the maximal packing fraction (unity). The second method (second line) is based on the calculation of ρ Z (d, T ) (Appendix A).
Both scaling-method outcomes are plotted in Fig. 5. We fixed the density in 4d to be ρ d=4 = 0.75, an intermediate value. This is not so low that meaningful measurable differences with respect to the d → ∞ analytic (R, γ) curve are detected, and not so high so the computed virial-potential energy coefficient spans a large interval (compare Fig. 4). Both scalings yield qualitatively similar results. One observes a convergence to the d → ∞ analytic prediction going from d = 2 to d = 4. The one-dimensional values are however significantly closer to the large-d result than the two-dimensional ones. We conclude that, similarly to the non-scaled plots of Fig. 4, already for low values of dimension, a convergence towards the large-dimensional (R, γ) diagram is observed.

Conclusion
We have introduced a set of simplified equations to compute the virial-potential energy correlation coefficient R and the density-scaling exponent γ, valid in any dimension and for any pair potential in the isotropic liquid phase. They are obtained from a low-density virial expansion. As such they are exact in two limits : the low-density limit (in a given dimension) and the infinite-dimensional limit (for any density), in the isotropic liquid.
We have specialized these results to the case of LJ(m , n ) systems. The interest of such kind of equations is that both R and γ can be computed straightforwardly, through numerical integrations of a few one-dimensional integrals, for any state point in the phase diagram. We showed through molecular dynamics simulations that this approximation applies qualitatively and almost quantitatively to get the (R, γ) diagram, if the system does not phase separatea case in which density scaling does not apply, even approximately. It appears that density corrections are weak for the LJ(m , n ) potentials, irrespective of the value of the exponents or temperature. The analytical shape γ(R) allows then to make the following prediction, robust in the whole fluid phase: if the measured exponent γ decreases from one state point to the other along an isochore, then we are probing a strongly-correlating regime, i.e. a regime where density scaling is satisfied to a good approximation [25]. This observation is useful in practice only for a ratio between respectively the repulsive and attractive exponents X 10.
We do not expect in general the monotonicity of γ to yield a direct indication of good scaling for any potential. Yet, the simplified low-density limit expressions from (9) can be helpful for other potentials, as they provide the relation γ(R), allowing to assess R from the possible measurement of γ [17,23], or the R-γ relation if the latter function is multivalued (this occurs for instance for a potential consisting in a sum of IPL potentials with different exponents). As an example, preliminary data shows that the function γ(R) for WCA potentials [24,51] gives the correct qualitative behaviour while significantly differing in the functional form from the LJ(m , n ) one.
It was recognized in Ref. [26] that perfect density scaling is achieved for many non-trivial potentials -i.e., potentials which are not necessarily Euler homogeneous -in the d → ∞ limit. LJ(m , n ) potentials do not display such a perfect scaling in the high-d limit, and are instead characterized by a (R, γ) diagram not restricted to a single point. We varied the dimension from d = 1 to d = 4 in order to check the convergence to the high-d (R, γ) diagram for intermediate densities where there are corrections to the low-density expressions. Using several possible scalings of densities with dimension leads to the same conclusion of a monotonic shrinking of fluctuations from d = 2 to 4. We interpret this as yet another instance of the large-dimensional being qualitatively and sometimes even quantitatively good [36] for low dimensions.
In this paper we have considered only first-order virial expansions; these apply at low density and/or high dimension. In finite dimensions one expects the first-order approximation to break down for densities at which the next term is relevant, i.e. whenever B 2 (T )ρ 2 ≈ B 3 (T )ρ 3 . This defines a density ρ Z (d, T ) = −B 2 (T )/B 3 (T ). At this density, the compressibility factor becomes βP/ρ Z = 1 + O(ρ 3 Z ), meaning that ρ Z (T ) coincides with the low-density Zeno 5 line [52,53], which in d = 3 lies inside the supercritical region well above the liquid-vapor critical point.
We computed numerically ρ Z (T ) for d = 1 to 4 ( Fig. 7(a)). Both virial coefficients are negative at low T , and become positive above it. B 2 (T ) vanishes for a dimension-independent value, the Boyle temperature T Boyle 3.4 (as can be readily seen from the definition in Eq. (18) from the same manipulations as the ones in Sec. 3.1) [53,54], whereas B 3 (T ) vanishes close to T = 1 depending on the dimension. This explains the observed divergence of ρ Z (d, T ) in this region. The negativity of these coefficients indicates phase separation at small enough density, i.e. here a liquid-vapor coexistence. Indeed from Eq. (18) at small density this negativity implies dP/dρ < 0, which signals a thermodynamic instability. As mentioned in Sec. 4 we wish to stay away from this regime in which our analytical results are no longer justified. The low-density equation of state Eq. (18) shows no sign of phase separation above T 1.7 for any dimension d = 1 to 4. Thus focusing on these higher temperatures, corresponding to the fluid phase above coexistence (at small enough densities), one realizes that all the density ratios at fixed temperature ρ Z (d + 1, T )/ρ Z (d, T ) are approximately constant, compare Fig. 7(b). Note that as we are interested in an order of magnitude for density ratios between different dimensional systems for which the first virial truncation breaks down, the sign of ρ Z (d, T ) does not matter. Consequently, we can extract a temperature-independent meaningful scaling of density by averaging the value of this ratio over the whole temperature range T ∈ [T 0 , 5]. We took T 0 = 1.96 > 1.7 as the choice of T 0 , which does not modify considerably the ratio values while we must at the same time consider enough statistics, as displayed in Figs. 7(c)-(d). For this value of T 0 , indeed, fluctuations of the computed density ratio are below 5% with respect to the average in all dimensions. This procedure provides the numbers indicated in the last line of Tab. 2, which are appreciably above the density ratios defined by the hardsphere dynamical transition densities (first line). The choice of T 0 close to 1.7 maximizes the ratio values with respect to higher temperature (see Figs. 7(b)-(c)). As higher densities are associated with stronger deviation from the virial approximation (compare Fig. 2), this is the most unfavorable situation in order to see smaller deviations to the large-d curve γ ∞ (R); yet we find in Fig. 5 a good convergence for d = 2 → 4 using such ratio values.

B Large-exponent ratio limit
In this Appendix we investigate the large-X limit (X = m/n) of Eq. (13). As argued in Sec. 3.1, the n dependence is rather mild, so that we send first n → ∞. From Eq. (15) we thus need to study the large X limit of I U (X, β ) and I W (X, β ). Consider the latter: From the behaviour of the integrand at large X one finds that only y X 1 contributes and the integral can be roughly approximated by  I W (X, β ) is calculated from the contribution of its saddle point y sp > 1 defined by the competition between y 2X and e −β y X /X , i.e. y sp = 2X β 1 X . From the Laplace method [55] at large X one gets, up to a numerical prefactor, a scaling that is well verified numerically for all temperatures. From (15) we arrive at R ∝ 1 √ X e −β /2 1 + (β − 1)e β and γ m = 1 X which yields γ m ∝ e β R 2 . This holds for β not scaled with X, smaller than any power of X, while a scaling such that β 1 is still described by this saddle point. Therefore, since in this limit β 1 (i.e. T → ∞) one has γ = m and R = 1 as the repulsive IPL dominates (see Eq. (10)), we expect that the relation approximates well the curve and should be a lower bound to all finite X curves. Since both γ and R go to zero for X → ∞, the non-trivial part of the finite-X curves gets shifted to higher temperatures, i.e., closer to R = 1, and washed out in the large-X limit.