Freezing of He-4 and its liquid-solid interface from Density Functional Theory

We show that, at high densities, fully variational solutions of solid-like type can be obtained from a density functional formalism originally designed for liquid 4He. Motivated by this finding, we propose an extension of the method that accurately describes the solid phase and the freezing transition of liquid 4He at zero temperature. The density profile of the interface between liquid and the (0001) surface of the 4He crystal is also investigated, and its surface energy evaluated. The interfacial tension is found to be in semiquantitative agreement with experiments and with other microscopic calculations. This opens the possibility to use unbiased DF methods to study highly non-homogeneous systems, like 4He interacting with strongly attractive impurities/substrates, or the nucleation of the solid phase in the metastable liquid.


I. INTRODUCTION
Helium crystals represent a fascinating system in which general properties of crystalline surfaces such as the equilibrium crystal shape, surface phase transitions (roughening) and elementary mechanisms of the crystal growth can be studied over a wide temperature (T ) range, in principle down to absolute zero. Helium crystals can be grown chemically and isotopically pure, and with very few lattice defects. Helium crystal growth is fast in comparison to that of other solids and it is characterized by the unique -quantum-phenomenon of crystallization waves, i.e., melting-freezing waves which can easily propagate on the liquid-solid interface at low temperatures. This allows an accurate measurement of the surface stiffness γ and surface tension σ and of their anisotropy [1]. The research on the surface of helium crystals has been recently reviewed [2].
A microscopic approach to the study of solid 4 He and liquid-solid coexistence at low temperature needs a fully quantum approach. Accurate Path Integral Monte Carlo [3] (PIMC) in the 5-30 K temperature range and Green Function Monte Carlo [4] (GFMC) calculations on solid 4 He have been reported in the past. More recently, the liquid-solid interface has been adressed within a Variational Monte Carlo (VMC) approach using shadow wavefunctions [5].
Density Functional (DF) methods [6] represent a useful computational tool to study the properties of quantum inhomogeneous fluids, especially for large systems where DF provides a good compromise between accuracy and computational cost. Indeed, the T = 0 properties of inhomogeneous liquid 4 He can be accurately described within DF theory (DFT) by using the phenomenological functional proposed in Ref. [7] and later improved in Ref. [8], and similar ones. They have been widely used in a variety of problems involving inhomogeneous 4 He systems like, e.g., liquid-vapor interface [7,8], pure and doped clusters [8,9], layering and prewetting transitions in films [8], alkali atom adsorption on the surface of liquid 4 He and droplets [10], vortices in 4 He clusters [11], etc.
In view of its conceptual interest on the one hand, and of the potential applications on the other hand, we have undertaken a fully variational study of the liquid-solid transition and coexistence of 4 He at zero temperature within DFT. Previous attempts on 4 He [12,14,15,16] or quantum hard spheres [17,18] have been reported. They adapted the techniques initially designed to address the liquid-solid transition in classical systems (for a review see, e.g., Ref. [13]). They use trial functions (sums of gaussians) to parametrize the solid density, and resort to a second-order expansion of the energy around the liquid density (Ramakrishnan-Yussouff (RY) approach) [12,14,15,16], or to a modified weighted density approximation exact to second order for small perturbations of the liquid [17,18]. These methods need the density-density static linear response function χ of the liquid as input. In helium, these approaches have been found to be, at most, a reasonable starting point to qualitatively describe the freezing transition [12]; moreover, some of these attempts always give the liquid [14] or the solid [15,16] as the stable phase. Additional problems arise in treating the liquid-solid interface, because it is found to be too narrow to allow the use of coarse graining methods [16]. Besides, the use of trial function density profiles makes all these approaches of limited applicability for systems characterized by the presence of both solid and liquid phases.

