Solid–liquid coexistence of neon, argon, krypton, and xenon studied by simulations

The thermodynamic and transport properties of condensed matter systems are determined by their potential-energy functions.¹ For a class of systems, the potential-energy function exhibits “hidden scale invariance,”² 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.²⁴ 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.²⁵ 

Specifically, we investigate below the simplified ab initio atomic potential (SAAP) recently suggested by Deiters and Sadus.²⁴ 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.²⁸ This approach is considered the gold standard of quantum chemistry²⁹ 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₁, r₂, r₃, …, 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) $≡$$v$(xσ)/ε is given by

The six a parameters²⁴ are determined by fitting to the results of the above-mentioned ab initio calculations on dimers.²⁶ 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.³⁴ Pressure corrections are estimated to range from −0.3ε/σ³ to −0.1ε/σ³ 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.²⁴ 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⁻¹² − x⁻⁶] (green dashed line; simulation details for the LJ system can be found in Ref. 25). As noted by Thiel and Alder,³⁵ 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⁵ exp(−12x)^{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³⁹ (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.

In order to illustrate the theory, this paper starts by treating in detail the case of argon (Secs. II–IV). Sections V–VII report the analogous results for neon, krypton, and xenon. Section VIII gives a brief summary.

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₀ = 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₀ = 22.591(4)ε/σ³ at T₀ 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.⁵⁰ The required slopes dp/dT are computed from isobaric simulations of a solid and a liquid using the Clausius–Claypeyron relation:⁵¹ 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⁵¹ 

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.⁵¹ 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⁴⁰ 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,⁵² 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 $ρ¯=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), ν/εσ⁹ = 0.068 753 6 is the Axilrod–Teller–Muto parameter⁵³ of Ref. 52, and λ is a parameter obtained from analyzing two-body and three-body simulation data.⁵⁴ 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.

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_BT.^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³ 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⁸ 

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₀, ρ₀) given in Table I.

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.¹⁸ 

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”⁷ 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⁻ⁿ), 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.⁶⁶ For the dense phases (liquid and solid), this results in the expression

where $v$⁽ⁱ⁾(r) is the ith 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₀. 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.⁶⁸ 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¹⁵ 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³⁷ located where γ = 2.33(3).⁷⁰ 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.⁷⁴ 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 $βVex≡(∂W/∂T)V/V$ and the excess isochoric heat capacity per volume $cVex≡(∂U/∂T)V/V$⁠,⁸ 

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–77c_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 ≃ 3k_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.⁷⁹ 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.⁶⁷ (this EOS is implemented into the CoolProp software library by Bell et al.⁸⁰). 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²⁵ to SAAP argon.

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₀. This leads²⁵ 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₀ = 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_BT 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.²⁵ 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⁸² (⟨u²⟩ = Δr²/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.0365k_B/ε needed in the theoretical prediction.²⁵ 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)].

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.⁵² 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,⁸⁴ 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).

As an aside, we investigate the validity of the Simon–Glatzel equation for the coexistence pressure,⁸⁶ 

We first use T_ref = T₀ and p_ref = p₀. 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$≅$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.

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.⁷ 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.

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.²⁵ 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.

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₂ molecule based on accurate ab initio dimer calculations.⁸⁷ Empirical data suggest that the isomorph theory applies to molecular systems,^{5,9,55,64,88,89} but the isomorph theory of melting²⁵ has not yet been tested for such systems.

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

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).

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