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). This paper investigates neon, argon, krypton, and xenon by molecular dynamics using the simplified atomic potentials recently proposed by Deiters and Sadus [J. Chem. Phys. **150**, 134504 (2019)], which are parameterized using data from accurate *ab initio* quantum-mechanical calculations by the coupled-cluster approach at the single-double-triple level. We compute the fcc freezing lines and find good agreement with the empirical values. At low pressures, predictions are improved by including many-body corrections. Hidden scale invariance of the potential-energy function is established by showing that mean-squared displacement and the static structure factor are invariant along the lines of constant excess entropy (isomorphs). The isomorph theory of melting [Pedersen *et al.*, Nat. Commun. **7**, 12386 (2016)] is used to predict from simulations at a single state point the freezing line’s shape, the entropy of melting, and the Lindemann parameter of the crystal at melting. Finally, our results suggest that the body-centered cubic crystal is the thermodynamically stable phase at high pressures.

## I. INTRODUCTION

The thermodynamic and transport properties of condensed matter systems are determined by their potential-energy functions.^{1} For a class of systems, the potential-energy function exhibits “hidden scale invariance,”^{2} making the phase diagram effectively one dimensional; thus, density and temperature collapse into a single parameter.^{2–13} This scaling, referred to as density scaling, has been demonstrated for simple model potentials,^{7,8,14–16} atomic liquids,^{17–22} and molecular liquids.^{4–6,9,23} In this paper, we investigate the noble elements neon (Ne), argon (Ar), krypton (Kr), and xenon (Xe) using a potential proposed recently from accurate *ab initio* calculations.^{24} We find that the energy surface obeys hidden scale invariance in the investigated part of the phase diagram and show how this fact can be used to predict the shape of the melting lines.^{25}

Specifically, we investigate below the simplified *ab initio* atomic potential (SAAP) recently suggested by Deiters and Sadus.^{24} This potential is parameterized for the noble elements Ne, Ar, Kr, and Xe from quantum-mechanical calculations using the coupled-cluster approach^{26,27} on the coupled-cluster single-double-triple [CCSD(T)] theoretical level.^{28} This approach is considered the gold standard of quantum chemistry^{29} and has been shown to give accurate prediction for the noble elements.^{27,30–33}

