The EXP pair-potential system. IV. Isotherms, isochores, and isomorphs in the two crystalline phases

This paper studies numerically the solid phase of a system of particles interacting by the exponentially repulsive pair potential, which is a facecentered cubic (fcc) crystal at low densities and a body-centered cubic (bcc) crystal at higher densities [U. R. Pedersen et al., J. Chem. Phys. 150, 174501 (2019)]. Structure is studied via the pair-distribution function and dynamics via the velocity autocorrelation function and the phonon density of states. These quantities are evaluated along isotherms, isochores, and three isomorphs in both crystal phases. Isomorphs are traced out by integrating the density-temperature relation characterizing configurational adiabats, starting from state points in the middle of the fcc-bcc coexistence region. Good isomorph invariance of structure and dynamics is seen in both crystal phases, which is notable in view of the large density variations studied. This is consistent with the fact that the virial potential-energy correlation coefficient is close to unity in the entire fcc phase and in most of the bcc phase (basically below the re-entrant density). Our findings confirm that the isomorph theory, developed and primarily studied for liquids, applies equally well for solids. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5144871., s


I. INTRODUCTION
The EXP pair potential is the purely repulsive exponentially decaying function, v EXP (r) = ε e −r/σ . (1) The system of particles interacting via this pair potential has been studied much less than, e.g., the Lennard-Jones or inverse power-law (IPL) pair-potential systems. Nevertheless, papers reporting studies of the EXP system or, more generally, systems that involve an EXP term in the pair potential have appeared regularly ever since 1929. [1][2][3][4][5][6][7][8][9][10][11][12][13][14] As argued in Ref. 15, the EXP system deserves a closer study for two reasons. First, real-world systems, e.g., the low-density limit of the Yukawa (screened Coulomb) potential system, are well described by the EXP pair potential, which also plays an important role in most potentials describing metals. Second, the EXP pair potential may be regarded as the "mother of all pair potentials" in the sense that it explains the quasiuniversality of the structure and dynamics of simple liquids that applies for a large class of pair potentials. 16,17 Thus, any pair potential, which can be written as a finite sum of EXP terms with coefficients much larger than kBT, defines a system in the same quasiuniversality class as the hard-sphere and Lennard-Jones systems. 16,17 The EXP system may, consequently, replace the discontinuous hard-sphere system as the generic system in liquid-state theory. The first paper in this series studied the structure and dynamics along the EXP system's fluid-phase isotherms and isochores. 18 Paper I 18 also provided an example of quasiuniversality by showing that the radial distribution function (RDF) of the Lennard-Jones system is close to that of the EXP system at state points where the two systems have the same reduced diffusion coefficient. Paper II studied the EXP system's fluid phase isomorphs, 19 demonstrating shows the phase diagram in the densitytemperature plane. At low pressures/densities, the solid phase is a face-centered cubic (fcc) crystal, while at higher pressures/densities, the solid phase is body-centered cubic (bcc). The black dot marks the triple point at which the three phases meet. The star marks the re-entrant state point above which no solid phase exists. The discrete gray dashed lines in (b) are the three isomorphs studied in Sec. V.
invariance of the reduced-unit structure and dynamics along isomorphs. 19 Note that while there, as for any purely repulsive system, is just a single fluid phase, the EXP system has nevertheless typical gas and typical liquid state points. Moreover, in contrast to the Lennard-Jones system, the EXP system has strong virial potential-energy correlations at all low-temperature fluid-phase state points, including the system's typical gas-phase state points. The region delimiting the liquid and gas phases defines the so-called Frenkel line, the different definitions of which [20][21][22] are isomorph invariant.

