# Freezing of Simple Liquid Metals

## 1 Introduction

A microscopic description of freezing/melting phenomena presents a formidable challenge to present-day condensed matter theory . Among the simple liquids, liquid metals have received limited attention in this field, stemming from the relative complexity of their interatomic potentials and thermodynamic properties . Previous related studies include those of Stroud and Ashcroft , who addressed the melting of simple metals using variational principles; Iglói et al. , who studied the freezing of liquid aluminum and magnesium by means of density-functional (DF) methods; and Moriarty and coworkers , who calculated the structural stabilities of third-period simple metals using first-principles techniques.

In this contribution, we report on a new study of the freezing phenomena of simple metals that exploits recent advances in thermodynamic perturbation theory for solids. In line with Ref. , and most modern theories of freezing, our approach incorporates classical DF theory  in calculating the Helmholtz free energy. Although formally rigorous, DF theory is limited in practice by insufficient knowledge of the free energy functional for most except hard-body systems. To circumvent this limitation, we apply Weeks-Chandler-Andersen (WCA) perturbation theory , which splits the effective interatomic pair potential into reference and perturbative parts and maps the reference system onto an effective hard-sphere (HS) system. For the inhomogeneous solid phase, DF theory provides an adequate approximation for the HS free energy, while an accurate model of the HS radial distribution function yields directly the first-order perturbation free energy. The homogeneous liquid phase is treated within the same perturbative framework, using essentially exact expressions for the HS free energy and pair distribution function. This treatment of the liquid and solid by the same theoretical approach constitutes an important requirement for a consistent description of phase coexistence.

The next section outlines the theoretical methods used to calculate effective pair potentials and free energies of simple metals. In Sec. 3 we present and discuss predictions of the theory for freezing transitions and relative stabilities of competing crystal structures. Finally, in Sec. 4 we summarize and conclude.

## 2 Theoretical Methods

Restricting attention to the simple metals, we construct effective density ($`\rho `$)-dependent interatomic pair potentials, $`\varphi (r;\rho )`$, and a volume-dependent, but structure insensitive, contribution to the binding energy, $`E_V(\rho )`$, via pseudpotential theory, assuming linear response of the conduction electrons to the ions. In practice we use empty-core pseudopotentials and the Ichimaru-Utsumi electron dielectric function, which are known to give reliable results for bulk liquid properties . Taking $`\varphi (r;\rho )`$ and $`E_V(\rho )`$ as the basic microscopic input, we proceed to calculate the thermodynamics of the liquid and solid phases by means of WCA perturbation theory. Thus, we separate $`\varphi (r;\rho )`$ at its first minimum into a short-range, purely repulsive, reference part, $`\varphi _\mathrm{o}(r)`$, and a long-range, weakly oscillating, perturbative part, $`\varphi _\mathrm{p}(r;\rho )`$, and map the reference system onto an effective HS system with a HS diameter prescribed by the well-known WCA criterion. This separation allows the ionic contribution to the Helmholtz free energy to be calculated within the coupling-constant formalism . To first order in $`\varphi _\mathrm{p}(r;\rho )`$, the resulting general (inhomogeneous) expression for the total Helmholtz free energy reads

$$F[\rho (𝐫)]=F_{\mathrm{HS}}[\rho (𝐫)]+\frac{2\pi N^2}{V}_0^{\mathrm{}}dr^{}r^2g_{\mathrm{HS}}(r^{};[\rho (𝐫)])\varphi _\mathrm{p}(r^{};\rho )+E_V(\rho ),$$
(1)   
where $`N`$ is the number of particles, $`V`$ the volume, and $`F_{\mathrm{HS}}[\rho (𝐫)]`$ and $`g_{\mathrm{HS}}(r;[\rho (𝐫)])`$ the free energy and radial distribution function (RDF), respectively, of the HS reference system, both functionals of the equilibrium one-particle number density $`\rho (𝐫)`$. The RDF of the solid is defined as an orientational and translational average of the two-particle density . For the liquid we use the homogeneous version of Eq. (1), taking for the HS free energy and RDF the Carnahan-Starling and Verlet-Weis parametrizations of computer simulation data .