We consider monatomic systems of *N* particles of mass *m* confined to the volume *V* with periodic boundaries with number density *ρ* = *N*/*V*. Let **R** ≡ (**r**_{1}, **r**_{2}, **r**_{3}, …, **r**_{N}) be the collective coordinate vector. The potential-energy function is the sum of pair-potential contributions,

in which the dimensionless SAAP pair potential *ϕ*(*x*) $\u2261$$v$(*xσ*)/*ε* is given by

The six *a* parameters^{24} are determined by fitting to the results of the above-mentioned *ab initio* calculations on dimers.^{26} In the simulations reported below, the pair potential is truncated and shifted at *r*_{c} = 4*σ*. The truncation leads to tail corrections of the full potential.^{34} Pressure corrections are estimated to range from −0.3*ε*/*σ*^{3} to −0.1*ε*/*σ*^{3} depending on the state points (see the supplementary material).

An advantage of the SAAP is that it is computationally efficient and, at the same time, accurately represents the underlying *ab initio* calculations.^{24} Figure 1(a) shows the SAAP of Ar (in MD units defined by *ε* = *σ* = 1). The SAAP is compared to the Lennard-Jones (LJ) pair potential 4[*x*^{−12} − *x*^{−6}] (green dashed line; simulation details for the LJ system can be found in Ref. 25). As noted by Thiel and Alder,^{35} the LJ potential is very steep at short distances [Fig. 1(b)] but overall describes well the SAAP, though being a bit broader. In comparison, the exponential repulsive (EXP) pair potential 4 · 10^{5} exp(−12*x*)^{11,12,16,36,37} (red dashed lines) approximates well the SAAP at short distances. This is consistent with the interpretation of high-pressure compression experiments, the so-called shock Hugoniots.^{35,38}

Simulations were conducted using the RUMD software package v3.4^{39} (the SAAP was implemented in v3.4 by hand but is now included in v3.5). We studied systems of *N* = 5120 particles in an elongated orthorhombic simulation cell. The box lengths in the *y* and *z* directions are identical, while the box length in the *x* direction is 2.5 times longer to accommodate the below described *interface-pinning* simulations. More details of simulation algorithms are given in the supplementary material.

## II. THE COEXISTENCE LINE OF ARGON

The Ar coexistence line between the liquid and the fcc solid phases is determined as follows (the supplementary material includes details of the pair potential and the simulation algorithm): First, we use the interface-pinning method^{41–49} to compute the solid–liquid chemical potential difference Δ*μ* at temperature *T*_{0} = 2*ε*/*k*_{B} (287 K) for a range of fcc lattice constants corresponding to different pressures. Δ*μ* is computed from the thermodynamic force on a solid–liquid interface in an equilibrium simulation with an auxiliary potential that biases the system toward two-phase configurations. To account for slow fluctuations of the solid–liquid interface, two-phase simulations were eight times longer than the bulk simulations. We find that the coexistence pressure (Δ*μ* = 0) is *p*_{0} = 22.591(4)*ε*/*σ*^{3} at *T*_{0} in which the number in parentheses indicates the estimated statistical uncertainty on the least significant digit. Other coexistence points are then determined by numerical integration along temperatures on the coexistence line using the fourth-order Runge–Kutta (RK4) algorithm.^{50} The required slopes *dp*/*dT* are computed from isobaric simulations of a solid and a liquid using the Clausius–Claypeyron relation:^{51} if Δ*V*_{m} = *V*_{liquid} − *V*_{solid} is the volume difference between liquid and solid and Δ*S*_{m} = (*U*_{liquid} − *U*_{solid} + *p*Δ*V*_{m})/*T* is the entropy difference, the local slope is computed from^{51}

The RK4 algorithm requires four slope evaluations; thus, each integration step involves four simulations of a liquid and four simulations of a solid. This recipe for computing solid–liquid coexistence lines was suggested in Ref. 37. As a consistency check, we have confirmed that the gradient of the central difference of the computed melting line agrees with Δ*S*_{m}/Δ*V*_{m} [Eq. (3)]. This technique is in the class of *Gibbs–Duhem integration* methods suggested by Kofke.^{51} Details are provided in the supplementary material.

The coexistence line for Ar is shown as a black solid line in Fig. 2(a) together with the empirical (experimental) values^{40} shown with a red dashed line. The agreement is good; however, as seen in the inset, the SAAP slightly overestimates the coexistence temperature at a given pressure. This may be due to the missing many-body interactions of the SAAP. To investigate this, following the suggestion by Deiters and Sadus,^{52} we apply a mean-field correction that only depends on the average density of a bulk phase. For the correction of the coexistence pressure, we take the different densities of the phases into account. Let $\rho \xaf=2/(vs+vl)$ be the average density between the two phases (where $v$_{s} is the solid phase volume per particle) at a given coexistence point (*T*^{SAAP}, *p*^{SAAP}) computed with the SAAP. The corrected state point is then

in which

is an energetic correction parameter of the *ε* in Eq. (1), *ν*/*εσ*^{9} = 0.068 753 6 is the Axilrod–Teller–Muto parameter^{53} of Ref. 52, and *λ* is a parameter obtained from analyzing two-body and three-body simulation data.^{54} For simplicity, we set *λ* to unity. The blue solid line on Fig. 2 shows the corrected melting line. The correction is small, but it explains the deviations from the experimental melting line at low pressures [inset of Fig. 2(a)]. At high pressures, however, the uncorrected melting line is better than the corrected one. This suggests that many-body interactions are not important at high pressures, or that the above correction is inaccurate. For the remainder of the paper, we ignore many-body corrections, but note that inclusion of many-body effects will give some minor quantitative changes to our result. We also note that tail corrections to the truncated potentials will play a role (see the supplementary material). A study of these effects is left to future studies.

## III. HIDDEN SCALE INVARIANCE OF ARGON

The following gives a brief introduction to the theory of systems with hidden scale invariance, known as isomorph theory,^{2,7,8,10,12} and applies it to the Ar parametrization of the SAAP.

Consider two same-density configurations **R**_{a} and **R**_{b} where

If the energy surface is scale invariant, it follows that