FIG. 2.
Excess isochoric specific heat per particle in the solid phase in units of k B , denoted byc ex V . Most state points are green, showing that the specific heat is close to that of a perfect harmonic crystal (3k B per particle, implying c ex V = 1.5). The arrow marks the zero-temperature coexistence density found in Paper III: ρ = 1.747 64 ⋅ 10 −2 in the EXP unit system. The line separating fcc and bcc state points found from simulations is full red; its extrapolation to the zero-temperature coexistence density 15 is dashed.
Paper III 15 established the thermodynamic phase diagram of the EXP system, which involves two crystalline phases: a facecentered cubic (fcc) phase at low densities and a body-centered (bcc) at higher densities. There is a re-entrant liquid phase at the highest densities studied, which are of order unity in the "EXP unit system" defined by σ and ε in Eq. (1). This means that above a certain temperature there is no stable solid phase, a property the EXP system has in common with other pair-potential systems such as the Gaussian-core model, 23 which do not have a diverging energy at zero particle separation. Figure 1, which is reproduced from Paper III, 15 shows the thermodynamic phase diagram of the EXP system in (a) the pressure-temperature representation and (b) the densitytemperature representation. The phase transformations between the three phases are all of first order, i.e., associated with a density change and a latent heat. The coexistence regions in the densitytemperature plot are barely visible in Fig. 1(b), however, except as a slight thickening of the phase-boundary lines. In particular, the fccbcc transition is only weakly first order; the relative density change is merely 0.3% at the liquid-fcc-bcc triple point and even lower (0.01%) at zero temperature. 15 As an example of a quantity monitored in the two crystalline phases, Fig. 2 shows the excess isochoric specific heat per particle in units of kB, denoted byc ex V , in a diagram in which the colors reflect different values. The prevalent color is green in both the fcc and bcc phases, showing that the excess specific heat is close to the 1.5kB per particle expected for a perfect harmonic crystal according to the Dulong-Petit rule. 24

II. METHODS
All simulations were carried out using the double-precision version of the Roskilde University Molecular Dynamics (RUMD) code, which is optimized for graphics-processing-unit (GPU) computing. 25 The Newtonian molecular dynamics (MD) simulations were standard NVT simulations using the leap-frog integrator and the Nosé-Hoover thermostat. Brownian dynamics was employed to check the Nosé-Hoover NVT MD results. All findings relating to structure and dynamics reported below are based on Newtonian MD simulations at state points for which Brownian dynamics gives the same structure and specific heat. This is the case except at very low of Chemical Physics (a) Potential-energy distributions for three different simulation methods at a low-temperature state point. The distribution of potential energies using Brownian dynamics is given by blue symbols and that of an all-particle-move Monte Carlo algorithm by the red symbols. These two methods result in the same potential energy distribution, which as expected is Gaussian. In contrast, a bi-modal distribution is observed for Nosé-Hoover (NVT) MD simulations (black symbols). This is caused by the near-harmonic state of the system, resulting in extremely slow equilibration. NVT MD simulations are of little use at such low temperatures. (b) The friction coefficients λ ≡ 1/(mμ) used in the Brownian simulations, where m is the particle mass. The symbols mark simulated state points. (c) The Brownian dynamics cutoffs used. In order to arrive at cutoff-independent results, it was necessary to use a cutoff larger than 6σ at the highest densities and at certain low-temperature state points. (d) State points simulated by standard Newtonian Nosé-Hoover NVT MD dynamics. All results for structure and dynamics reported below were obtained by this dynamics.
temperatures. As in Paper III, 15 a cut-off at 6σ was used for all MDsimulated state points; a larger cut-off was occasionally used for the Brownian dynamics simulations (see below). Brownian dynamics simulations were carried out using an integrator implemented in RUMD as part of this work. The integrator is given 26 by (in which n is the discrete-time index, μ is the mobility, and R is the vector of all particle coordinates) Here, Δt is the time step and N(0, 1) for each time step is a vector of random numbers drawn from a zero-mean unit-variance Gaussian distribution. Note that changing the friction coefficient affects both the magnitude of the displacements and the relative importance of the random forces compared to the forces from the gradient. The with R depending primarily on the density, decreasing as density increases. Only at the highest densities studied, close to the re-entrant density, does the correlation coefficient go below the threshold of 0.9 that defines an R-simple system. 27,28 (b) The quantity 1 − R plotted as a function of the pressure along the solid-liquid phase boundary. The green curve represents the liquid phase at coexistence, and the red and blue curves represent the fcc and bcc phases at coexistence, respectively. Note that correlations are somewhat stronger in the solid phases than in the co-existing liquid phase. At the highest pressure, which is above the re-entrant state point, R becomes negative challenge is to find a friction coefficient for which the exploration of the phase space is reasonably fast. Figure 3(a) shows the equilibrium potential-energy distributions obtained by three different simulation methods: NVT MD (black points), a standard small-step Monte Carlo method implemented for an extra consistency check (red points), and Brownian dynamics (blue points). The latter two methods result in the same Gaussian distribution, while this is not the case for the MD simulation. The reason is that MD is highly inefficient at equilibrating an almost harmonic crystal, making it virtually impossible to use MD at very low temperatures. Depending on the quantity of interest, we therefore used Brownian dynamics at lowtemperature state points. Figures 3(c) and 3(d) show which state points were simulated by Brownian dynamics respectively MD, while Fig. 3(b) reports the friction coefficients used.
In Brownian dynamics, the configurational part of the isochoric specific heat calculated from the fluctuations should have the same value as when calculated from MD. This provides a consistency check ensuring that both methods reproduce the correct canonical ensemble probabilities. Based on this, we found it safe to use MD whenever the temperature is above 1.5 ⋅ 10 −5 or the density is less five times the melting density at the temperature in question. All structure and dynamics data reported below refer to this region of phase space. At selected state points here, we compared the structure, the density-scaling exponent (Sec. III), the virial potential-energy correlation coefficient (Sec. III), and cV ex with the same quantities obtained from Brownian Dynamics and found very good agreement.

