Density-scaling exponents and virial potential-energy correlation coefficients for the (2n,n) Lennard-Jones system

TL;DR

This study explores the relationship between the density-scaling exponent and virial potential-energy correlation coefficient in generalized Lennard-Jones systems across various dimensions and parameters, revealing a linear relation at low densities and a transition near high correlation values.

Contribution

It establishes a roughly linear relation between the density-scaling exponent and the correlation coefficient for generalized Lennard-Jones systems, extending understanding across multiple dimensions and parameters.

Findings

Linear relation $oxed{ ext{γ} oughly 3nR/d}$ at low densities.

Transition in $ ext{γ}$ versus $R$ plot around $R oughly 0.9$.

As $R ightarrow 1$, $ ext{γ} ightarrow 2n/d$, indicating dominance of repulsive interactions.

Abstract

This paper investigates the relation between the density-scaling exponent $\gamma$ and the virial potential-energy correlation coefficient $R$ at several thermodynamic state points in three dimensions for the generalized $(2n,n)$ Lennard-Jones (LJ) system for $n=4, 9, 12, 18$, as well as for the standard $n=6$ LJ system in two, three, and four dimensions. The state points studied include many low-density states at which the virial potential-energy correlations are not strong. For these state points we find the roughly linear relation $\gamma\cong 3nR/d$ in $d$ dimensions. This result is discussed in light of the approximate "extended inverse power law" description of generalized LJ potentials [N. P. Bailey et al., J. Chem. Phys. 129, 184508 (2008)]. In the plot of $\gamma$ versus $R$ there is in all cases a transition around $R\approx 0.9$, above which $\gamma$ starts to decrease as $R$ approaches unity. This is consistent with the fact that $\gamma\rightarrow 2n/d$ for $R\rightarrow 1$, a limit that is approached at high densities and/or temperatures at which the repulsive $r^{-2n}$ term dominates the physics.

Figures

Equations

$$R\left(\rho ,T\right)=\frac{⟨\Delta W\Delta U⟩}{\sqrt{⟨{\left(\Delta W\right)}^2⟩⟨{\left(\Delta U\right)}^2}⟩}.$$

$$U\left({\mathbf{R}}_{\mathrm{a}}\right)<U\left({\mathbf{R}}_{\mathrm{b}}\right)\Rightarrow U\left(\lambda {\mathbf{R}}_{\mathrm{a}}\right)<U\left(\lambda {\mathbf{R}}_{\mathrm{b}}\right).$$

$$U\left({\mathbf{R}}_{\mathrm{a}}\right)=U\left({\mathbf{R}}_{\mathrm{b}}\right)\Rightarrow U\left(\lambda {\mathbf{R}}_{\mathrm{a}}\right)=U\left(\lambda {\mathbf{R}}_{\mathrm{b}}\right).$$

$${\mathbf{R}}_{\mathrm{a}}\cdot \nabla U(\lambda {\mathbf{R}}_{\mathrm{a}})={\mathbf{R}}_{\mathrm{b}}\cdot \nabla U(\lambda {\mathbf{R}}_{\mathrm{b}}).$$

$$W({\mathbf{R}}_{\mathrm{a}})=W({\mathbf{R}}_{\mathrm{b}}).$$

$${\rho }_1^{1/d}{\mathbf{R}}_1={\rho }_2^{1/d}{\mathbf{R}}_2\equiv \tilde{\mathbf{R}}.$$

$$S_{\text{ex}}(\mathbf{R})\equiv S_{\text{ex}}\left(\rho ,U\left(\mathbf{R}\right)\right).$$

$$U\left(\mathbf{R}\right)=U\left(\rho ,S_{\text{ex}}\left(\mathbf{R}\right)\right)$$

$$S_{\text{ex}}(\mathbf{R})=S_{\text{ex}}(\tilde{\mathbf{R}}).$$

$$U\left(\mathbf{R}\right)=U(\rho ,S_{\text{ex}}(\tilde{\mathbf{R}})).$$

$$\frac{d^2\tilde{\mathbf{R}}}{d{\tilde{t}}^2}=\frac{\mathbf{F}\left(\mathbf{R}\right)}{{\rho }^{1/d}{k}_{B}T}\equiv \tilde{\mathbf{F}}(\mathbf{R}).$$

$$\mathbf{F}=-\nabla U=-{\left(\frac{\partial U}{\partial S_{\text{ex}}}\right)}_{\rho }\nabla S_{\text{ex}}(\tilde{\mathbf{R}}).$$