where *λ* determines the magnitude of an affine scaling of all particle positions (and thus also the change of density). Scale invariance is trivial if *U* is a sum of inverse power-laws with the same exponent, so-called conformal potentials.^{4,5,7,14,15,55} However, the SAAP, LJ, and EXP potentials are not conformal. Nevertheless, depending on the state point in question, the above scaling can be a good approximation for relevant configurations. In that case, we refer to the scaling as *hidden* because it reflects an approximate property that is not obvious from the mathematical expression for the potential-energy function. From the definition of hidden scale invariance, it may be shown that there are lines in the thermodynamic phase diagram, referred to as isomorphs, along which structure, dynamics, and certain thermodynamics quantities are invariant in “reduced” units.^{2,8} These units are state-point dependent and defined from a combination of particle mass *m*, the number density *ρ*, and the thermal energy *k*_{B}*T*.^{2,8} Isomorphs are defined as lines where the excess entropy *S*_{ex} is constant, i.e., an isomorph is a configurational adiabat. Here, “ex” refers to the entropy *in excess* of the ideal gas entropy at the same density and temperature,^{1,3,22,56}

We note that the above equation in chemical thermodynamics is referred to as the *residual* entropy^{57,58} and should not be confused with the *excess entropy* defined in relation to ideal mixtures.^{3,59,60}

Scaling by reference to the excess entropy was first suggested by Rosenfeld^{3} and has recently gained renewed interest.^{13,20,56,57,61,62} As mentioned, a configurational adiabat is referred to as an *isomorph* for state points with hidden scale invariance, which implies invariant reduced structure and dynamics (*isomorph* is the greek word for same-shape). The slope of a configurational adiabat (and thus an isomorph) in the double logarithmic temperature-density plane,

can be computed from the fluctuations of virial *W* and potential energy *U* in the *NVT* (canonical) ensemble as^{8}

Here, ⟨⋯⟩ is the expectation value and Δ denotes the deviation from the mean (in practice, ⟨⋯⟩’s are computed as time-averages). The dashed lines in Fig. 2(b) show a liquid and a solid isomorph computed by numerical integration of Eq. (9) using the RK4 method from the reference state point (*T*_{0}, *ρ*_{0}) given in Table I.

. | ε/k_{B} (K)
. | σ (Å)
. | p_{0} (ε/σ^{3})
. | V_{l}/N (σ^{3})
. | V_{s}/N (σ^{3})
. | ΔV_{m}/N (σ^{3})
. | ΔS_{m}/N (k_{B})
. |
---|---|---|---|---|---|---|---|

Ne | 42.36 | 2.759 | 21.911(2) | 0.934 82(2) | 0.880 69(2) | 0.054 128(8) | 1.097 07(13) |

Ar | 143.5 | 3.355 | 22.591(4) | 0.926 12(3) | 0.874 89(3) | 0.051 227(11) | 1.085 15(13) |

Kr | 201.1 | 3.580 | 23.079(2) | 0.919 99(2) | 0.870 75(2) | 0.049 2320(9) | 1.073 43(16) |

Xe | 280.2 | 3.901 | 23.423(4) | 0.916 23(3) | 0.868 28(3) | 0.047 951(7) | 1.067 53(11) |

LJ | 20.8270(8) | 0.940 160(8) | 0.882 777(7) | 0.057 388(7) | 1.097 27(13) |

. | ε/k_{B} (K)
. | σ (Å)
. | p_{0} (ε/σ^{3})
. | V_{l}/N (σ^{3})
. | V_{s}/N (σ^{3})
. | ΔV_{m}/N (σ^{3})
. | ΔS_{m}/N (k_{B})
. |
---|---|---|---|---|---|---|---|

Ne | 42.36 | 2.759 | 21.911(2) | 0.934 82(2) | 0.880 69(2) | 0.054 128(8) | 1.097 07(13) |

Ar | 143.5 | 3.355 | 22.591(4) | 0.926 12(3) | 0.874 89(3) | 0.051 227(11) | 1.085 15(13) |

Kr | 201.1 | 3.580 | 23.079(2) | 0.919 99(2) | 0.870 75(2) | 0.049 2320(9) | 1.073 43(16) |

Xe | 280.2 | 3.901 | 23.423(4) | 0.916 23(3) | 0.868 28(3) | 0.047 951(7) | 1.067 53(11) |