FIG. 6.
The density-scaling exponent γ defined in Eq. (5). γ depends primarily on the density and decreases with increasing density. The same behavior is observed in the condensed-liquid phase of the EXP system. In contrast, at low densities and high temperatures where the EXP system is gas-like, γ depends primarily on the temperature [compare Fig. 3(c) in Ref. 19]. The simplest version of isomorph theory predicts that γ depends only on the density. 35 That treatment focused on the dense liquid phase in which each molecule interacts strongly and continuously with several nearest neighbors, which does not apply in the gas phase.

III. THE VIRIAL POTENTIAL-ENERGY CORRELATION COEFFICIENT AND THE DENSITY-SCALING EXPONENT
Isomorph theory applies for systems with strong correlations between the constant-volume equilibrium fluctuations of virial FIG. 7. Density-scaling exponent γ and virial potential-energy correlation coefficient R of the two crystalline phases in different representations. Red symbols and lines represent fcc state points, and blue symbols and lines represent bcc state points. (a) γ (left) and R (right) plotted vs density. The dashed and full curves are the predictions of the approximate analytical theory for the T → 0 limit of perfect crystals 28 (Appendix A). The predictions work best at low densities. (b) R as a function of temperature along two fcc isochores, showing that the approximate analytical theory (horizontal lines) works best at low density and low temperature. (c) γ plotted vs R for all simulated state points. The data almost collapse to a single curve common to both crystal structures. For comparison, the black dashed line gives the exact analytical prediction for the zero-density limit of the gas phase. 19 The red and blue lines are the fcc and bcc predictions of the approximate analytical theory 28 (Appendix A). The inset provides a blow up of the data with strongest correlations. (d) Log-log plot of γ vs 1 − R. The black dashed line is the gas-phase theory. 18 Green lines are data for liquid-phase isochores and purple lines for gas-phase isochores. In the left part of the figure, the identical red and blue lines give the theoretical predictions, while the triangles are numerical data for different isochores with the upper ones giving data for lower densities. The arrows give the direction of decreasing temperature on a given isochore. The isochores connect to liquid-phase isochores of virtually same density (green lines, data from Ref. 18). The dashed-dotted lines represent coexistence state points. The blue (bcc) data points are below the theoretical prediction while the fcc data are above it.

