The EXP pair-potential system. II. Fluid phase isomorphs

This paper continues the investigation of the exponentially repulsive EXP pair-potential system of Paper I with a focus on isomorphs in the low-temperature gas and liquid phases. As expected from the EXP system's strong virial potential-energy correlations, the system's reduced-unit structure and dynamics are isomorph invariant to a good approximation. Three methods for generating isomorphs are compared: the small-step method that is exact in the limit of small density changes and two versions of the direct-isomorph-check method that allows for much larger density changes. Results from approximate methods are compared to those of the small-step method for each of three isomorphs generated by 230 one percent density increases, covering in all one decade of density variation. Both approximate methods work well.


I. INTRODUCTION
This paper and its companion (Paper I, [1]) present an investigation of the EXP pair-potential system consisting of identical particles interacting via the purely repulsive pair potential, v EXP (r) = ε e −r/σ . (1) Both papers investigate the region of the thermodynamic phase diagram where temperature is so low that the finite value v EXP (0) plays little role for the physics, i.e., where k B T ε. The focus is, moreover, on the low-density gas and liquid phases, i.e., where ρσ 3 1. While the densities considered are low in relation to the σ parameter of Eq.
(1), it should be emphasized that at low temperature the densities are not, in fact, low relative to the effective hard-sphere radius, but actually typical for usual studies of simple fluid gas-solid systems.
The EXP pair potential has not been studied much on its own right, in fact even less than other purely repulsive pair potentials like the family of inverse-power law pair potentials [2][3][4][5][6][7][8][9]. In most cases, the exponential function appears as a term in mathematically more involved potentials, for instance: 1) giving the repulsive part of the Born-Meyer pair potential from 1932 [10] or in embedded-atom models of metals [11,12], 2) multiplied by a Coulomb term to give the Yukawa (screened Coulomb) potential [13,14], or 3) giving the attractive long-ranged part in a model that rigorously obeys the van der Waals equation of state in one dimension [15].
To a good approximation the EXP system conforms to the following "hidden-scale invariance" condition for uniform scaling of same-density configurations R a and R b [16] (R is the vector of all particle coordinates and U (R) the system's potential-energy function): Here λ > 0 is a scaling parameter. Equation (2) expresses that if the potential energy of configuration R a is lower than that of configuration R b , this applies also after a uniform scaling of both configurations. For most systems, including the EXP system, Eq. (2) is approximate and not obvious from the mathematical expression for the potential energy, hence the term "hidden scale invariance" [17,18]. Paper I demonstrated one of the consequences of Eq. (2) for the EXP pair potential, namely strong virial potentialenergy correlations in the constant-volume thermal equilibrium fluctuations. The EXP pair-potential system has stronger such correlations than, e.g., the Lennard-Jones system [19,20]. In fact, the EXP system's virial potentialenergy Pearson correlation coefficient R is larger than 99% in a large part of its phase diagram (see Paper I and Refs. 20 and 21). A system is termed "R-simple" if it has better than 90% correlation.
The fact that the EXP pair-potential system obeys Eq. (2) implies that it has isomorphs, which are curves in the phase diagram along which structure and dynamics are invariant to a good approximation when given in reduced units [16,22]. Isomorphs are defined as curves of constant excess entropy, i.e., they are the system's configurational adiabats. While all systems have configurational adiabats, only R-simple systems, i.e., those obeying the hidden-scale-invariance condition Eq. (2) to a good approximation, have invariant physics along their configurational adiabats [16,[22][23][24][25].
A derivation of simple liquids' quasiuniversality based on the EXP pair potential was given in Refs. 21 and 25; an alternative proof utilizing constant-potential-energy (N V U ) dynamics [26] was presented in Paper I. Both proofs are based on the fact that under certain conditions a sum of two EXP pair potentials describes a system, which to a good approximation has the same physics as that of a single EXP pair-potential system. From this one easily shows that quasiuniversality applies for any system with a pair potential that to a good approximation may be written as a sum of EXP terms with numerically large prefactors in reduced units (see below) [21]. The EXP pair-potential system is thus central for explaining the physics of simple liquids, which justifies a closer investigation of the properties of the EXP system itself.
The isomorph theory is based on the use of reduced units [2,22,23,27], which are different from those usually applied for presenting simulation data involving the potential-energy function's characteristic energy and length. Instead reduced units utilize macroscopic parameters that vary with the thermodynamic state point. Consider a state point of temperature T and density ρ = N/V (N is the number of particles and V the system volume). If the average particle mass is m, reduced units make quantities dimensionless by scaling with the length l 0 = ρ −1/3 , the energy e 0 = k B T , and the time t 0 = ρ −1/3 m/k B T [22,24,25]. Reduced quantities are denoted by a tilde, for instancẽ R ≡ R/l 0 = ρ 1/3 R is the reduced configuration vector. We can now make precise the statement that the isomorphs of an R-simple system are lines of virtually identical physic; this refers to the system's reduced-unit structure and dynamics.
Besides demonstrating approximate isomorph invariance of the EXP system's structure and dynamics, the present paper discusses methods for generating the EXP system's isomorphs in computer simulations. Before doing this we demonstrate numerically in Sec. II examples of the system's hidden scale invariance, and Sec. III presents results for the density-scaling exponent's variation throughout the thermodynamic phase diagram (the isomorph slope in the log-log density-temperature phase diagram). The next sections report results from computer simulations along isomorphs traced out in different ways. Sec. IV presents the "small-step" method which, in the limit of infinitely small density changes, rigorously identifies configurational adiabats; Sec. V discusses two versions of the so-called direct isomorph check method, both allowing for much larger density changes. Finally, Sec. VI gives a brief discussion.

II. THE EXP SYSTEM'S HIDDEN SCALE INVARIANCE
The excess entropy S ex of a thermodynamic state point is defined as the entropy minus that of an ideal gas at the same density and temperature [28] (note that S ex < 0 since no system is more disordered than an ideal gas). S ex is the non-trivial part of a system's entropy, the contribution deriving from interactions. In general, excess thermodynamic quantities are calculated by leaving out the momentum degrees of freedom in the partition function [28]; they obey all standard thermodynamic relations like T = (∂U/∂S ex ) ρ , etc.
If one defines the microscopic excess entropy function S ex (R) as the thermodynamic equilibrium excess entropy of the state point with the density ρ of the configuration R and with average potential energy equal to U (R), it is straightforward to show that the hidden-scale-invariance condition Eq. (2) implies a scale-invariant entropy function, i.e., S ex (R) = S ex (λR) in which λ is a scaling parameter [16]. This means that S ex (R) depends only on the configuration's reduced coordinates,R = ρ 1/3 R. Inserting this into the identity U (R) = U (ρ, S ex (R)) that defines S ex (R) (the function U (ρ, S ex ) is the thermodynamic average potential energy as a function of density and excess entropy), one arrives [16] at This is the basic identity characterizing an R-simple system from which the isomorph invariance of structure and dynamics follow [16]. This does not imply, however, that all of the excess thermodynamics is isomorph invariant since there is also a density dependence in Eq. (3). For instance, the Helmholtz and Gibbs free energies are not isomorph invariant. A recent paper detailing how to use the isomorph theory for predicting how a number of quantities vary along the melting line gives an example of how the Gibbs free energy variation along a Lennard-Jones system isomorph is determined [29]. Physically, the hidden-scale-invariance condition Eq.
(2) states that if configurations at some density are ordered according to their potential energy, this ordering is maintained if the configurations are scaled uniformly to a different density. Note that Eq. (2) -and thus Eq. (3) and the entire isomorph theory -are approximate except for systems with an Euler-homogeneous potential-energy function (plus a constant).
A consequence of Eq. (3) is that an R-simple system has strong correlations between its constant-density thermalequilibrium fluctuations of virial and potential energy. This property was documented for the EXP system in Ref. [21] and in Paper I, and previously also for many other systems, including molecular and polymeric systems [19,20,30,31]. Recall that the microscopic virial W (R) is defined from the change of potential energy upon a uniform scaling of all particle coordinates, i.e., W (R) ≡ (∂U (R)/∂ ln ρ)R since a uniform scaling leavesR unchanged. Substituting Eq. (3) into this expression one finds In other words, W (R) = W (ρ, S ex (R)) in which the function W (ρ, S ex ) is the thermodynamic virial, i.e., the average of the microscopic virial at the state point with density ρ and excess entropy S ex . In conjunction with Eq. (3), the identity W (R) = W (ρ, S ex (R)) implies perfect correlation between the virial and the potential energy at fixed density in the sense that one of these two quantities uniquely determines the other. Note that this one-to-one relation between W and U is predicted to apply whether or not configurations are selected from an equilibrium simulation, it applies for instance also during aging [32]. We proceed to demonstrate numerically that the EXP pair-potential system obeys the hidden-scale-invariance condition Eq. (2) to a good approximation. While most quantities below are reported in reduced units, density and temperature are by definition constant in reduced units. This makes it impossible to specify a state point using reduced units for density and temperature so numerical values of the density are reported below in units of 1/σ 3 and numerical temperatures in units of ε/k B . This is referred to as the "EXP unit system" (Paper I). (2) for the EXP pair-potential system. R is the virial potential-energy Pearson correlation coefficient. (a) Potential energies of 20 statistically independent configurations from an equilibrium simulation of a system with 1000 particles at density ρ = 10 −3 and temperature T = 10 −4 (red dashed line), scaled uniformly to different densities and plotted as a function of density. (b) Same data with the average potential energy subtracted at each density, making it easier to investigate the implication of Eq. (2) that no crossings take place. This is seen to apply to a good approximation. (c) Same as in (b) for configurations selected from a simulation at ρ = 10 −3 and T = 10 −2 . There are here more crossings, meaning that Eq. (2) is less accurately obeyed, consistent with the lower R.  From the simulations we selected 20 statistically independent configurations (separated by 5 · 10 5 time steps). Each configuration was scaled uniformly to a different density ρ in the range 0.25 · 10 −3 < ρ < 1.75 · 10 −3 , i.e., a factor of seven density variation is involved. Figure 1(a) plots the potential energies of the scaled configurations as a function of their density. Note that no new simulations were performed to generate this figure, we merely scaled the 20 configurations uniformly and then evaluated their potential energies. Not surprisingly, for all configurations the potential energy increases strongly with increasing density. While they visually follow each other closely in Fig. 1(a), the figure does not allow for checking Eq. (2). To do this, Fig. 1(b) plots the same data by subtracting at each density the average potential energy of the 20 scaled configurations, making it possible to use a much smaller unit on the potential-energy axis. There are few crossings of the curves. This confirms that the EXP system obeys Eq. (2) to a good approximation, i.e., that it is R-simple in this region of the thermodynamic phase diagram. Figure 1(c) shows a plot like (b) at the same density but the higher temperature T = 10 −2 . Here crossings are more common, implying that Eq. (2) is less accurately obeyed. This is consistent with the finding of Paper I that the virial potential-energy correlation coefficient decreases as temperature increases. Moving in the opposite direction, Fig. 1(d) shows data for T = 10 −6 where the pattern is similar to that of (b), but with even fewer crossings. In summary, Eq. (2) works well at low temperatures, but breaks down gradually as higher temperatures are approached.
The absence of crossings upon a uniform scaling of all particle coordinates is not trivial. Hidden scale invariance is exact only for inverse power-law pair potentials, but it applies to a good approximation for many pair potential systems, e.g., Lennard-Jones type systems [16]. In contrast, there are many crossings, e.g., for the Lennard-Jones Gaussian pair potential for a density change of merely 20% [16].

III. THE DENSITY-SCALING EXPONENT
At any given state point one defines the so-called density-scaling exponent γ [22] by For an R-simple system, since isomorphs are curves of constant S ex [22], γ gives the slope of the isomorph through a given state point in the log-log density-temperature phase diagram. If γ were constant, according to the isomorph theory there would be invariance of the reduced-unit structure and dynamics along the phase diagram's lines of constant ρ γ /T . This is the origin of the name "density-scaling exponent". Before isomorph theory was developed density scaling was demonstrated experimentally for many glass-forming liquids [33]. By means of standard thermodynamic fluctuation theory γ may be calculated from the canonical-ensemble constantdensity (N V T ) equilibrium virial and potential-energy fluctuations [22], This identity makes it straightforward to determine γ in simulations. Figure 2 shows results for the density-scaling exponent with (a) giving γ's density variation along isotherms and (b) giving γ's temperature variation along isochores.
In the figure we mark three phases: gas, liquid, and solid, although the EXP system like any purely repulsive system has no liquid-gas phase transition, but merely a single fluid phase. The rationale for distinguishing between gas and liquid states -depending on whether or not the radial distribution function has a well-defined minimum above its nearest-neighbor peak, see Paper I -is the quasiuniversality of simple liquids' structure and dynamics that includes the EXP system (Sec. VII of Paper I). There is a large transition region between typical gas and typical liquid states; the least dense isomorph studied below is located in this region and referred to as the "gas-liquid" isomorph. The data in Fig. 2 reveal two regimes: At high densities and low temperatures (liquid and solid phases) γ is mainly density dependent, at low densities (gas phase) γ is mainly temperature dependent. It is possible to construct approximate analytical theories for the two limiting behaviors. Consider first the highdensity case corresponding to the liquid and solid phases, both of which are characterized by strong interactions between a given particle and its several nearest neighbors. To explain the strong virial potential-energy correlations of the Lennard-Jones (LJ) pair-potential liquid, Ref. 34 developed an approach based on the "extended inverse powerlaw" (eIPL) approximation that works as follows. At typical liquid or solid densities, within the first coordination shell the LJ pair potential is very well approximated [34] by the eIPL pair potential v eIPL (r) = Ar −n + B + Cr .
If one imagines a particle being displaced within its nearest-neighbor cage, some distances increase and some decrease, but the sum of all nearest-neighbor distances remains almost unchanged. This means that the B + Cr term is almost constant, so the eIPL pair potential may effectively be replaced by v eIPL (r) ∼ = Ar −n +D. When added over all particle pairs, this implies that the hidden-scale-invariance condition Eq. (2) applies to a good approximation. The above explains the strong virial potential-energy correlations found for the LJ and similar pair-potential systems [20,35,36]. It also suggests a means for approximately calculating the density-scaling exponent γ. Note first that for the IPL system v(r) = Ar −n +D; this follows from Eq. (6) and the definition of the pair virial w(r) ≡ (−1/3)ru (r) [34]. How to identify an effective IPL exponent n at a given state point? For the eIPL pair potential one has v eIPL (r) = −nAr −(n+1) + C, v eIPL (r) = n(n + 1)Ar −(n+2) , and v eIPL (r) = −n(n + 1)(n + 2)Ar −(n+3) . This implies n = −2 − rv eIPL (r)/v eIPL (r). Thus n is given by n = n 2 (r) if one for any pair potential v(r) defines the r-dependent effective IPL exponent n p (r) [34] by (v (p) (r) is the pth derivative of v(r))      Realistic values of n are arrived at by using for r a typical nearest-neighbor distance, of course. For the EXP unit system this means putting in which Λ ∼ = 1 is a numerical constant. Using this value of r for the EXP pair potential leads for p = 2, via Eq. (8), Eq. (9), and Eq. (10), to the following estimate of the density-scaling exponent in the EXP system's condensed liquid and solid phases: This is the black dashed line in Fig. 2(a) in which Λ = 1.075 was determined to get the best fit to data. At high densities the different isotherms collapse onto the line, a collapse that at low temperatures takes place earlier than at high temperatures. Next we study the dilute gas limit. Here the system is characterized by much longer typical distances between the particles than the interaction range σ, i.e., ρ 1. In this limit particle interactions predominantly take place via two-particles collisions. As shown in Paper I, it is possible to construct an analytical theory for the correlation coefficient R in the gas phase by assuming that particle collisions are random and uncorrelated. It is likewise possible to calculate γ analytically in the gas phase via Eq. (6) (compare Appendix I of Paper I). The relevant equation is in which v is the EXP pair potential treated as an independent variable, w = (−1/3)rv (r) = ln(1/v)v/3 is the pair virial and averages are taken over the non-normalizable v-probability distribution p(v) in which The integrals may be worked out analytically in terms of π, Euler's constant 0.577..., and the Riemann zeta function evaluated at 3 (Appendix I of Paper I). Numerically, the result is At low temperatures (β 1) the dominant term is γ ∼ = ln β/3. This corresponds to the p = 0 effective IPL exponent in Eq. (9) evaluated at the distance at which the pair potential equals k B T , the typical distance of nearest approach in a collision. Our numerical data indicate that Eq. (15) becomes exact at low densities, compare Fig. 2(b) that shows Eq. (15) as the dashed black line. Figure 3(a) summarizes our numerical findings for R and γ in a single phase diagram [21]. The color coding gives R, the line-segment slopes give γ. The line segments mark how the isomorphs run in the phase diagram, basically parallel to the melting line that is itself an approximate isomorph [21,22]. Figure 3(b) plots all values of (R, γ) with the dashed line marking the gas-phase analytical prediction. An important conclusion from this figure is that as γ → ∞ one has R → 1. Since large γ corresponds to an effectively very strongly repulsive pair potential on the k B T energy scale, this means that one expects R → 1 in the hard-sphere limit of the EXP system obtained by following an isomorph to zero temperature.

IV. STRUCTURE, DYNAMICS, AND SPECIFIC HEAT ALONG THREE ISOMORPHS
The "small-step" method for tracing out an isomorph in the phase diagram is based directly on Eq. (5) and Eq. (6). This section investigates predicted invariances along three isomorphs, which serve as reference "true" isomorphs in the next section dealing with faster, approximate methods for generating isomorphs.
It is straightforward to calculate in an N V T computer simulation the canonical averages on the right-hand side of Eq. (6). Typical values of γ for the EXP system are between 0.5 and 5 (Fig. 2). If for instance γ = 3, upon a 1% density increase the temperature is to be increased by 3% in order to keep S ex constant (Eq. (6)). After this change of density and temperature one recalculates the right-hand side of Eq. (6), and so on. In this way an isomorph is traced  FIG. 3. (a) shows the density-temperature phase diagram with line slopes given by the density-scaling exponent and color coding indicating the virial potential-energy correlation coefficient R. Diamonds mark gas-phase state points, circles mark liquid or solid state points. The line segments give the isomorph slopes, compare Eq. (5). The black dashed line is the approximate melting-line isomorph (covering the entire coexistence region). (b) Density-scaling exponent γ versus the virial potential-energy correlation coefficient R for all state points simulated. Red symbols are gas state points, blue symbols are liquid and solid state points. The dashed line is the prediction of the analytical gas-phase theory, which is obtained by combining Eq. (13) with R = A3/ √ A2A4 derived in Appendix I of Paper I. For a given value of R the gas phase has the highest γ, for a given value of γ the gas phase has the lowest R.
out with a method that is, in principle, exact in the limit of infinitely small density changes. The method is tedious since many steps are needed if 1% density changes are used to generate an isomorph covering a large density variation. The many steps required also mean that, in order to avoid accumulation of errors, each state-point simulation must be long enough to provide accurate data. Initial simulations for 5%, 2%, and 1% density changes from various state points showed that the latter two give virtually indistinguishable results. We concluded that a 1% density change is small enough to be reliable.
Three isomorphs were traced out for systems of 1000 particles using 1% density changes to cover one decade of density based on 230 simulations, each involving 10 million time steps (Fig. 4). One isomorph is located in the gas-liquid transition region, one is in the dilute liquid phase, and a third one is the liquid phase near the melting line. Note that in the log-log phase diagram the isomorphs are virtually parallel to the melting line. This is because in a simplified version of the isomorph theory [22] the melting line is an isomorph and isomorphs are given by an expression of the form h(ρ)/T = Const [37]. Different isomorphs correspond to different constants, so the isomorphs are parallel to one another in the log-log density-temperature phase diagram. A more accurate melting theory is now available [29], but the corrections to the older understanding of Ref. 22 are relatively small.
In the following, data are presented for each isomorph for the starting state point and ten more state points evenly spaced on the logarithmic density axis (the coordinates of the selected isomorph state points are given in the Appendix). We first investigate how the structure changes along an isomorph in order to see whether structure is invariant. Figure 5 shows the structure probed by the reduced radial distribution function (RDF) at the selected state points for each of the three isomorphs. The density change along each isomorph spans one decade, the temperatures span more than two decades. These variations are large compared to the first isomorph-theory simulations covering density variations of just a few percent [22]. Nevertheless, deviations from collapse are small; these are mainly observed around the first peak [22,30].
For each of the three isomorphs, Fig. 6 shows the reduced mean-square displacement (MSD) as a function of reduced time at the same eleven state points. Good collapse is observed. This may be compared to the lack of collapse of the  reduced MSD along the EXP system's isotherms and isochores (Paper I). Figure 7(a) demonstrates that the reduced diffusion constant (derived from the long-time MSD) is invariant along each of the three isomorphs. The diffusion (in reduced units) is faster the more gas-like the isomorph is. Figure 7(b) shows a contour plot of the reduced diffusion constant that again illustrates its isomorph invariance. Figure 8 gives the reduced excess isochoric specific heat per particlec ex V calculated from the fluctuations in potential energy in the N V T (canonical) ensemble via the Einstein expressionc ex V = (∆U ) 2 /N k 2 B T 2 . The more gas-like the structure is, the lower is the excess specific heat because interactions become infrequent. In the original (2009) version of isomorph theory [22] the excess specific heat was predicted to be an isomorph invariant. Figure 8 shows that this is not the case; in factc ex V ∝ ρ 1/3 in the gas phase. This confirms the need for the 2014 revision of the isomorph theory [16] in which it was shown that the originally predicted isomorph invariance ofc ex V results from a first-order approximation to a simpler and more correct theory, which starts from the hidden-scale-invariance condition Eq. (2).
Along the gas-phase isomorph the density variation of C V is determined as follows. For individual particle collisions -in particular at low temperature -the EXP pair potential can be approximated by an effective inverse power-law pair potential ∝ r −n with exponent given by n = −d ln v/d ln r = −r v EXP (r)/v EXP (r) evaluated at the distance r = r 0 of closest approach (this corresponds to the p = 0 case of Eq. (9)). The distance r 0 is estimated from k B T = v EXP (r 0 ). In the EXP unit system this leads to γ = − ln T /3 since the density-scaling exponent is given by γ = n/3 (Eq. (8)). Substituting γ = − ln T /3 into Eq. (5) leads to a simple first-order differential equation for ln T as a function of ln ρ. Its solution is ln T = A(S ex )ρ −1/3 in which A(S ex ) is an integration constant. This implies that (∂ ln T /∂S ex ) ρ = A (S ex )ρ −1/3 , which via the identity C V = (∂S ex /∂ ln T ) ρ implies that C V ∝ ρ 1/3 along the gas-phase isomorphs. The dashed line in Fig. 8 has slope 1/3; it fits well to the density variation of C V along the gas-liquid isomorph. Figure 9 shows that the density-scaling exponent γ decreases significantly with increasing density along each isomorph. Notably, the system has R > 0.9, i.e., is strongly correlating even with γ values as low as 1.5.  Fig. 4 probed via the reduced radial distribution function (RDF). (a) Gas-liquid isomorph; this isomorph is located where the system has little structure but this is approximately invariant. (b) Dilute-liquid isomorph; the structure here exhibits more structure and is still approximately invariant. (c) Dense-liquid isomorph; the system has here a typical liquid-like structure that is well maintained along the isomorph.

V. DIRECT-ISOMORPH-CHECK APPROXIMATE ISOMORPH
In this section we use the same three starting state points as above to trace out approximate isomorphs using two versions of the so-called direct isomorph check. This method, which allows for much larger density jumps than 1%, is justified as follows.
Suppose simulations at the state point (ρ 1 , T 1 ) generate a series of configurations. For a different density ρ 2 we wish to determine the temperature T 2 for which the state point (ρ 2 , T 2 ) is on the same isomorph as (ρ 1 , T 1 ). Each of the generated configurations R 1 is scaled uniformly to density ρ 2 via R 2 = (ρ 1 /ρ 2 ) 1/3 R 1 . We denote the potential energies of R 1 and R 2 by U 1 and U 2 , respectively. Because R 1 and R 2 have the same reduced coordinates and S ex (R) for any R-simple system depends only on the configuration's reduced coordinate [16], one has S ex (ρ 1 , U 1 ) = S ex (ρ 2 , U 2 ). Consequently, Eq. (3) implies Since (∂U/∂S ex ) ρ = T and ρ 1 and ρ 2 are both fixed, this implies for ratio of the potential-energy variations among different R 1 configurations, ∆U 1 , to the variation among the scaled configurations' potential energy, ∆U 2 , that In other words, the slope of a U (R 2 ) versus U (R 1 ) scatter plot is T 2 /T 1 , which allows for an easy way to determine T 2 . This is the direct isomorph check (DIC) [16,22]. For it to work properly it is important that the potential energies of the original and the scaled configurations are well correlated [22]; for the EXP system we find correlation coefficients above 99.8% when the density is doubled. We proceed to compare isomorphs generated from ten successive DIC jumps to what is termed "single state point DIC" (SSDIC)-generated isomorphs that start from the same density ρ 1 . The latter method ultimately increases density by one order of magnitude, whereas the DIC-generated isomorphs cover one decade of densities via ten jumps that are equally spaced on the logarithmic density axis, each increasing density by 25.9%. Figure 10 shows the relative temperature differences between approximate and "exact" isomorphs as a function of density with (a) giving the DIC method results and (b) the SSDIC method results. As the step size increases, the SSDIC method systematically overshoots the temperature. In this light it may seem surprising that good collapse is still maintained for the physics (Fig. 11 and Fig. 12). This is because the deviations in temperature from the small-step method temperatures obtained in Sec. IV are below 7%, even for temperature changes spanning two decades. Figure 13 shows the reduced diffusion constant and the reduced specific heats along the DIC-and SSDIC-generated approximate isomorphs, respectively, in both cases showing results close to the those of the exact isomorphs.

VI. DISCUSSION
As shown in Paper I, the EXP system has strong virial potential-energy correlations in the low-temperature part of its phase diagram. The present paper has demonstrated the existence of isomorphs in this part of the phase diagram. Three isomorphs were studied, one on the gas-liquid border, one in the dilute-liquid phase, and one in the condensedliquid phase. Each isomorph covers one decade of density; they were generated by 230 simulations using Eq. (6) in a step-by-step numerical integration of Eq. (5). Two approximate methods were studied for generating isomorphs numerically, the direct-isomorph-check (DIC) method and its single-state-point version (SSDIC). Both methods work well.
The virial potential-energy correlation coefficient R is very close to unity at the lowest temperatures studied T = 10 −6 , in fact here above 99.8% at all densities simulated (Paper I). From this one may be tempted to conclude that R → 1 for T → 0 at fixed density. We do not think this is the case, though. At any fixed density, the system  Not surprisingly, the isomorph with least interparticle interactions -the gas-liquid isomorph -has the lowest excess specific heat. A systematic increase with density is observed for all three isomorphs. This is at variance with the original (2009) version of isomorph theory in which CV is an isomorph invariant [22], but it is consistent with the 2014 version [16]. The dashed line has slope 1/3, which is the prediction for the density variation along gas-phase isomorphs (see the text).
eventually crystallizes upon cooling. Even though the crystal obeys R ∼ = 1 at low temperatures, there is no reason to expect R → 1 for T → 0 at a fixed density. The point is that the finite range of the EXP pair potential means that in a harmonic approximation there are not just nearest-neighbor springs, but also next-nearest neighbor springs, etc, and the different spring constants are not proportional to one another when density is changed to investigate whether or not Eq. (2) applies rigorously. We conjecture, however, that R → 1 for T → 0 along any isomorph. Here the physics is invariant, and as the temperature is lowered towards zero, the density-scaling exponent γ increases towards infinity because the density   shows results where each of the ten jumps were generated by the SSDIC method, i.e., jumping from the same low-density starting point. Not surprisingly, the latter method is less accurate than the DIC method, in particular for the largest density jumps, but in both cases deviations from the exact isomorph temperatures are small.
decreases. The EXP system becomes more and more hard-sphere (HS) like, and Fig. 3(b) shows that as γ diverges, R → 1 because R was empirically found always to be larger than its gas limit (for which one rigorously has R → 1 as ρ → 0, the dashed curve of Fig. 3(b)).
For future work it would be interesting to undertake a systematic investigation of the crystal phase. We have observed that the stable crystal structure is bcc at relatively high temperatures, but that it transforms into the fcc structure below certain temperatures (roughly when γ goes above two).

ACKNOWLEDGMENTS
We thank Lorenzo Costigliola for helpful discussions. This work was supported by the VILLUM Foundation's Matter grant (16515).