LJ | 20.8270(8) | 0.940 160(8) | 0.882 777(7) | 0.057 388(7) | 1.097 27(13) |

Figure 3(a) shows that the liquid structure is indeed invariant along the isomorph. This is done by investigating the static structure factor

where

For a liquid, the structure factor depends only on the wave-vector length denoted by *q*. Figure 3(b) shows *S*(*q*) for state points along an isochore starting near the triple point. Figures 3(c) and 3(d) show *S*(*q*) for the fcc solid along state points of an isomorph and an isotherm, respectively. As for the liquid, the structure of the solid is isomorph invariant to a good approximation. We note that the long-wavelength (small *q*) limit of the structure factor does not scale well [inset of Fig. 3(a)]. This limit is proportional to the isothermal compressibility, which is not isomorph invariant even in reduced units.^{18}

Figure 4(a) shows that the dynamics is invariant along the liquid isomorph. This is done by investigating the mean-squared displacement and the diffusion coefficient (inset) computed from its long-time limit (dashed line). Figure 4(b) shows the same along an isochore.

The microscopic virial function is defined by

where $v$′ is the first derivative of the pair potential. The virial is referred to as the “potential part of the pressure” since the pressure *p* is given by the relation

where *W* = ⟨*W*(**R**)⟩. Systems with hidden scale invariance are sometimes referred to as “strongly correlating”^{7} because the fluctuations of virial and potential energy are strongly correlated in the *NVT* (canonical) ensemble. Figure 5 shows the Pearson correlation coefficient *R* between *W*(**R**) and *U*(**R**) of the canonical ensemble,

for the liquid (blue points) and the crystal isomorph (red points). The correlation is strong as expected from the invariant structure and dynamics (Figs. 3 and 4), and we find *R* > 0.94 for all investigated state points. The correlation increases with increasing temperature (and density). This is consistent with the fact that the structure is more invariant on the high-temperature part of the configurational adiabat (inset of Fig. 5).

If the pair potential is an inverse power-law (*r*^{−n}), the isomorphs are given by *ρ*^{γ}/*T* = const. in which *γ* = *n*/3.^{8,14} For this reason, *γ* is referred to as the “density-scaling exponent.” In general, *γ* is state-point dependent.^{7,8,63–65} It has been demonstrated that for many systems with pair interaction, including the LJ and the EXP systems, the exponent can be approximated by fitting an *effective* (state-point dependent) inverse power-law to the pair potential as some distance.^{66} For the dense phases (liquid and solid), this results in the expression