ARTICLE
scitation.org/journal/jcp W and potential energy U. [27][28][29][30][31] Recall that the average virial ⟨W⟩ gives the contribution to the pressure p deriving from particle interactions according to the general equation of state pV = NkBT + ⟨W⟩, in which N is the number of particles and T is the temperature.
Correlations between W and U are investigated numerically by standard Nosé-Hoover NVT MD simulations. For any configuration R defined as the 3N-dimensional vector of all particle positions, the potential energy U(R) is of course known from the simulation. For the EXP pair-potential system, the virial W(R) is given by the following sum over all pairs: 32 Note that the virial is always positive. This applies, in fact, throughout the phase diagram of any purely repulsive pair potential because the pressure is always larger than that of an ideal gas. Figure 4 shows plots of virial and potential energy equilibrium NVT fluctuations as a function of time with both quantities normalized to unit variance in order to make it easier to estimate the degree of correlation. 27 (a) and (b) show results for typical fcc and bcc crystalline state points, while (c) shows the fluctuations at a bcc state point close to that of the re-entrant state point. The latter case shows much weaker correlations.

ARTICLE scitation.org/journal/jcp
The Pearson correlation coefficient R quantifying the correlation between W and U is defined by (in which the angular brackets denote constant-volume canonical NVT averages and Δ the deviation from the equilibrium value) (3) Figure 5 gives data for the virial potential-energy correlation coefficient in which (a) shows by color coding how R varies. Interestingly, the fcc-bcc phase transition does not visibly change R. In most of the phase diagram investigated the EXP system obeys the R > 0.9 criterion designating an R-simple system. 17,27,28,33 Figure 5(b) shows how R varies along the liquid-solid phase boundaries, plotted as a function of the pressure. Below the horizontal dashed line, R > 0.9; this applies for most of the coexistence state points, involving both fcc and bcc phases. A more recent definition of an R-simple system is one for which most configurations obey the "hidden scale-invariance" condition 34 This expresses the property that if one configuration Ra has a lower potential energy than another one at the same density R b , this is also   (4) only applies for all configurations in the case of perfect virial potential-energy correlations (R = 1), 34 which only happens in the non-physical situation of a potential-energy function that is a constant plus an Euler-homogeneous function of the coordinates. Thus, in realistic situations, hidden scale invariance is merely an approximate symmetry. 17,33 Isomorphs are defined as curves in the phase diagram of constant excess entropy Sex for an R-simple system (Sex is the entropy minus that of an ideal gas at same density and temperature, 36 a quantity that is always negative). All systems have curves of constant Sex, obviously, but only for R-simple systems do these curves have approximately invariant reduced-unit structure and dynamics. 30,33,34,37 Hence, the name "isomorphs" for such curves. 30 The slope of an isomorph in the logarithmic densitytemperature phase diagram, γ, is defined 30 by For an R-simple system, if γ were constant, there would be invariance of the reduced-unit structure and dynamics along the lines of constant ρ γ /T. For this reason, γ is referred to as the densityscaling exponent. 30,33,38 In a computer simulation, one evaluates γ from the constant-density canonical ensemble equilibrium virial and potential-energy fluctuations via the following general statisticalmechanical identity: 30 Figure 6 shows how γ varies throughout the two solid phases. We see that γ depends primarily on the density. In fact, the original 2009 version of isomorph theory 30 predicted that γ is a function of the density independent of temperature. 35 It was later found, however, that the generic version of isomorph theory based on Eq. (4) allows for a more general variation of γ. 34 Interestingly, in the gas phase of the EXP system, γ depends more on temperature than on density. 19 It is appears that γ depends primarily on the density when a system is in the condensed liquid or solid phases in which any given particle interacts strongly with several nearest neighbors. 28,39 For many system, e.g., the Lennard-Jones and related systems, this corresponds to the main R-simple part of the thermodynamic phase diagram because (except at very high temperatures) the gas phase usually does not exhibit strong virial potentialenergy correlations. The EXP system is an interesting exception to this that exhibits strong correlations even in the gas phase (whenever kBT ≪ ε). 18 Figures 7(a) and 7(b) show γ and R as functions of density and temperature. In these figures, "predictions" refer to the theory of Ref. 28, which is summarized in Appendix A. This theory is not   18 derived analytically the black dashed line for the dilute gas phase and showed from simulations that all gas and liquid state points fall below this line. We now see that this also applies for the two solid phases. Interestingly, there is an almost one-to-one relation between R and γ for the solid state points, independent of the crystal phase. These data are consistent with an approximate theory (Appendix A) according to which a single analytical expression may be derived for R and γ, which applies for both crystalline phases. Figure 7(d) shows γ vs 1 − R in a log-log plot. The arrow indicates the direction of lower temperature for data of the same density. The theoretical prediction is approached in this limit for the fcc data (red) whereas the high-density bcc data fall below the prediction. The right part of Fig. 7(d) reproduces data for the fluid phase reported in Ref. 18.