To calculate the free energy of the HS solid we invoke classical DF theory . The DF approach is based on the existence of a functional $`[\rho (𝐫)]`$ of the density $`\rho (𝐫)`$ that satisfies a variational principle, according to which $`[\rho (𝐫)]`$ is minimized – for a given average density and external potential – by the equilibrium density, its minimum value equaling the Helmholtz free energy $`F`$. In practice, $`[\rho (𝐫)]`$ is split into an (exactly-known) ideal-gas contribution $`_{\mathrm{id}}[\rho (𝐫)]`$ and an excess contribution $`_{\mathrm{ex}}[\rho (𝐫)]`$, the latter depending entirely upon internal interactions. Here we approximate $`_{\mathrm{ex}}[\rho (𝐫)]`$ by the modified weighted-density approximation (MWDA) . This approximation maps the excess free energy per particle of the solid onto that of a corresponding uniform fluid of effective density $`\widehat{\rho }`$, according to

$$_{\mathrm{ex}}^{\mathrm{MWDA}}[\rho (𝐫)]=Nf_{\mathrm{HS}}(\widehat{\rho }),$$
(2)   
where the effective density, defined as

$$\widehat{\rho }=\frac{1}{N}d𝐫d𝐫^{}\rho (𝐫)\rho (𝐫^{})w(|𝐫𝐫^{}|;\widehat{\rho }),$$
(3)   
is a self-consistently determined weighted average of $`\rho (𝐫)`$. The weight function, $`w(r)`$, is specified by normalization and by the requirement that $`_{\mathrm{ex}}^{\mathrm{MWDA}}[\rho (𝐫)]`$ generates the exact two-particle direct correlation function in the uniform limit (see Ref.  for details).

Practical calculation of $`_{\mathrm{HS}}[\rho (𝐫)]`$ and $`g_{\mathrm{HS}}(r;[\rho (𝐫)])`$ requires specifying the solid density, i.e., the coordinates of the lattice sites and the shape of the density distribution about these sites. Here we consider the hcp, fcc, and bcc crystals with the density distribution modeled by the usual Gaussian ansatz, introducing a parameter $`\alpha `$ determining the width of the distribution. This parametrization allows the ideal contribution to the free energy functional to be accurately approximated by $`_{\mathrm{id}}/N=(3/2)k_\mathrm{B}T\mathrm{ln}(\alpha \mathrm{\Lambda }^2)5/2`$, where $`\mathrm{\Lambda }`$ is the thermal de Broglie wavelength, and yields the HS free energy as the minimum with respect to $`\alpha `$ of the approximate functional $`_{\mathrm{HS}}[\rho (𝐫)]=_{\mathrm{id}}[\rho (𝐫)]+_{\mathrm{ex}}^{\mathrm{MWDA}}[\rho (𝐫)]`$.

For the RDF appearing in the perturbation contribution to the free energy \[Eq. (1)\] we adopt the approach of Rascón et al. . Expressing $`g_{\mathrm{HS}}(r;[\rho (𝐫)])`$ as a sum over coordination shells, this scheme approximates the distributions of the second and higher coordination shells in simple mean-field fashion, but parametrizes the first-shell distribution in such a way as to incorporate nearest-neighbour correlations. The first-peak parameters are determined by imposing sum rules, e.g., the virial theorem and coordination number, so as to specify the contact ($`r=d`$) value and shape (area and first moment) of the first peak (for details see Ref. ). For a given solid structure, the lattice distances, the coordination number, the $`\alpha `$ that minimizes $`_{\mathrm{HS}}[\rho (𝐫)]`$, and the HS pressure, derived directly from $`F_{\mathrm{HS}}[\rho (𝐫)]`$, then combine to determine $`g_{\mathrm{HS}}(r;[\rho (𝐫)])`$ and hence the perturbation free energy.

Following the above procedure, we obtain the ionic free energy, $`F_{\mathrm{ion}}`$, which is the only structure-dependent contribution to the total free energy, and which suffices for assessing relative stabilities of different solid structures. To fully describe freezing/melting, however, it is essential to consider the total free energy, $`F=F_{\mathrm{ion}}+E_\mathrm{V}`$, where $`E_\mathrm{V}`$ is an additional contribution from the conduction electrons. The electronic free energy (or volume energy), $`E_\mathrm{V}`$, although independent of structure, depends on the average density and thus influences the densities of coexisting phases. Here we take for $`E_\mathrm{V}`$ the usual expression following from second-order perturbation theory  and consistent with the effective pair potential.

## 3 Results and Discussion

Figure 1 illustrates the density dependence of the effective HS diameter, $`d`$, for three different crystal structures of Al near its melting point. According to the WCA prescription, $`d`$ depends both on the reference pair potential and on the first peak of the HS pair distribution function. Since the first peak is identical for hcp and fcc crystals, $`d`$ is the same for these two close-packed structures. Differences in the first peak for the bcc crystal result in a smaller diameter for that more open structure. Because of the sensitivity of $`F_{\mathrm{HS}}`$ to the HS packing fraction, these structural variations in $`d`$, though only of the order of 1 %, can amount to significant variations in the free energy. In passing, we note that the widely-used Barker-Henderson prescription for $`d`$  neglects such structural dependencies and that only knowledge of the HS RDF permits consistent application of the more accurate WCA prescription .