$$\mathbf{F}=-{\rho }^{1/d}T\tilde{\nabla }{S}_{\text{ex}}(\tilde{\mathbf{R}})$$

$$\tilde{\mathbf{F}}=-\frac{\tilde{\nabla }{S}_{\text{ex}}(\tilde{\mathbf{R}})}{k_B}.$$

$$\Delta W\left(t\right)\simeq \gamma \Delta U\left(t\right).$$

$$\gamma =\frac{⟨\Delta W\Delta U⟩}{⟨{\left(\Delta U\right)}^2⟩}.$$

$$\gamma ={\left(\frac{\partial \ln T}{\partial \ln \rho }\right)}_{S_{\text{ex}}}.$$

$${\gamma }_2=\sqrt{\frac{⟨(\Delta W{)}^2⟩}{⟨{\left(\Delta U\right)}^2⟩}}.$$

$$\gamma ={\gamma }_{2}R.$$

$$v_{m,n}^{\text{LJ}}\left(r\right)=\frac{\epsilon }{m-n}\left[n{\left(\frac{\sigma }{r}\right)}^m-m{\left(\frac{\sigma }{r}\right)}^n\right].$$

$$v_{2n,n}^{\text{LJ}}\left(r\right)=\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{2n}-2{\left(\frac{\sigma }{r}\right)}^n\right]\text{.}$$

$$\frac{h\left(\rho,S_{\text{ex}}\right)}{T}\,=\,{\rm Const.}\$$

$$h\left(\rho ,S_{\text{ex}}\right)=\left(\frac{d}{n}{\gamma }_0-1\right){\left(\frac{\rho }{{\rho }_0}\right)}^{2n/d}-\left(\frac{d}{n}{\gamma }_0-2\right){\left(\frac{\rho }{{\rho }_0}\right)}^{n/d}.$$

$$\gamma =\frac{d\ln h}{d\ln \rho }.$$

$$\rho >{\left(\frac{d{\gamma }_0-2n}{d{\gamma }_0-n}\right)}^{d/n}{\rho }_0.$$

$$\gamma /n\simeq R.$$

$$v(r)\cong A{r}^{-p}+B+Cr.$$

$$\gamma =\frac{p}{d}.$$

$$p=m+n.$$

$$\gamma \simeq n.$$

Full text

I Introduction

In the past decade a class of systems has been identified that is well described by the isomorph theory Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a, b); Thomas B. Schrøder and Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Jeppe C. Dyre (2009); Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009); Thomas B. Schrøder and Nicoletta Gnan and Ulf R. Pedersen and Nicholas Bailey and Jeppe C. Dyre (2011); Trond S. Ingebrigtsen and Thomas B. Schrøder and Jeppe C. Dyre (2012); Jeppe C. Dyre (2014, 2016), sometimes referred to as R (Roskilde) simple systems Malins et al. (2013); Abramson (2014); Fernandez and Lopez (2014); Flenner et al. (2014); Prasad and Chakravarty (2014); Schrøder, Thomas B. and Dyre, Jeppe C. (2014); Buchenau (2015); Harris and Kanakubo (2015); Heyes et al. (2015); Schirmacher et al. (2015). This class is believed to include most van der Waals bonded and metallic liquids and solids, as well as most weakly dipolar or ionic systems Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a, b); Thomas B. Schrøder and Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Jeppe C. Dyre (2009); Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009); Thomas B. Schrøder and Nicoletta Gnan and Ulf R. Pedersen and Nicholas Bailey and Jeppe C. Dyre (2011); Trond S. Ingebrigtsen and Thomas B. Schrøder and Jeppe C. Dyre (2012); Jeppe C. Dyre (2014, 2016). The class does not include systems with strong directional bonds like hydrogen or covalently bonded systems Jeppe C. Dyre (2016).

To determine whether the isomorph theory holds for a given system one calculates the strength of the virial potential-energy correlations at the state points in question in $NVT$ simulations, i.e., in the canonical ensemble Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a, b). Isomorphs are curves of invariant structure and dynamics in the thermodynamic phase diagram, see Sec. II. A system has isomorphs if and only if it has strong virial potential-energy correlations Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009). Virial potential-energy correlations are quantified by the state-point dependent Pearson correlation coefficient $R$ defined as

[TABLE]

