Solid-liquid coexistence of the noble elements. I. Theory illustrated by the case of argon

The noble elements constitute the simplest group of atoms. At low temperatures or high pressures they freeze into the face-centered cubic (fcc) crystal structure (except helium). We perform molecular dynamics using the recently proposed simplified ab initio atomic (SAAP) potential [Deiters and Sadus, J. Chem. Phys. 150, 134504 (2019)] . This potential is parameterized using data from accurate ab initio quantum mechanical calculations by the coupled-cluster approach on the CCSD(T) level. We compute the fcc freezing lines for Argon and find a great agreement with the experimental values. At low pressures, this agreement is further enhanced by using many-body corrections. Hidden scale invariance of the potential energy function is validated by computing lines of constant excess entropy (configurational adiabats) and shows that mean square displacement and the static structure factor are invariant. These lines (isomorphs) can be generated from simulations at a single state-point by having knowledge of the pair potential. The isomorph theory for the solid-liquid transition is used to accurately predict the shape of the freezing line in the pressure-temperature plane, the shape in the density-temperature plane, the entropy of melting and the Lindemann parameters along the melting line. We finally predict that the body-centered cubic (bcc) crystal is stable at high pressures.


INTRODUCTION
This paper we investigating the solid-liquid coexistence of noble elements.The companion paper, Ref. [1], presented results for the argon (Ar) parametrization of the simplified ab initio atomic (SAAP) potential recently suggested by Deiters and Sadus [2].Here, we extend the investigation to include the noble elements Ne, Kr, and Xe.The parameters for the SAAP potential are determined from ab initio quantum-mechanical calculations using the coupled-cluster approach [3,4] on the CCSD(T) theoretical level [5,6].
We perform molecular dynamics as presented in paper I [1]: We consider monatomic systems of N particles with mass m confined in a volume V with periodic boundaries for the number density ρ = N/V .Let R = (r 1 , r 2 , r 3 , . . ., r N ) be the collective coordinate vector.The potential-energy function is a sum of pair potential contributions, U (R) = N i>j ε v(|r i −r j |/σ) where the SAAP pair potential is v(r) = a 0 exp(a 1 r)/r + a 2 exp(a 3 r) + a 4 1 + a 5 r 6 . ( For each of the noble elements the six a i parameters can be found in Reference [2]; they are determined by fitting to results of the above mentioned ab initio calculations on dimers [3].The SAAP pair potential is truncated and shifted at r c = 4 in units of σ. Figure 1(a) shows the pair potentials of Ne, Ar, Kr, Xe in units of σ and ε (Table I) and, for comparison, of the Lennard-Jones (LJ) potential truncated and shifted at r c = 6.The SAAP pair potentials are parameterized to have the same minimum in MD units as the LJ pair potential.Note that the Ne pair potential appears to be quite "hard" while the Xe pair potential is more "soft" at short distances.This difference is reflected in the shapes of the solid-liquid coexistence lines, as shown below.We use the RUMD software package [7] to study systems of N = 5120 particles in an elongated orthorhombic simulation cell where the box lengths in the y and z directions are equal, and the box length in the x direction is 2.5 times longer.We performed molecular dynamics simulations for 2 22 ≃ 4 × 10 6 steps after equilibration (also 2 22 steps) using a leap-frog time-step of 0.004σ m/ε.This result in a simulation time of roughly 1.7 × 10 4 σ m/ε, corresponding to 33 ns in argon units.The temperature T and pressure p is kept constant using the Langevin type dynamics suggested by Grønbech-Jensen et al. [8].

THE SOLID-LIQUID COEXISTENCE LINES
The solid-liquid coexistence lines of Ne, Kr, and Xe are computed as for Ar in Paper I [1]: The interface pinning method [9] is used to compute the coexistence state-point at the reference temperature T 0 = 2ε/k B .Table I shows the estimated coexistence pressures p, volume per particle in the liquid state V l /N , volume per particle in the solid-state V s /N , the volume change per particle at melting ∆V m /N , and the entropy of melting per particle ∆S m /N .Other coexistence points are subsequently determined by numerical integration of Clausius-Clapeyron identity, dP/dT = ∆S m /∆V m , using the standard fourth-order Runge-Kutta algorithm.The slope dP/dT is evaluated from thermodynamic information derived from N pT simulations.The numerical integration is carried out with temperature as the independent variable, starting from the reference temperature T 0 = 2ε/k B and moving down to temperature 0.65ε/k B and up to 4ε/k B .Figures 2(a)-(c) show the resulting coexistence lines for Ne, Kr and Xe, respectively.For reference, the gray lines on each panel show the Ar coexistence line.The dotted red lines are empirical coexistence lines [10,11].The SAAP potential systematically overestimates the coexistence temperature at a given temperature.This is likely due to missing many body interactions, as discussed for the case of Argon in Ref. [1].
Deiters and Sadus investigated the gas-liquid coexistence lines for the SAAP potentials of the noble elements [12].Here we compute the coexistence between the liquid and the face-centred cubic (fcc) solid.In all, this information allows us to compute the gas-liquid-fcc triple points (Table II).The triple point temperatures of the elements are very similar, ranging from 0.64251(35)ε/k B for Xe to 0.66054(5)ε/k B for Ne.This is not surprising, given the similar shapes of the pair potentials (Fig. 1).The pair-potential parameters were chosen to have the same minimum in reduced units as the LJ model.However, the triple point temperature of the LJ model is somewhat higher, 0.694ε/k B [13], which we interpret as an effect of the broader range of attraction of the LJ pair potential compared to that of the SAAP potentials.
As an aside, we investigate the validity of the Simon-Glatzel equation for the coexistence pressure [14], We first use T ref = T 0 and p ref = p 0 .Figure 3(a) show fit to the SAAP coexistence lines where the a and c parameters are determined by the least square method (see Table III).The accuracy of the fit is within a few MPa (Fig. 3(b) show residuals).The triple point temperature is often used when fitting empirical data.Table IV gives parameters using p ref = 0 and T ref as a third fitting parameter (in addition to a and c).The accuracy of the fit is compatible for the two approaches (see Figs. 3(b) and 3(c)).With the latter procedure, the reference temperature almost identical to the triple point temperature: T ref ≃ T tp (since the triple point pressure is nearly zero for the relevant pressure scale).Table IV compare SAAP parameters with parameters from empirical data.The agreement is in good.The a parameter and T ref of the SAAP fit is systematically lower than the parameters determined from empirical data.This is likely due to missing many body interactions of the SAAP potential, as discussed for the case of Ar in Ref. [1].In the remainder of the paper we do not use the Simon-Glatzel approximate.

HIDDEN SCALE INVARIANCE
In the companion paper [1] we show that the Ar parameterization of the SAAP potential has hidden scale invariance (for the investigated state-points).This fact was used to make an accurate prediction of the shape of the coexistence line as well as of the variation of several properties along the coexistence line.Below, we apply this theoretical framework to Ne, Kr and Xe.
Hidden scale invariance implies the existence of lines in the phase diagram along which structure, dynamics and some thermodynamic properties are invariant in given reduced units to a good approximation.These lines, referred to as "isomorphs", are configurational adiabats, i.e., lines of constant excess entropy.They can be computed by numerical integration in the logarithmic density-temperature plane of the "density-scaling exponent" γ ≡ ∂ ln T ∂ ln ρ Sex .The density-scaling exponent can be computed in an N V T simulation from the fluctuations in potential energy and virial: γ = ∆W ∆U / (∆U ) 2 [15].Isomorphic state points can then be found by numerical integration, using for instance the recently introduced fourth-order standard Runge-Kutta method [16].The initial state point for the integration is chosen as the coexistence points at the reference temperature T 0 = 2ε/k B (Table I).The dots on Fig. 6 show the isomorphs in the density-temperature plane of solid (red) and liquid (blue) states for the four noble elements under study.For all the elements, including Ar [1], the isomorphs follow the boundary of the coexistence region (solid lines) with minor deviations.The largest deviations are consistently found for the solid phase near the triple point.Figure 2 show the same information in the temperature-pressure plane.
Figure 4(a) shows the density-scaling exponents of the elements along the liquid isomorphs as a function of the temperature.At low temperatures, near the triple point, the exponents are 5.6±0.3.This value is close to that of the LJ potential [17].For the SAAP potential, γ decrease at higher temperatures and densities as it does for the LJ model; however, the γ variation is larger for the SAAP el-TABLE I. Thermodynamic data for estimated coexistence state points at the (reference) temperature T0 = 2ε/kB obtained by the interface-pinning method [9].ements and γ even goes below the LJ infinite-temperature limit of four.We conclude that the LJ potential is insufficient in describing the configurational adiabats of the noble elements.The γ's decrease with increasing atomic number.The value of γ can be estimated from the pair potential [1], and the decrease of γ with increasing atom number is directly related to the softness of the pair interaction: the softer pair interactions of Xe explain why its γ is lower than that of Ne. Figure 4(b) shows the density-scaling exponents γ of the elements along the solid isomorphs.The conclusions are the same for the liquids.
Figure 5(a) shows the Pearson correlation coefficients between the virial W and the potential energy U in the N V T ensemble defined by R ≡ ∆W ∆U / (∆W ) 2 (∆U ) 2 .The correlation coefficient is close to unity, R > 0.92, demonstrating that the potential energy function has hidden-scale invariance.Thus the structure, dynamics, and certain thermodynamics quantities are expected to be nearly invariant (in reduced units).This was demonstrated for Ar in Paper I [1].The W U -correlation is slightly weaker than for  IV. the LJ model (green dashed line on Fig. 5(a)).Figure 5(b) shows that R > 0.98 for the isomorphs of the solid phases.

ISOMORPH THEORY OF THE SOLID-LIQUID COEXISTENCE LINE
We have established that the potential-energy functions of the SAAP elements obey hidden scale invariance.This fact allows one to use the framework presented in Reference [18] to make theoretical predictions of the shape and property variation along the coexistence line.The basic idea is to make a first-order Taylor expansion from the isothermal state-points along the isomorphs, utilizing the fact that the coexistence lines are almost isomorphs.Property variations along an isomorph can be predicted from a single state point [1], and thus, information is in principle only needed at a single coexistence  state point.We use information along a liquid and a solid isomorph (computed by numerical integration of γ).The accuracy of the theory was demonstrated for Ar in Paper I [1].The green dashed lines in Figs.2(a)-(c) show the theoretical predictions for Ne, Kr, Xe, respectively.The agreement is excellent to the point that the lines are barely visible as they fall on top of the true coexistence lines (shown in solid black).A comparison with the results of the LJ mode (Fig. 2(d)) shows, however, that these results are slightly worse than for the LJ model.This is to be anticipated since the W U correlation coefficient R is weaker (Fig. 5), and the density scaling exponent γ is changing more rapidly (Fig. 4) for the SAAP potentials than for the LJ potential.
Figures 6(a)-(c) show the theoretical prediction of the coexistence region's boundaries in the densitytemperature plane as green dashed lines.The agreement is sound, but some deviations are noticeable at lower temperatures near the triple point temperature.For these temperatures, the density of the solid isomorphs are sev- eral percent lower than the density of melting.These deviations may come from the fact that only the firstorder terms in the Taylor expansion were included in the analysis.We hope to investigate this in the future.
The isomorph theory of the coexistence lines predicts both dynamical, structural and thermodynamic properties along the melting line.In Paper I [1] this was used for SAAP Ar to predict: i) the value of the diffusion constant along the liquid freezing line, ii) the Lindemann ratio of the solid along the melting line, and iii) the entropy of fusion ∆S m .As an example, Figs. 7(a)-(c) show the latter (∆S m ) for the remaining elements; Ne, Kr and Xe, respectively.The accuracy is comparable to that of Ar, but again slightly worse than that of the LJ model (Fig. 7(d)).For comparison, we note that hard-sphere based melting models predict ∆S m to be constant.Thus, the theoretical predictions of the isomorph framework are encouraging.