where $v$^{(i)}(*r*) is the *i*th derivative of the pair potential with respect to *r* and Λ(*S*_{ex}) is a free parameter for any given isomorph, expected to be close to unity. Under the assumption that Λ is the same for different isomorphs, *γ* is only a function of density to a good approximation. Figure 6 shows the true *γ* of Eq. (9) (dots) and the *γ* estimated from the pair interactions (solid lines) via Eq. (16), using the Λ that yields the correct *γ* at the reference temperature *T*_{0}. The agreement is excellent. Thus, an isomorph can be computed from a single reference state point by integrating Eq. (9) using Eq. (16) with Λ determined at the reference state point by calculating *γ* from the fluctuations via Eq. (10).In this paper, we investigate the fcc solid noting, however, that it is well known that noble elements can form other crystal structures such as hexagonal closed packing (hcp) at higher pressures than investigated here.^{68} At very high pressures (high densities and temperatures), the density-scaling exponent *γ* decreases and approaches that of the EXP potential. As such, the pair potential becomes long ranged compared to the interparticle distance. When *γ* < 2.3, it is expected that the body-centered cubic (bcc) crystal becomes thermodynamically stable compared to the closed packed crystals (fcc and hcp).^{15,17,69} Hoover and co-workers^{15} explained the fcc–bcc–fluid triple point for many metals by the fact that the effective pair potential becomes soft. It was recently shown that the EXP pair potential has a fcc–bcc–fluid triple point^{37} located where *γ* = 2.33(3).^{70} Since the noble-elements’ SAAP is approximated by the EXP potential at high densities, we expect that these elements have an fcc–bcc–fluid triple point where *γ* ≃ 2.3. For Ar, we estimate the fcc–bcc–fluid triple point to be at *T* = 21*ε*/*k*_{B} = 3000 K (see the inset of Fig. 6). Belonoshko and co-workers argued that this triple point exists for Xe, based on theoretical calculations and re-interpretations of experiments.^{71,72} The experiments presented in Ref. 73 were, however, unsuccessful in detecting such a triple point. We note that *γ* for the LJ potential is 4 in the high-pressure limit, and thus, an fcc–bcc–fluid triple point is not expected to exist.^{74} The LJ potential, however, does not describe the noble elements at high pressures since it is too harsh. Moreover, we note that the EXP high-pressure limit of the pair interactions also suggests a re-entrance temperature above which no crystalline phase is stable (see *T*^{⋆} in Ref. 70).

The density-scaling exponent *γ* can be determined from thermodynamic data as the ratio between the excess pressure coefficient $\beta Vex\u2261(\u2202W/\u2202T)V/V$ and the excess isochoric heat capacity per volume $cVex\u2261(\u2202U/\u2202T)V/V$,^{8}

These two thermodynamic quantities are usually not directly available from experiments. Using standard thermodynamic relations, *γ* can be rewritten, however, as

where

is the thermodynamic Grüneisen parameter,^{38,75–77}*c*_{V} is the isochoric heat capacity, *α* is the thermal expansion coefficient, and *K*_{T} is the isothermal bulk modulus. Within the Dulong–Petit approximation, *c*_{V} ≃ 3*k*_{B}, one finds (compare Ref. 78)

The value of the thermodynamic Grüneisen parameter is *γ*_{G} = 2.9^{13,79} near the gas–liquid–solid triple point of Ar. This corresponds to *γ* = 5.1, which is in good agreement with the value obtained by the SAAP (Fig. 6). Amoros *et al.*^{79} noticed that *γ*_{G} is only a function of density to a good approximation. This is explained by the fact that Λ in Eq. (16) is close to unity for all *S*_{ex}, making *γ* and *γ*_{G} functions solely of the density. Thus, *γ*_{G} is also only a function of density. This is only expected to be true in dense phases, i.e., for the liquid and the solid. In the gas limit, only temperature is expected to be relevant, as illustrated for the EXP potential in Ref. 16. The reason is that the typical collision distance of gas particles only depends on the temperature. In Fig. 6 (open circles), we compare the *γ* along the liquid isomorph of the SAAP potential to that of the empirical equation of state (EOS) obtained by Tegeler *et al.*^{67} (this EOS is implemented into the CoolProp software library by Bell *et al.*^{80}). The agreement is good.

To summarize, we have established that (i) the SAAP gives a good representation of Ar, (ii) the SAAP posseses hidden scale invariance near the coexistence line in both the liquid and solid phases, and (iii) the isomorphs can be computed from a reference coexistence point by means of Eq. (16) and Eq. (9).

Next, we apply the isomorph theory of the melting line^{25} to SAAP argon.

## IV. THEORY OF THE MELTING LINE

Using the LJ system as an example, Ref. 25 showed how the freezing and melting lines, as well as the variation of several properties along these lines, can be calculated from simulations carried out at a single coexistence state point by knowing the solid and liquid isomorphs through this state point. In particular, the coexistence pressure as a function of temperature, *p*_{m}(*T*), can be computed from the liquid and the solid isomorphs through a reference state point on the coexistence line with temperature *T*_{0}. This leads^{25} to

in which

and

Here, the superscript “I” indicates values along the isomorphs, the “0” indicates values at the reference state point, the subscript “s” indicates the solid value, and the subscript “l” indicates the liquid value. In Fig. 7(a), we test the prediction of the melting theory for Ar using the liquid and solid isomorphs defined by the reference temperature *T*_{0} = 2*ε*/*k*_{B} [shown in Fig. 7(b)].

Not only is the pressure variation along the melting line predicted but so are the variations of a number of other properties along the freezing and melting lines—again using only information obtained by simulations at the reference state point. The density of the liquid on the freezing line, *ρ*_{l}, is given by