IV. ISOTHERMS AND ISOCHORES
This section presents results for the variation of structure and dynamics along isotherms and isochores of the fcc and bcc phases of the EXP system; the same analysis is carried out along isomorphs in Sec. V. The structure is probed by the radial distribution function (RDF), while the dynamics is probed by the velocity autocorrelation function as well as its Fourier transform giving the phonon density of states. Figure 8 plots the RDF along four isotherms as a function of the reduced particle distancer ≡ ρ 1/3 r. (b) and (c) are for the same isotherm in the fcc and bcc phases, respectively; the same applies for (d) and (e). The structure changes significantly along all isotherms, with growing peak heights as density is decreased. Not surprisingly, the lowest temperature isotherm reported in (a) shows the largest and most narrow peaks, reflecting that here the thermal vibrations are smallest. Figure 9 shows data along the same four isotherms for the normalized velocity autocorrelation as a function of the reduced time defined byt ≡ ρ 1/3 √ kBT/mt. 30 The insets show the phonon density of states in reduced units, calculated from the velocity autocorrelation function. We conclude that the dynamics also changes significantly along isotherms.
We proceed to present a similar study of how structure and dynamics varies along the isochores, i.e., at constant volume. Figure 10 shows the reduced-unit RDFs along three isochores with (b) and (c) giving data for the fcc and bcc phases of the same isochore. At any given density, the lower the temperature is, the higher are the peaks. Figure 11 shows the normalized velocity autocorrelation function's variation with reduced time at the same three densities. Just as along the isotherms, there is a significant variation along any given isochore of both structure and dynamics.

V. ISOMORPHS
Section IV reported the variation of structure and dynamics along the EXP system's solid phase isotherms and isochores. Given the orders of magnitude of variation of density and temperature, not surprisingly both structure and dynamics vary a lot, even when given in reduced units in order to eliminate variations deriving from trivial scaling of time and length. As mentioned, the isomorph theory predicts the existence of lines in the phase diagram of approximately invariant reduced-unit structure and dynamics whenever virial potential-energy correlations are strong, which is the case for almost all state points of the two solid phases [compare Fig. 5(a)].
This section investigates to which degree this prediction holds. At a given temperature and pressure, the average of the fcc and bcc coexistence densities is chosen as the starting point for tracing out isomorphs in the fcc phase (going to lower densities) as well as in the bcc phase (going to higher densities). The isomorphic state points are determined by means of Eqs. (5) and (6) using a fourthorder Runge-Kutta integration in log space. 41 Four simulations are needed to determine a new state point in this approach. One of these four simulations refers to the given (starting) state point, the other three are used as correctors to the prediction. 41 The simulation at the state point in question was based on 2 × 10 6 time steps, while each of the three corrector simulations used just 50 000 time steps. The state point simulations were long because these were used also for calculating the RDF, velocity autocorrelation function, R, and γ-the Runge-Kutta integration scheme would work fine with only the short simulation time of the corrector simulations. Three isomorphs were generated by numerically integrating Eq. (5) starting from a state point on the fcc-bcc coexistence line as described above and moving, respectively, to lower and higher densities. Each direction involve one decade of density variation. The Runge-Kutta integration steps involved a density increase of ∼25% when going to higher density, i.e., into the bcc phase, and a density decrease of ∼20% when going to lower densities, i.e., into the fcc phase. The three isomorphs generated in this way and studied below are shown as gray dashed lines in Fig. 1(b), one of which is very close to the melting line. The precise state points are listed in Appendix B.
Before presenting data for the variation of structure and dynamics along selected isomorphs, we investigate the relation between melting line(s) and isomorphs. In the original paper on isomorphs from 2009, 30 it was predicted that the melting line is an isomorph. This statement was qualified in 2016, 42 and the melting line is now predicted to be an approximate isomorph, a fact that may nevertheless be utilized for calculating the melting-line variation of several quantities. It is generally a good approximation, however, to assume that the melting line is an isomorph. This is illustrated in Fig. 12 in which (a) shows a temperature-pressure plot of the liquid-phase isomorph through the triple point (black dashed line) and the two melting lines (fcc and bcc). The isomorph follows these lines closely, with significant deviations only at the highest densities. Here, the re-entrant point is approached, the correlation coefficient drops significantly (Fig. 5), and isomorph theory breaks down. Figure 12(b) gives the log-log density-temperature slopes of the fcc freezing, bcc freezing, and liquid melting lines, plotted as a function of the pressure. For comparison, Fig. 12(c) gives the log-log density-temperature slopes of the isomorphs (the density-scaling exponents) in the three phases emanating from the triple point. At the fcc-bcc transition, the slope is 2.33. This is consistent with expectations from the inverse-power-law pair potential where the