FIG. 2 .
FIG.2.The solid-liquid coexistence line in the pressure-temperature plane of (a) neon, (b) krypton, (c) xenon and (d) the LJ model.For reference, each figure shows the coexistence line of argon (gray line)[1].The dots represent solid and liquid isomorphic state points (generated from the reference state point at T0 = 2ε/kB), the green dashed line is the theoretical prediction of the isomorph theory (see below).Empirical melting lines are shown as red dotted lines[10,11].

FIG. 3 .
FIG. 3. (a) Fit of the Simon-Glatzel equation (Eq.2; dashed lines) to the SAAP solid-liquid coexistence line (dots) of Ne (green), Ar (red), Kr (purple) and Xe (blue).The reference points (T0, p0) are indicated with arrows.Parameters for Eq. 2 are given in Table III.(b) Residuals using T0 as the reference temperature (indicated with arrows).(c) Residuals using p ref = 0 and with the reference temperature T ref as a fitting parameter (T ref ≃ Ttp is indicated with arrows).Parameters for Eq. 2 are given in TableIV.

FIG. 5 .
FIG. 5. Person's correlation coefficient R between the virialW and the potential energy U .R is computed along the liquid isomorphs (a) and the solid isomorphs (b) and is shown as a function of temperature T .The fact that R is close to unity implies that the potential-energy function has hidden scale invariance along these lines[15].

FIG. 7 .
FIG. 7. The entropy of fusion ∆Sm of Ne (a), Kr (b), Xe (c) and LJ (d).The solid black lines are the true ∆Sm and the green dashed lines are the theoretical predictions of the isomorph theory.Hard-sphere based melting theories predict a constant entropy of fusion.

TABLE II .
Thermodynamic data for the estimated triple points.

TABLE III .
Parameters for the Simon-Glatzel equation (Eq.2) of SAAP coexistence state points using T ref = T0 and p ref = p0 as reference state-point.The fit is shown on and Fig. 3(a).

TABLE IV .
Parameters for the Simon-Glatzel equation (Eq. 2 ) of SAAP coexistence state points and of empirical data using the triple point temperature as reference: T ref = Ttp.For the SAAP results, T ref is treated as a fitting parameter and p ref = 0.