where *W*_{m} = *p*_{m}/*ρ*_{l} − *k*_{B}*T* is the virial at melting. The density of the solid, *ρ*_{s}, has an analogous expression. The theoretical prediction is shown in Fig. 7(b) as dots. The prediction is good, though deviations are noted at low temperatures. Similar deviations (but small) are seen also for the LJ system.^{25} The deviation from the theory are not unexpected since R [Eq. (15)] is lower near the triple point.

Figure 8(a) shows the entropy of fusion Δ*S*_{m} per particle along the coexistence line (solid line). The theoretical prediction, shown as red dots, is quite good. For comparison, we note that the hard-sphere picture predicts a constant entropy of fusion. Figure 8(b) shows the Lindemann parameter of the fcc crystal,

where

is the root-mean-squared displacement of particles in the crystal at long times^{82} (⟨*u*^{2}⟩ = Δ*r*^{2}/3 in Ref. 82). The +’s in the inset show the Lindemann parameter along the *ρ* = 1.14 isochore. The dashed line is a linear fit that yields *∂δ*_{L}/*∂T*|_{ρ} = 0.0365*k*_{B}/*ε* needed in the theoretical prediction.^{25} The theoretical predictions are good for both Δ*S*_{m} and *δ*_{L}; however, they both become less accurate at low temperatures. This is related to the relatively poor prediction for the solid density near the triple point [Fig. 7(b)].

## V. COEXISTENCE LINES OF NEON, KRYPTON, AND XENON

The solid–liquid coexistence lines of Ne, Kr, and Xe are computed in the same way as for Ar. Figures 9(a)–9(c) show the results. For reference, the light gray lines in each panel show the Ar coexistence line [Fig. 7(a)] and Fig. 9(d) shows results for the LJ model. The dotted red lines are empirical coexistence lines.^{68,83} The SAAP potential systematically slightly overestimates the coexistence temperature at a given pressure. This is likely due to missing many-body interactions, as discussed for the case of Ar.

Deiters and Sadus investigated the gas–liquid coexistence lines for the SAAP potentials of the noble elements.^{52} Here, we compute the coexistence between the liquid and the face-centered cubic (fcc) solid. In all, this information allows one to compute the gas–liquid–fcc triple points (Table II). The triple-point temperatures of the elements are very similar, ranging from 0.642 51(35)*ε*/*k*_{B} for Xe to 0.660 54(5)*ε*/*k*_{B} for Ne. This is not surprising, given the similar shapes of the pair potentials (Fig. 1). However, the SAAP triple-point temperatures are 10%–12% above the empirical values, and the SAAP triple-point densities of the solids are 2%–5% higher than the empirical values. These systematic deviations are likely due to the missing many-body interactions, as discussed for the case of argon. The triple-point temperature of the LJ model is above that of the SAAP potential, 0.694*ε*/*k*_{B},^{84} which is possibly an effect of the broader range of attraction of the LJ pair potential compared to that of the SAAP (see Fig. 1).

. | . | T_{tp} (ε/k)
. | T_{tp} (K)
. | ρ_{s} (σ^{−3})
. | ρ_{s} (g/cm^{3})
. | ρ_{l} (σ^{−3})
. | ρ_{l} (g/cm^{3})
. | ρ_{g} (σ^{−3})
. |
---|---|---|---|---|---|---|---|---|

Ne | SAAP empirical | 0.650 54(5) | 27.558(2) | 0.952 39(5) | 1.519 37(8) | 0.825 92(16) | 1.317 6(3) | 0.004 21(25) |

0.579 77^{a} | 24.56 | 0.905 14 | 1.444 | 0.782 28 | 1.248 | |||

Ar | SAAP empirical | 0.646 79(34) | 92.80(5) | 0.949 53(5) | 1.667 72(9) | 0.825 32(35) | 1.449 6(6) | 0.004 25(25) |

0.584 04 | 83.81 | 0.924 07 | 1.623 | 0.802 80 | 1.410 | |||

Kr | SAAP empirical | 0.645 77(54) | 129.85(9) | 0.947 50(2) | 2.873 53(6) | 0.826 24(26) | 2.505 8(8) | 0.003 51(28) |

0.575 73 | 115.78 | 0.931 83 | 2.826 | 0.808 18 | 2.451 | |||

Xe | SAAP empirical | 0.642 51(35) | 180.02(10) | 0.945 34(4) | 3.471 25(15) | 0.823 78(15) | 3.024 88(55) | 0.004 38(11) |