II. THE DENSITY FUNCTIONAL FOR SOLID 4 HE
We attempt here a fully variational DF description of both liquid and solid phases on the same footing. Our starting point is a simplified version of the Orsay-Trento (OT) functional [8], where the energy of the 4 He system is written as The first term is the usual quantum kinetic energy term, the second term contains a Lennard-Jones He-He pair potential V LJ (r) screened at distances shorter than a characteristic length h. In the third and fourth terms, the weighted densityρ is the average of the density over a sphere of radiush. These terms account phenomenologically for short range correlations. The parameters h, c 2 and c 3 (in the original formulation, h =h) are fixed to reproduce the experimental density, energy per atom and compressibility for the liquid at zero temperature and pressure [8]. The original non-local kinetic energy term proportional to ∇ρ(r) · ∇ ′ ρ(r ′ ) [8] has been dropped, because, as one might expect, it leads to unavoidable instabilities in the solid phase. Its neglecting causes the static response function χ of the liquid to be only qualitatively described. We should mention the existence of an alternative parametrization of this term that still reproduces the experimental static function, while being free from instabilities [19]. As the resulting DF stems from a somewhat different strategy than current DF's do, we have preferred to keep using an OT-like functional, introducing the changes mentioned below. The equilibrium density ρ(r) is obtained by minimizing E[ρ] with respect to density variations, subject to the constraint of a constant number of 4 He atoms N . This is achieved by evolving in the imaginary time a non-linear Schrödinger equation for the order parameter Ψ(r) ≡ ρ(r), where the Hamiltonian is given by . The effective potential U is defined in terms of the variational derivative of the last three terms in the energy functional Eq. (1) with respect to Ψ [8].
The calculations are performed in a periodically repeated supercell containing a fixed number of 4 He atoms, and the starting configuration is a superposition of gaussian profiles with arbitrary but small width, placed in the positions of the hcp structure, which is the experimental solid structure of 4 He. Interestingly, for average densities corresponding to the stable liquid phase, the equilibrium density obtained is homogeneous, whereas for higher densities -those roughly corresponding to the solid-a solid 4 He structure, characterized by a periodic density distribution with helium atoms at the lattice sites of an hcp crystal, is found to be the stable phase. 'Solid' means here a highly inhomogeneous configuration characterized by regions, called 'atoms' for brevity in the following, where the density is very large. Such atoms are only slightly overlapping with the density tails of neighboring atoms. As an example, we show in Fig. 1 the calculated density profile for one such solid structure, computed at the experimental freezing density ρ f = 0.0287Å −3 .
We have studied other ordered structures as well. An fcc lattice is also found to be a stable phase at typical solid densities, its energy being only slightly higher than that of the hcp phase (by 0.02 K per atom). A simple cubic structure, at the same densities, is found instead to be unstable towards the homogeneous (liquid) phase. We want to point out that the correct solid structure is found even if the initial gaussian superposition has not the correct parameters. For instance, we found that by starting with gaussians, placed on a hcp lattice with a lattice constant twice the equilibrium one, but such that the total normalization gives the correct 4 He density, the equilibrium hcp structure is eventually recovered at the end of the minimization process.
In spite of this encouraging result, a close analysis of the density dependence of the total energy of the hcp solid 4 He, calculated using Eq. (1), shows that it largely deviates from the experimental result. The reason is an unphysical, large density pile-up in the core region of the 4 He atoms, such that the resulting solid is too dense. To fix this unrealistic behavior, we need a mechanism  [21]; the experimental EOS for the liquid is not shown since it coincides, by construction, with that of the OT functional, whose agreement with experiment is excellent [8].
that makes energetically costly such density pile-up. The simplest one is to add a 'penalty' energy term to the functional, which has the following form: Here f [ρ(r)] is a 'switch' function which becomes appreciably different from zero only when the density is larger than a predefined value ρ m . One such function is where C, β and ρ m are DF parameters. The effect of E p [ρ] is to add to the effective potential U [ρ] the term If C > 0 and when ρ(r) ρ m , this term acts as a repulsive barrier, which forbids extra pile-up of the density. We can thus use C, β, and ρ m as adjustable parameters in such a way to get agreement with the experimental equation of state (EOS) for the solid phase. Note that there is a large freedom in the choice of the penalty term since, by construction, it has no effect whatsoever on the liquid structure, even in a highly structured liquid, because in practice ρ m is a factor 10 than typical liquid densities. We have found that the choice C = 0.1 Hartree, β = 40 A 3 , and ρ m = 0.37Å −3 [20], yields an EOS for the solid in good agreement with experiments [21], as shown in Fig. 2.