ARTICLE
scitation.org/journal/jcp transition takes place around the exponent n = 7 (corresponding to the slope γ = n/3). 43,44 Finally, Fig. 12(d) plots the density-scaling exponent γ in the solid phases along the melting lines as a function of the density. The overall conclusion of Fig. 12 is that both the fcc-liquid and the bcc-liquid melting lines are isomorphs to a good approximation and that, consequently, the melting-line slopes are well approximated by the density-scaling exponent. Significant deviations only occur close to the re-entrant point. Figure 13 demonstrates invariance of the reduced-unit RDFs along six isomorphs generated by Runge-Kutta integration of Eq. (5) into the fcc and bcc phases as described above. We observe a nice collapse along all three isomorphs. A comparison with Figs. 8 and 10 shows that the invariance of structure does not arise simply because structure is invariant throughout the phase diagram. Figure 14 shows approximate invariance of the dynamics along the same three isomorphs. The collapse is not as good as for the structure, but it should be kept in mind that a much larger density variation is considered here than in most previous studies. Thus the first data published demonstrating isomorph invariance of the dynamics were for a 3% density change of the highly viscous Kob-Andersen binary Lennard-Jones mixture, 30  comparison is with the isotherms of Fig. 9 or the isochores of Fig. 11, revealing a considerably better collapse along the isomorphs. Figure 15 shows how γ and R vary along the three isomorphs studied in Figs. 13 and 14. There is a substantial variation of γ, which as we have already seen is primarily a function of density. This shows that the EXP system cannot be approximated by an inverse powerlaw (IPL) pair-potential system. Since isomorphs are only exact for IPL systems, it is all the more remarkable that isomorph invariance of both structure and dynamics applies in the fcc and bcc phases of the EXP system. The reason is that, despite the large variation of γ, the virial potential-energy correlations are very strong (except at the largest densities studied).