0.575 94 | 161.37 | 0.925 66 | 3.399 | 0.837 70 | 3.076 | |||

LJ | 0.694^{b} | 0.96 | 0.84 |

. | . | T_{tp} (ε/k)
. | T_{tp} (K)
. | ρ_{s} (σ^{−3})
. | ρ_{s} (g/cm^{3})
. | ρ_{l} (σ^{−3})
. | ρ_{l} (g/cm^{3})
. | ρ_{g} (σ^{−3})
. |
---|---|---|---|---|---|---|---|---|

Ne | SAAP empirical | 0.650 54(5) | 27.558(2) | 0.952 39(5) | 1.519 37(8) | 0.825 92(16) | 1.317 6(3) | 0.004 21(25) |

0.579 77^{a} | 24.56 | 0.905 14 | 1.444 | 0.782 28 | 1.248 | |||

Ar | SAAP empirical | 0.646 79(34) | 92.80(5) | 0.949 53(5) | 1.667 72(9) | 0.825 32(35) | 1.449 6(6) | 0.004 25(25) |

0.584 04 | 83.81 | 0.924 07 | 1.623 | 0.802 80 | 1.410 | |||

Kr | SAAP empirical | 0.645 77(54) | 129.85(9) | 0.947 50(2) | 2.873 53(6) | 0.826 24(26) | 2.505 8(8) | 0.003 51(28) |

0.575 73 | 115.78 | 0.931 83 | 2.826 | 0.808 18 | 2.451 | |||

Xe | SAAP empirical | 0.642 51(35) | 180.02(10) | 0.945 34(4) | 3.471 25(15) | 0.823 78(15) | 3.024 88(55) | 0.004 38(11) |

0.575 94 | 161.37 | 0.925 66 | 3.399 | 0.837 70 | 3.076 | |||

LJ | 0.694^{b} | 0.96 | 0.84 |

As an aside, we investigate the validity of the Simon–Glatzel equation for the coexistence pressure,^{86}