In Eq. (1) $\rho$ is the (number) density ($\rho=N/V$ with $N$ particles in volume $V$), $T$ is the temperature, $W$ and $U$ are the virial and potential energy, the brackets denote $NVT$ averages, and $\Delta$ is the instantaneous deviation from equilibrium mean values.

The isomorph theory is exact only for systems with a potential that is a constant plus an Euler-homogeneous function, in which case the correlation coefficient $R$ in Eq. (1) is unity for all $\rho$ and $T$ Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009). An example is the inverse-power-law (IPL) pair-potential systems. For non-Euler-homogeneous potentials, which includes all realistic systems, the quantity $R$ is less than unity. In that case the applicability of the isomorph theory is restricted to a certain region of the phase diagram, usually the condensed-phase region encompassing the solid and “ordinary” liquid phase.

A system is defined to be R simple whenever $R>0.9$ Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a), because this ensures that isomorphs exist to a good approximation Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009). The $R=0.9$ threshold is somewhat arbitrary, however. The present paper investigates what happens when $R$ falls below this threshold. From the isomorph theory there is little help, because it generally breaks down when correlations are no longer strong, i.e., when the correlation coefficient goes significantly below $0.9$ Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008b); Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009); Jeppe C. Dyre (2014). In 2014, however, an interesting paper appeared by Prasad and Chakravarty Prasad and Chakravarty (2014) establishing that Rosenfeld’s excess entropy scaling as well as density-scaling Roland et al. (2005), which may both be derived from the isomorph theory Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009), need not break down even if the virial potential-energy correlations are weak. Reference Prasad and Chakravarty, 2014 studied the transition to simple liquid behavior in computer simulations of modified water models. We have taken inspiration from this work to systematically study the region of the phase diagram where virial potential-energy correlations are not strong for a class of models that – in contrast to water models – do have sizable regions of strong correlations in the thermodynamic phase diagram. This is done by studying Lennard-Jones type systems at lower densities than have previously been in focus.

This paper presents a study of the LJ($2n,n$) class of potentials defined as single-component systems with LJ-type radially symmetric pair potentials with a repulsive term proportional to $r^{-2n}$ and an attractive term proportional to $r^{-n}$ ($r$ being the interparticle distance). Results are presented for the cases $n=4,9,12,18$ in three dimensions, as well as for the standard LJ potential ($n=6$) in two, three, and four spatial dimensions (Sec. III). Before doing this the isomorph theory is briefly reviewed (Sec. II). Section IV summarizes the paper’s findings.

II Theoretical/Computational

II.1 Aspects of the isomorph theory

R simple systems – previously termed “strongly correlating” which, however, led to confusion with strongly correlated quantum systems – were defined in 2009 by reference to strong virial potential-energy correlations Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a). This section presents the alternative definition of the same class of systems given in 2014 Schrøder, Thomas B. and Dyre, Jeppe C. (2014) that refers to their “hidden scale invariance” Jeppe C. Dyre (2014). The original formulation of the isomorph theory Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009) is recovered via a first-order Taylor expansion Schrøder, Thomas B. and Dyre, Jeppe C. (2014).

Consider an $N$-particle system in $d$ spatial dimensions. The system is defined to be R simple if for any two configurations corresponding to the same density that obey $U({\bf R}_{\rm a})<U({\bf R}_{\rm b})$, this ordering is preserved after a uniform scaling of the two configurations Schrøder, Thomas B. and Dyre, Jeppe C. (2014) ($\mathbb{R}$ is the $d\,N$ dimensional vector describing a configuration of $N$ particles in $d$ spatial dimensions). Formally, this condition is written Schrøder, Thomas B. and Dyre, Jeppe C. (2014)

[TABLE]

In computer simulations periodic boundary conditions are applied, and it is understood that the box size scales with $\mathbb{R}$. Equation (2) need not be obeyed for all configurations; as long as it applies for most of the physically relevant configurations, the predictions of isomorph theory are obeyed to a good approximation. In Eq. (2) the scaling factor $\lambda$ can be any positive real number, and the equation holds for scaling ”both ways”. Because of this an R simple system obeys Schrøder, Thomas B. and Dyre, Jeppe C. (2014)

[TABLE]

Thus if two configurations have the same potential energy, their scaled potential energies are also identical.

R simple systems are characterized by strong correlations between virial and potential-energy fluctuations in the $NVT$ ensemble. This can be derived from Eq. (3) in the following way. For two configurations, ${\bf R}_{\rm a}$ and ${\bf R}_{\rm b}$, at the same density with the same potential energy, Eq. (3) states that $U(\lambda{\bf R}_{\rm a})=U(\lambda{\bf R}_{\rm b})$. Taking the derivative of this identity with respect to $\lambda$ one obtains