III. THE FREEZING TRANSITION AND THE LIQUID-SOLID INTERFACE
Once the EOS of the homogeneous liquid and the EOS of the solid are determined, one can easily locate the freezing transition by means of a double-tangent Maxwell construction. We find a freezing pressure P f = 26.1 bar and a chemical potential at coexistence µ f = 0.69 K. Our calculated values for P f and the solid and liquid densities at coexistence are reported in Table I, together with those obtained by other authors.
We now turn to the liquid-solid interface of 4 He, characterized by the interfacial tension σ and width ξ. 4 He is an example of a system for which σ has been directly measured, thanks to the possibility of propagating crystallization waves at its surface. The value of σ is 0.17 × 10 −3 N/m [1], but ξ is not known experimentally.
The 4 He liquid-solid interface is an example of a system whose direct simulations by means of exact Monte Carlo methods can be computationally very expensive. Recently, it has been tackled by means of shadow wavefunctions VMC [5], and also by the rescaled DF we have referred to before [16]. We have used the modified DF, Eqs. (1)(2), to study the liquid-solid interface. To describe the coexisting phases, we start from a slab of hcp crystal 4 He in contact with a region of liquid 4 He, and minimize the energy functional to find the equilibrium configuration, imposing that during the minimization the chemical potential is constant and equal to the value µ f = 0.69 K found from the double-tangent construction illustrated in Fig. 2. We consider explicitly the interface between liquid 4 He and the (0001) surface of the crystal (basal plane). Fig. 3 shows a view of the equilibrium density configuration at the liquid-solid interface. A related quantity, the 4 He density averaged over a plane parallel to the basal plane, is shown in Fig. 4. The thickness ξ of the liquid-solid interface has been estimated by evaluating the 10-90% width of the curve interpolating the maxima in the density profile shown in Fig. 4. We find 9.3Å, in close agreement with the VMC results. Our findings are compared in Table II with other experimental and theoretical results.
From the calculated density profiles we estimate the liquid-solid surface tension as σ = (E −µ f N +P f V )/(2A). Here V is the volume of the supercell, and A is the surface area of the exposed (0001) face. A factor 1/2 appears in the equation to account for the two free surfaces delimit- ing the solid film in our slab geometry. Our calculations, which are three dimensional in nature, have been done by using a supercell accomodating ∼ 6 layers of solid 4 He and a thick liquid layer in contact with it. We find that the convergence with respect to the amount of liquid present in the supercell is rather slow (to obtain σ accurately, one should have two very wide liquid and solid regions in contact). Our estimated value is σ = 0.1×10 −3 N/m [23], which is a rather good result, of quality similar to the VMC value [5]. Method σ (10 −3 N/m) ξ (Å) Experiment [1] 0.17 -VMC [5] 0.25± 0.1 10-12.5 DFT [16] 0.47 5.6 Present work 0.1 9.3

IV. CONCLUSIONS
We have shown that DF's currently used for liquid 4 He admit fully variational solutions of highly inhomogeneous type that may represent localized helium 'atoms'. With a simple modification, one such DF is found to accurately describe both liquid and solid phases of 4 He and their coexistence. Our scheme does not rely on any specific trial function density profile for the solid phase, but gives an unbiased description of both phases on the same footing. This broadens the range of applicability of current DF, permitting to study in an unbiased way highly non-homogeneous 4 He systems, like droplets doped with strongly attractive impurities, as for instance alkali ions [24], or 4 He on strongly attractive substrates, such as graphite or carbon nanotubes [25]. Another potentially useful application is the study of the nucleation of the solid phase in the metastable superfluid [26].
We thank Franco Dalfovo for useful discussions. F. A. aknowledges funding from CESCA-CEPBA, Barcelona, through the program HPC-Europa Transnational Access.