VI. OUTLOOK
This paper has investigated the two crystalline phases of the EXP system. Except at the highest densities studied, close to the density of the re-entrant transition, the virial potential-energy correlations are very strong in both the fcc and the bcc phases (Fig. 5). Consistent with this, we find excellent invariance of structure along isomorphs and, in view of the large density range simulated, also a quite good invariance of the dynamics.
Only a few studies have been performed on isomorphs in crystals, 42,[45][46][47][48] but the findings of this paper confirm that hidden scale invariance and the concept of isomorphs is not limited to liquids. of Chemical Physics There is a much larger variation of the density-scaling exponent γ than found for Lennard-Jones type systems; at the same time R is close to unity at most state points.
Compared to other pair potentials, the EXP pair potential is interesting because it exhibits a large variation of the density-scaling exponent γ and has isomorphs in the gas, liquid, and solid phases.
In the latter case, the large γ variation results in two stable crystal structures. Based on the study presented in this paper, it is now possible to include the solid phase in the proposition that the EXP system may be regarded as the mother of all simple pair-potential systems. 17 Indeed, an exponential function often plays an important role in empirical as well as ab initio based pair potentials for metals. An interesting question for future work is what happens at densities above unity, i.e., above the re-entrant point. For the Gaussiancore model, unpublished results of ours show that R may become fairly close to −1, which gives rise to "anti-isomorphs" that have negative slope but still good invariance of the reduced-unit physics. Whether a similar thing happens for the EXP system is not known. If the answer is yes, this may be utilized for interpreting observations on real systems, e.g., dense solutions of star-shaped polymers as well as of metals. In relation to the latter, interestingly it has recently has been proposed that all metals may have a reentrant transition. 49

ACKNOWLEDGMENTS
This work was supported by the VILLUM Foundation's Matter Grant (No. 16515).
APPENDIX A: APPROXIMATE ANALYTICAL THEORY FOR γ AND R OF CRYSTALS Figure 7 presents predictions based on an approximate analytical theory developed some time ago. 28 The theory, which assumes that all interactions beyond nearest-neighbor pairs may be ignored and that transverse and longitudinal fluctuations of nearestneighbor vectors are uncorrelated, works as follows (compare Sec. II B.2 in Ref. 28). Let the index b represent nearest-neighbor bonds. If the nearest-neighbor bond b has length r b and the equilibrium bond length is ac, one defines At any time, the potential energy U and the virial W can both be expressed as linear combinations of S 1 , S 2 , S 3 , . . .. At low temperatures, S 1 /ac and S 2 /a 2 c have similar variance while the fluctuations of all the other terms are insignificant. 28 For the fluctuating parts of U and W, this leads to U = C U 1 S 1 /ac + C U 2 S 2 /a 2 c and 3W = C W 1 S 1 /ac + C W 2 S 2 /a 2 c . Here, C U p ≡ kpa p c /p! and C W p ≡ −kpa p c / (p − 1)! − k p+1 a p+1 c /p! in which kp is defined by Taylor expanding the pair potential in question as follows, v(r) = ∑ ∞ p=0 kp(r − ac) p /p! . For the EXP pair potential, if α ≡ ac/σ and β ≡ exp(−ac/σ), one finds kpa p c = ε(−α) p β. Thus, If the relative displacement of the pair of particles involved in bond b, u b , is decomposed into a term parallel to and a term perpendicular to the bond vector, u b = u b,∥ + u b, , it follows from Eq. (A1) that for small displacements one has S 1 = ∑b|ub, | 2 /(2ac) and S 2 = ∑b|ub,∥| 2 . Making the Einstein-crystal-like assumption that u b,∥ and u b, are uncorrelated and assuming that the length squared of the latter is twice that of the former because there are two transverse degrees of freedom, the theory predicts 28 that and that [in Ref. 28, the below factor of 1/3 is erroneously missing from the first line of Eq. (47)] From this, we find for the EXP pair-potential system and γ = 1 3 Nearest-neighbor interactions dominate in the low-density limit, i.e., when α ≫ 1; in this limit, the theory predicts R ≅ 1 and γ ≅ α/3. Equations (A5) and (A6) give the predictions shown as full curves in Fig. 7. For the fcc crystal, the α parameter is given by ρα 3 = √ 2, and for the bcc crystal, α is given by ρα 3 = 3 √ 3/4. According to the theory, R and γ are functions of α only. The theory thus predicts the same relation between these R and γ for both crystal phases. This is confirmed qualitatively in Figs. 7(c) and 7(d). As expected, the theory works best at low densities, though the theory is not exact even here.

APPENDIX B: STATE POINTS OF THE THREE ISOMORPHS STUDIED
The three isomorphs studied involve the state points listed in Tables I-III. The state points were determined by integrating Eqs. (5) and (6)