[TABLE]

Using the definition of virial $W(\bf{R})\equiv{\bf{R}}\cdot\nabla U(\bf{R})/{\it d}$ Hansen and McDonald (2013), one gets for $\lambda=1$

[TABLE]

Thus $U(\mathbb{R})$ determines $W(\mathbb{R})$, implying perfect correlation. This applies for systems that satisfy Eq. (2) or, equivalently Eq. (3), for all configurations; it holds for realistic R simple systems to a good approximation. As a consequence, for such systems the correlation coefficient $R$ of Eq. (1) is close to unity, but not exactly unity. As mentioned, the threshold defining R simple systems has usually been taken to be $R=0.9$ Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a).

Two configurations of density $\rho_{1}$ and $\rho_{2}$, which scale uniformly into one another, i.e., $\mathbb{R}_{2}=\lambda\mathbb{R}_{1}$ for some $\lambda$, are related through the following equality (where $d$ is the dimension):

[TABLE]

The last equality defines the “reduced” (dimensionless) configuration vector $\tilde{\mathbb{R}}$. Henceforth reduced quantities are marked by a tilde. The reduced version of a physical quantity is obtained by using the “macroscopic” unit system in which the length unit is $\rho^{-1/d}$, the energy unit is $k_{B}T$, and the time unit is $\rho^{-1/d}/\sqrt{mk_{B}T}$ ($m$ is the particle mass) Yaakov Rosenfeld (1999).

The entropy of any physical system can be written as an ideal-gas term plus an excess term $S_{\text{ex}}$. The ideal-gas term is the entropy of an ideal gas at the density and temperature in question, the excess term is the contribution to entropy from the interactions between the particles. We define the microscopic excess entropy function, $S_{\text{ex}}({\bf R})$, to be the thermodynamic excess entropy Schrøder, Thomas B. and Dyre, Jeppe C. (2014); L.D. Landau and E.M. Lifshitz (1958) of the potential-energy surface the configuration belongs to, i.e.,

[TABLE]

The function on the right-hand side, $S_{\text{ex}}\left(\rho,U\right)$, is the thermodynamic excess entropy of the state point with density $\rho$ and average potential energy $U(\mathbb{R})$. In other words, the microscopic excess entropy of a configuration $\mathbb{R}$ is defined as the excess entropy of a thermodynamic equilibrium system of same density with average energy precisely equal to $U(\boldsymbol{R})$.

Inverting Eq. (7) we get

[TABLE]

in which the function on the right-hand side, $U\left(\rho,S_{\text{ex}}\right)$, is the thermodynamic average potential energy of the state point with density $\rho$ and excess entropy $S_{\text{ex}}$.

Equations (7) and (8) apply for any system L.D. Landau and E.M. Lifshitz (1958), but they are particularly significant for R simple systems for which they imply that the configurational adiabats are curves of invariant structure and dynamics. To prove this, we first note that, as shown in Ref. Schrøder, Thomas B. and Dyre, Jeppe C., 2014, Eq. (2) implies that the excess-entropy function is invariant under uniform scaling, i.e., it only depends on the reduced coordinate vector $\tilde{\mathbb{R}}$:

[TABLE]

The relation for the potential-energy function Eq. (8) consequently becomes

[TABLE]

Isomorphs are defined as the configurational adiabats, i.e., curves along which the excess entropy is constant, in the region of phase diagram where the system is R simple Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009). To demonstrate invariance of structure and dynamics along the isomorphs we show that Newton’s second law in reduced units is invariant along an isomorph. In complete generality, this law is in reduced units (note that the particle mass is absorbed into the reduced time)

[TABLE]

For R simple systems Eq. (10) implies

[TABLE]

Using $\nabla=\rho^{1/d}\tilde{\nabla}$ and $\left(\partial U/\partial S_{\text{ex}}\right)_{\rho}=T$, the above expression becomes

[TABLE]

or

[TABLE]

This expression reveals that for R simple systems the reduced force vector $\tilde{\boldsymbol{F}}$ is a function of the reduced coordinate vector, implying that Eq. (11) is invariant along an isomorph. Thus the dynamics is isomorph invariant in reduced units, which implies that the reduced-unit structure is also isomorph invariant Schrøder, Thomas B. and Dyre, Jeppe C. (2014).