Predicted structural free energy differences for Na, Mg, and Al are displayed in Fig. 2. Evidently the theory correctly predicts the observed equilibrium structures for these three metals, namely bcc-Na, hcp-Mg, and fcc-Al, although the density at crossover from fcc\- to bcc-Na is somewhat overestimated. Compared with predictions of first-principles approaches, our results agree qualitatively well. For Mg and Al the stability trends of Fig. 2 match those from generalized pseudopotential theory (GPT) and the linear muffin-tin orbitals (LMTO) method . As well, for Na we predict the same trend as GPT, while LMTO fails for this case. Compared with available experimental data, our structural free energy differences are in somewhat better quantitative agreement than those of either GPT or LMTO. For example, for Mg and Al we predict differences between hcp and fcc free energies of $`1.6`$ and $`3.2`$ mRy, respectively, compared with $`0.6`$ and $`1.7`$ mRy from GPT, $`0.7`$ and $`1.6`$ mRy from LMTO, and $`1.5`$ and $`4.2`$ mRy from experiment (extrapolations of thermochemical alloy data, as quoted in Ref. ). It should be emphasized, however, that while the first-principles methods calculate ground-state ($`T=0`$) energies, thermodynamic perturbation theory yields free energies of finite-temperature systems. We note further that application of perturbation theory is restricted to those structures for which the HS system is at least metastable, i.e., for which $`_{\mathrm{HS}}[\rho (𝐫)]`$ has a local minimum. This excludes, e.g., the diamond structure of Si, the HS diamond crystal being unstable.

Finally, we investigate the freezing/melting transition by comparing predicted free energies for the liquid and solid phases. Figure 3 exemplifies the comparison for Al near its observed melting temperature. Crossing of the ionic free energy curves for the liquid and fcc phases implies a distinct transition near the observed equilibrium solid density. Similar comparisons for Na and Mg also demonstrate liquid-solid phase transitions in qualitative agreement with observation. In this respect, the effective-pair-potential approach reasonably describes the many-body interacting system. Note that the liquid-solid crossover density serves as a useful upper (lower) bound on the liquid (solid) density at coexistence. Adding the volume-energy contribution – identical for the two phases – while not affecting the crossover density, does alter the shape of the curves. As a result, the densities of coexisting liquid and solid phases, as determined by a Maxwell common tangent construction, depend on $`E_\mathrm{V}`$. An accurate volume energy is known to be essential for calculating realistic pressures and bulk moduli of metals from pseudopotential theory at the level of effective pair interactions . Similarly, we find the liquid-solid coexistence analysis to be sensitive to the form of $`E_\mathrm{V}`$, particularly its curvature with respect to $`\rho `$. In fact, the sensitivity is such as to preclude here a reliable determination of coexistence densities. Given the reasonable predictions for ionic free energies, this apparent limitation of the theory may reflect a need to go beyond the linear response approximation to include higher-order response effects in the volume energy (see also , p. 42ff). The relative significance of such effects, as well as distinctions  between the real-space representation of $`E_\mathrm{V}`$ used here and the reciprocal-space representation of Ref. , are interesting issues deserving of further attention.

## 4 Conclusions

Summarizing, working within an effective-pair-potential framework, we have implemented a consistent form of thermodynamic perturbation theory, one that is of comparable accuracy for both the liquid and solid phases, and used it to study structural stabilities and freezing behaviour of third-period simple metals. Predictions for structural free energy differences are at least as accurate as those from two first-principles methods and are in generally good qualitative agreement with experiment. On the basis of ionic free energies, freezing/melting transitions are predicted near the observed densities. As the densities of coexisting liquid and solid phases are found to depend sensitively on the electronic free energy, a complete coexistence analysis has not been attempted pending more accurate knowledge of this part of the free energy. The computationally practical approach demonstrated here should prove useful in future studies of the freezing of liquid metal mixtures, as well as in guiding more sophisticated ab initio simulations.

Acknowledgements

ARD acknowledges the Forschungszentrum Jülich for use of its computing facilities. GK and JH acknowledge the support of the Österreichische Forschungsfonds under Proj. Nos. P11194 and P13062 and of the Oesterreichische Nationalbank under Proj. No. 6241.