We first use *T*_{ref} = *T*_{0} and *p*_{ref} = *p*_{0}. Figure 10(a) shows fit to the SAAP coexistence lines where the *a* and *c* parameters are determined by least-squares fitting (see Table III). The accuracy of the fit is within a few MPa [Fig. 10(b) shows the 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 [Figs. 10(b) and 10(c)]. With the latter procedure, the reference temperature is almost identical to the triple-point temperature: *T*_{ref}$\u2245$*T*_{tp} (since the triple point pressure is nearly zero for the relevant pressure scale). Table IV compares the SAAP parameters with parameters from empirical data. The agreement is in general good. The *a* parameter and *T*_{ref} of the SAAP fit are systematically lower than the parameters determined from the empirical data. This is likely due to the missing many-body interactions of the SAAP, as discussed for the case of Ar.

. | T_{0} (K)
. | p_{0} (GPa)
. | a (GPa)
. | c . |
---|---|---|---|---|

Ne | 84.72 | 0.610 10 | 0.751 | 1.4983 |

Ar | 286.98 | 1.184 98 | 1.436 | 1.5473 |

Kr | 402.16 | 1.396 45 | 1.681 | 1.5699 |

Xe | 560.37 | 1.526 07 | 1.824 | 1.5924 |

. | T_{0} (K)
. | p_{0} (GPa)
. | a (GPa)
. | c . |
---|---|---|---|---|

Ne | 84.72 | 0.610 10 | 0.751 | 1.4983 |

Ar | 286.98 | 1.184 98 | 1.436 | 1.5473 |

Kr | 402.16 | 1.396 45 | 1.681 | 1.5699 |

Xe | 560.37 | 1.526 07 | 1.824 | 1.5924 |

. | . | T_{ref} (K)
. | a (GPa)
. | c
. |
---|---|---|---|---|

Ne | SAAP | 27.75 | 0.1409 | 1.4989 |

empirical, Ref. 83 | 24.55 | 0.1286 | 1.4587 | |

Ar | SAAP | 92.91 | 0.2501 | 1.5487 |

empirical, Ref. 68 | 83.81 | 0.2245 | 1.5354 | |

Kr | SAAP | 129.64 | 0.2835 | 1.5713 |

empirical, Ref. 68 | 115.77 | 0.2666 | 1.4951 | |

Xe | SAAP | 179.45 | 0.2966 | 1.5942 |

empirical, Ref. 68 | 161.40 | 0.2594 | 1.4905 |

. | . | T_{ref} (K)
. | a (GPa)
. | c
. |
---|---|---|---|---|

Ne | SAAP | 27.75 | 0.1409 | 1.4989 |

empirical, Ref. 83 | 24.55 | 0.1286 | 1.4587 | |

Ar | SAAP | 92.91 | 0.2501 | 1.5487 |

empirical, Ref. 68 | 83.81 | 0.2245 | 1.5354 | |

Kr | SAAP | 129.64 | 0.2835 | 1.5713 |

empirical, Ref. 68 | 115.77 | 0.2666 | 1.4951 | |

Xe | SAAP | 179.45 | 0.2966 | 1.5942 |

empirical, Ref. 68 | 161.40 | 0.2594 | 1.4905 |

## VI. HIDDEN SCALE INVARIANCE OF NEON, KRYPTON, AND XENON

Figure 11(a) shows the density-scaling exponents of the elements along the liquid isomorphs plotted as functions 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.^{7} For SAAP, *γ* decreases at higher temperatures and densities as for the LJ model; however, the *γ* variation is larger for the SAAP elements, and *γ* eventually goes below the LJ infinite-temperature limit of 4. Thus, 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, and the decrease of *γ* with increasing atomic 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 11(b) shows the density-scaling exponents *γ* of the elements along the solid isomorphs. The conclusions are the same as for the liquids.

Figure 12(a) shows the Pearson correlation coefficients *R* between the virial *W* and the potential energy *U* [Eq. (15)]. 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 approximately isomorph invariant (in reduced units). This was demonstrated for Ar above, and we find similar results for the three other elements (results not shown). The *WU*-correlation coefficient *R* for SAAP is somewhat smaller than that of the LJ potential [green dashed line in Fig. 12(a)]. Figure 12(b) shows that *R* > 0.98 for the isomorphs of the solid phases.

## VII. ISOMORPH THEORY OF THE SOLID–LIQUID COEXISTENCE LINE FOR NEON, KRYPTON, AND XENON

Figures 13(a)–13(c) show the theoretical prediction of the coexistence region’s boundaries in the density–temperature plane as green dashed lines. The agreement is quite good, though with some deviations at lower temperatures near the triple point temperature. For these temperatures, the density of the solid isomorphs is several percent lower than the density of melting. These deviations may come from the fact that only the first-order terms in the Taylor expansion were included in the analysis.^{25} We leave this for future investigations.

Figures 14(a)–14(c) show the entropy of fusion (Δ*S*_{m}) for: Ne, Kr, and Xe, respectively. The accuracy is comparable to that of Ar but slightly worse than that of the LJ model [Fig. 14(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.

## VIII. CONCLUSION AND OUTLOOK

In summary, we have investigated the solid–liquid coexistence of Ne, Ar, Kr, and Xe using the SAAP. We conclude that the isomorph theory of melting, which has no free parameters except those determined by simulations at the reference temperature, gives accurate predictions for the coexistence line and of variations of properties along this line. An obvious question is how well the isomorph melting theory describes molecular fluids. As a possible starting point, Hellmann provided an analytical potential for the CO_{2} molecule based on accurate *ab initio* dimer calculations.^{87} Empirical data suggest that the isomorph theory applies to molecular systems,^{5,9,55,64,88,89} but the isomorph theory of melting^{25} has not yet been tested for such systems.

## SUPPLEMENTARY MATERIAL

The supplementary material contains details about the simulation algorithms, raw data, and additional figures.

## ACKNOWLEDGMENTS

The authors thank Søren Toxværd, Lorenzo Costigliola, Thomas B. Schrøder, and Nicholas Bailey for their suggestions during the preparation of this manuscript. This work was supported by the VILLUM Foundation’s Matter Grant (No. 16515).

## DATA AVAILABILITY

The data used in this study are available at http://doi.org/10.5281/zenodo.3888373 and in the supplementary material.

## REFERENCES

_{2}O), and hydrogen (H

_{2})