II.2 The density-scaling exponent

When strong correlations between virial and potential energy are present, a constant of proportionality between the instantaneous fluctuations of these two quantities’ deviation from their equilibrium values, $\gamma$, can be introduced via

[TABLE]

The correlation coefficient associated with this linear regression is that given in Eq. (1). The exact definition of $\gamma$ is the following Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009):

[TABLE]

By applying a standard fluctuation relation and the volume-temperature Maxwell relation one can show Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009) that

[TABLE]

The number $\gamma$ is termed the density-scaling exponent Roland et al. (2005); Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009). It determines the configurational adiabats, which as mentioned are isomorphs in the R simple region of the phase diagram. If variations of $\gamma$ are insignificant, these curves are via Eq. (17) given by $\rho^{\gamma}/T=$ Const. In the general case, Eq. (17) can be used to trace out isomorphs step-by-step by repeatedly changing density by typically an amount of order 1%, calculating the temperature change via Eq. (17), recalculating $\gamma$ from Eq. (16) at the new state point, etc.

In early publications Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a, b) we identified the constant of proportionality $\gamma$ of Eq. (15) by the following symmetric fluctuation expression, which was later renamed $\gamma_{2}$ to distinguish it from $\gamma$ Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009):

[TABLE]

Whenever the correlation coefficient $R$ is close to unity, one has $\gamma\cong\gamma_{2}$ since the following applies

[TABLE]

Now that the main ingredients of the isomorph theory have been introduced, we proceed to present the simulation results.

III Results and discussion

III.1 Generalized Lennard-Jones pair potentials in three dimensions

If $r$ is the interparticle distance, the generalized LJ pair potential is defined as follows

[TABLE]

Here $m$ and $n$ are positive integers, and in order to ensure thermodynamic stability it is assumed that $m>n$. The constants $\sigma$ and $\varepsilon$ define the potential’s length and energy scales, respectively, and the normalization used in Eq. (20) ensures that the minimum pair potential energy is $-\varepsilon$ which is obtained at $r=\sigma$. Note that the normalization is different from that usually employed for the standard 12-6 LJ pair potential parametrized as $4\varepsilon\left[(r/\sigma)^{-12}-(r/\sigma)^{-6}\right]$.

The aim of our study is to investigate whether any relation between the correlation coefficient and the density-scaling exponent can be determined. This involves simulating several state points for which the virial potential-energy correlations are not strong. With this goal in mind, a particular case of the generalized LJ potential was simulated, the case where $m=2n$:

[TABLE]

The following four instances of this pair potential, henceforth denoted by LJ(2n,n), were simulated in three dimensions: $n=4,n=9,n=12$ and $n=18$ (the $n=6$ case corresponding to the standard LJ potential is considered later). Figure 1 shows that these four pair potentials are quite different as functions of the pairwise distance $r$.

The simulations were performed in the $NVT$ ensemble, i.e., for a constant number of particles ($N=864$ for $n=4,9,12$ and $N=4000$ for $n=18$) at constant temperature $T$ and constant volume $V$. The time step was $0.001$ in LJ units and the simulations were performed with periodic boundary conditions and a standard shifted potential cut-off at $2.5\sigma$. Using an FCC crystal as starting configuration, each system was equilibrated for $2\cdot 10^{6}$ time steps before the collection of data began. After equilibration, the simulation ran for $50\cdot 10^{6}$ time steps during which data were collected. For each of the four systems, six densities were considered, $\rho=0.25,0.50,0.75,1.00,1.25,1.50$ (LJ units); the temperature was varied within the range $0.25<T<5.00$ (LJ units). At each state point the correlation coefficient and the density-scaling exponent were calculated from Eq. (1) and Eq. (16), respectively.

What to expect for the behavior of $R$ when density is lowered? To answer this question we refer to the isomorph theory, which works well for the LJ and related pair potentials in the ordinary liquid phase. According to isomorph theory, if a reference state point $(\rho_{0},T_{0})$ and another state point $(\rho,T)$ are on the same isomorph, the ratio between $T$ and $T_{0}$ defines what may be termed the isomorph shape function $h\left(\rho,S_{\text{ex}}\right)$ via $T/T_{0}=h\left(\rho,S_{\text{ex}}\right)$ L. Bøhling and T. S. Ingebrigtsen and A. Grzybowski and M. Paluch and J. C. Dyre and T. B. Schrøder (2012). Thus each isomorph is mapped out in the phase diagram by

[TABLE]

The original isomorph theory predicted the function $h\left(\rho,S_{\text{ex}}\right)$ to be independent of $S_{\text{ex}}$, whereas the more correct theory from 2014 predicts that $h$ may vary slightly from isomorph to isomorph Schrøder, Thomas B. and Dyre, Jeppe C. (2014). In both cases, however, the analytical form of the density dependence is the same. For the $\text{LJ}(2n,n)$ pair potential it can be shown that the isomorph shape function takes the following form in $d$ dimensions L. Bøhling and T. S. Ingebrigtsen and A. Grzybowski and M. Paluch and J. C. Dyre and T. B. Schrøder (2012); Trond S. Ingebrigtsen and Lasse Bøhling and Thomas B. Schrøder and Jeppe C. Dyre (2012); Lorenzo Costigliola (2016); Schrøder, Thomas B. and Dyre, Jeppe C. (2014)

[TABLE]

Here $(\rho_{0},T_{0})$ is the isomorph’s reference state point, at which $\gamma_{0}$ is the density-scaling exponent – the ${S}_{\rm ex}$ (isomorph) dependence of the function $h$ is contained in $\gamma_{0}$. Combining Eqs. (22) and (23) one can map out an isomorph from information obtained by computer simulations at the reference state point, i.e., with no need to use Eq. (17) in a tedious step-by-step process.

Note that Eq. (22) and the definition of the density-scaling exponent Eq. (17) implies Ingebrigtsen et al. (2012)

[TABLE]

Thus when the density-scaling exponent is known at a single reference state point, the theory via Eq. (23) and Eq. (24) predicts how $\gamma$ varies with density along the isomorph in question.

In Eq. (23) there are two terms, one positive and one negative. At high densities the positive term dominates. Upon lowering the density a point is reached where $h(\rho,S_{\text{ex}})$ changes sign. Below this density the theory implies negative isomorph temperatures, compare Eq. (22), which shows that the isomorph theory must break down here. Since the theory works well whenever there are strong virial potential-energy correlations, one concludes that the correlation coefficient must decrease upon lowering density along a configurational adiabat, i.e., for generalized LJ systems in any dimension the isomorph theory predicts its own breakdown at sufficiently low densities.

To summarize, moving along a given isomorph it is possible to define a density below which the isomorph theory does not hold. The condition that the shape function is non-negative sets the following limit for the isomorph theory to work:

[TABLE]

We have seen that $R$ decreases as density is lowered along a configurational adiabat (isomorph). Moreover, according to Eq. (16) the scaling exponent must approach zero if $R\rightarrow 0$. Since $R=0$ implies $\gamma=0$, one expects that $\gamma$ also decreases upon lowering the density. If density is increased, on the other hand, the positive term of Eq. (23) dominates and the density-scaling exponent eventually approaches $2n/d$ due to the dominance of the repulsive $r^{-2n}$ term L. Bøhling and T. S. Ingebrigtsen and A. Grzybowski and M. Paluch and J. C. Dyre and T. B. Schrøder (2012); Trond S. Ingebrigtsen and Lasse Bøhling and Thomas B. Schrøder and Jeppe C. Dyre (2012); Lorenzo Costigliola (2016). At the same time $R$ tends to unity. After considering these two limits, we turn to the simulation results.

Figure 2 shows the state points studied, marking in each case whether it is a solid or liquid state point, etc. Most state points are supercritical in which case we marked them by an F indicating “fluid”. Figure 3 plots $\gamma$ versus $R$ along the isochores for the $\text{LJ}(2n,n)$ potentials with $n=4,9,12,18$. Our findings are consistent with the theoretically predicted result that as $R\rightarrow 1$, $\gamma$ approaches the high-density limiting values 8/3, 6, 8, and 12, respectively, although we never come really close to these values that apply whenever the attractive term of the potential may be completely ignored. The density-scaling exponent changes behavior and starts to decrease with increasing $R$ in the region where $R$ is around $0.9$ which, as mentioned earlier, is the threshold usually used to define R simple systems. In the opposite limit, $R\rightarrow 0$, we find $\gamma\rightarrow 0$, as expected.

Our data show that the change of behavior in $\gamma$ versus $R$ takes place at the boundary of the region in which the system is R simple. A more extensive study is needed to determine whether this holds also for other R simple systems. If confirmed, the relation between the change in the behavior of the density-scaling exponent and the value of the correlation coefficient could be useful in practice, because $R$ is not accessible in experiment while $\gamma$ is Gundermann et al. . A change in the density or temperature dependence of $\gamma$ might therefore be used for identifying in which regions of the phase diagram a liquid is expected to obey isomorph theory predictions like isochronal superposition, excess entropy scaling, etc Jeppe C. Dyre (2014).

Figure 4 summarizes the data of Fig. 3. There is a roughly linear behavior between $\gamma$ and $R$, but it does not extend above $R\approx 0.9$ where $\gamma$ as mentioned starts to decrease towards the value $2n/3$ predicted at high density (where the $r^{-2n}$ repulsive term dominates and the virial potential-energy correlations become virtually perfect). For $R>0.9$ the system is R simple and the density variation of $\gamma$ is well described by Eqs. (23) and (24). For this reason we henceforth focus on the region in which $R<0.9$.

Figure 4(b) plots the data differently and establishes for most of the $R<0.9$ data the following approximate proportionality:

[TABLE]

How to understand this? In Ref. Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre, 2008b it was shown that an LJ-type pair potential may be approximated by the “extended inverse power law (eIPL)” pair potential with an effective state-point-dependent exponent $p$ that is not simply the exponent of the repulsive term of the LJ-type potential:

[TABLE]

This approximation usually works very well within the entire first coordination shell Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008b). At typical condensed-matter low and moderate pressure state points the term $B+Cr$ contributes little to the forces on a given particle, because the sum of these terms over all nearest neighbors tends to be almost constant. This reflects the fact that if the particle is moved slightly, some nearest-neighbor distances increase and some decrease, but their sum stays almost constant Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008b). This implies that the total force contribution from the $B+Cr$ term is small and the physics is dominated by the $r^{-p}$ IPL term.

For an IPL pair potential $\propto r^{-p}$ the density-scaling exponent is given in $d$ dimensions by Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008b)

[TABLE]

As mentioned, the effective exponent $p$ of Eq. (27) varies with state point. In large parts of the phase diagram of the LJ(m,n) pair potential the effective exponent $p$ is considerably larger than $m$ Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008b). This is because even within the repulsive part of the pair potential, i.e., at distances below the potential minimum, there is a sizable contribution from the attractive $r^{-n}$ term making the repulsions significantly steeper than expected from the $r^{-m}$ repulsive term alone. At densities and temperatures dominated by the pair potential minimum, i.e., at typical low or moderate pressure condensed-matter state points, it may be shown from a curvature argument referring to the potential-energy minimum that the effective exponent of the generalized LJ pair potential Eq. (20) is given Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008b) by

[TABLE]

At sufficiently high density and/or temperature the repulsive $r^{-m}$ term dominates the physics, however, and here one finds that $p\rightarrow m$. These conditions arise as $R\rightarrow 1$.

For the standard LJ pair potential at typical state points the density-scaling exponent is between $5$ and $6$, corresponding via Eq. (28) to $p$ between $15$ and $18$ Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a, b). Based on the above, we expect that the LJ(2n,n) pair potential in a large part of its phase diagram is equivalent to an eIPL pair potential with $p\simeq 2n+n=3n$. Via Eq. (28), this corresponds in three dimensions to

[TABLE]

In summary, $R\rightarrow 1$ at state points of high density and/or temperature where the $r^{-2n}$ term dominates the physics. This behavior is hinted at in Fig. 4(b) above $R\simeq 0.9$ where $\gamma/n$ starts to decrease and systematically deviates from the black dashed line. The limit $\gamma\simeq 2n/3$ was not reached in our simulations, however, which shows that only at very high density or temperature the attractive term of the generalized LJ pair potential may be ignored. A linear extrapolation of most of the data to $R=1$ gives the limit value suggested by the eIPL approximation in conjunction with Eq. (29). We take this as an indication that the eIPL approximation works well at the state points simulated, whether or not these are characterized by strong virial potential-energy correlations.

Having identified the limits $\gamma\rightarrow 0$ for $R\rightarrow 0$ and $\gamma\rightarrow n$ for $R\rightarrow 1$ (as long as $R<0.9$, i.e., for most data), the next question one may ask is: why do the data follow the approximate linear relation Eq. (26)? A possible explanation refers again to the eIPL approximation. The virial is calculated in $d$ dimensions from the pair potential as a sum over terms of the form $rv^{\prime}(r)/d$ Nicholas P. Bailey and Ulf R. Pedersen and Nicoletta Gnan and Thomas B. Schrøder and Jeppe C. Dyre (2008a). For the eIPL pair potential Eq. (27) this gives $Apr^{-p}/d$ plus an approximate constant (compare the above argument). Thus for fluctuations away from the equilibrium value, if the linear term of Eq. (27) is uncorrelated to the IPL term, the eIPL approximation suggests that one may have effectively

[TABLE]

with $p\cong 3n$ for $d=3$. This must be, admittedly, a very rough approximation when the $W$ $U$ correlations are not strong, but if one nevertheless makes it, Eq. (1) and Eq. (16) imply

[TABLE]

This is Eq. (26) that summarizes our findings for most state points in three dimensions. A concise way of expressing this is that the exponent $\gamma_{2}$ of Eq. (18) is roughly constant throughout the non-R-simple region of the phase diagram, compare Eq. (19).

III.2 The Lennard-Jones pair potential in two, three, and four dimensions

We proceed to report results for the standard LJ system in two, three, and four dimensions. Details on how the simulations were performed can be found in Refs. Lorenzo Costigliola, 2016; Lorenzo Costigliola, Thomas B. Schrøder and Jeppe C. Dyre, 2016. Figure 5(a) shows the density-scaling exponent versus the correlation coefficient along three configurational adiabats in the liquid phase and two in the crystalline phase for the two-dimensional LJ system; (b) shows a similar plot along two configurational adiabats in the liquid phase of the three-dimensional LJ system.

Figure 6 shows the density-scaling exponent versus the correlation coefficient along the configurational adiabats defined by the critical points for the LJ system in two, three, and four dimensions. Data for the critical points are taken from Refs. Smit, B. and Frenkel, D., 1991; Potoff, Jeffrey J. and Panagiotopoulos, Athanassios Z., 1998; Hloucha, M. and Sandler, S. I., 1999. For the configurational adiabats in the liquid phase the relation between $R$ crossing the threshold value $0.9$ and the start of the decreasing behavior of $\gamma$ is observed. For the crystalline phase, the correlation coefficient is above $0.9$ for all the state points studied and no crossover is observed.

In $d$ dimensions the eIPL-justified relation Eq. (32) implies via $p=2n+n=3n$ that

[TABLE]

This is tested in Fig. 7(a) for the LJ simulations in $2-4$ dimensions. We see that Eq. (33) overall works well.

Finally, Fig. 7(b) tests all data presented in this paper versus Eq. (33). The data in the region where the system is not R simple conforms roughly to Eq. (33), unless $R<0$ which are state points of liquid-gas coexistence.

IV Conclusions

For the state points where the system is not R simple ($R<0.9$) the generalized Lennard-Jones systems have a roughly linear relation Eq. (33) between the density-scaling exponent and the virial potential-energy correlation coefficient of (unless $R<0$ corresponding to liquid-gas coexistence state points). This finding has been rationalized in terms of the extended inverse-power-law approximation. Above $R\cong 0.9$ the approximate linear relation $\gamma\propto R$ breaks down; when $R$ approaches unity, the density-scaling exponent starts to decrease with increasing $R$. $\gamma$ is expected to eventually obey $\gamma\rightarrow 2n/d$ as $R\rightarrow 1$ when the repulsive $r^{-2n}$ term dominates completely. This happens only at very large densities and/or temperatures, however.

Our numerical data show that the threshold value $R\cong 0.9$ has a hitherto not recognized significance. While it was previously known that this threshold delimits the state points for which the isomorph theory works well Nicoletta Gnan and Thomas B. Schrøder and Ulf R. Pedersen and Nicholas P. Bailey and Jeppe C. Dyre (2009); Jeppe C. Dyre (2014), it is a new finding that roughly the same threshold marks a crossover from the $\gamma\propto R$ behavior found for $R<0.9$ to the one of decreasing $\gamma$ with increasing $R$ predicted by Eq. (23) and Eq. (24). It would be interesting to find out how general the finding of a threshold at $R\cong 0.9$ is. If it is general, this opens for the possibility – at least in principle since the required experiments are challenging – that determining the experimental variation of the density-scaling exponent throughout the phase diagram may be used for determining where the system in question is R simple and where it is not.

Acknowledgements.

We are indebted to Daniele Coslovich for helpful discussions regarding the interpretation of our numerical findings. This work was supported by the Danish National Research Foundation’s grant DNRF61 and by the Villum Foundation’s grant VKR-023455.