The Solid Phase of 4He: A Monte Carlo Simulation Study

The thermodynamics of solid (hcp) 4He is studied theoretically by means of unbiased Monte Carlo simulations at finite temperature, in a wide range of density. This study complements and extends previous theoretical work, mainly by obtaining results at significantly lower temperatures (down to 60 mK) and for systems of greater size, by including in full the effect of quantum statistics, and by comparing estimates yielded by different pair potentials. All the main thermodynamic properties of the crystal, e.g., the kinetic energy per atom, are predicted to be essentially independent of temperature below ∼ 1 K. Quantum-mechanical exchanges are virtually non-existent in this system, even at the lowest temperature considered. However, effects of quantum statistics are detectable in the momentum distribution. Comparison with available measurements shows general agreement within the experimental uncertainties.


Introduction
Among naturally occurring substances, the solid phase of helium (which is stable only under applied pressure) represents a unique example of "quantum crystal", i.e., one in which quantum mechanics most considerably and observably affects its thermodynamic properties [1].The most direct manifestations include a kinetic energy per particle markedly above its classical value 3T/2 (we adopt a system of units in which the Boltzmann constant k B = 1); quantitative theoretical predictions for such excess kinetic energy have been made [2] for all rare gas solids, and there exist a number of experimental confirmations (see, for instance, Ref. [3]).For the rare gas solids of elements heavier than helium, however, the excess kinetic energy is essentially the only measurable manifestation of their quantummechanical nature.Other relevant physical properties (e.g., the momentum distribution) can be largely understood along classical lines.
In solid helium, on the other hand, the classical picture is quantitatively altered by zeropoint motion, which is nowhere near as important in other solids (with the sole exception, albeit not to the same extent, of molecular hydrogen [4].This includes enhanced atomic excursions away from the equilibrium (lattice) positions, with respect to what is predicted and observed in most solids, based on classical (thermal) arguments.Indeed, for a long time, zero point motion was believed to be the physical reason underlying the failure of liquid helium to solidify, under the effect of its own vapor, a view that has relatively recently been challenged, as quantum statistics (i.e., exchanges of indistinguishable particles) has been shown to play a pivotal role in the stabilization of the crystal phase, at least in 4 He [5].In solid 3 He, exchanges have been proposed as the physical agent responsible for the stabilization of the bcc over the hcp crystalline phase [6].
As helium can be regarded as the archetypal quantum crystal, it poses a challenge to many-body theorists; achieving an accurate microscopic description, based on realistic interatomic potentials, capable of making reliable predictions for experimentally cogent quantities, is a worthwhile goal in theoretical condensed matter and many-body physics.Of particular interest are calculations in which the fundamental quantum-mechanical equations are treated without any uncontrolled approximations.Ideally, one starts from a microscopic Hamiltonian in which the only external (independent) input is the potential energy of interaction among atoms.
Experimentally, the bulk of the available information on the solid phase of helium has been delivered by x-ray [7][8][9][10] and neutron [11][12][13][14][15][16][17] scattering measurements, whose main outcome is the thermodynamic equation of state, as well as single-atom dynamics through the momentum distribution, from which the mean atomic kinetic energy can be inferred [18].On the theoretical side, Path Integral Monte Carlo simulations have played a major role in shaping our current understanding of the physics of condensed helium in the condensed phase [19,20].These are unbiased calculations which involve no a priori physical assumptions.Indeed, only two main assumptions are built into the (almost) totality of PIMC simulations, namely (a) that He atoms, which are taken as the elementary constituents (see discussion below) can be regarded as distinguishable particles, as quantum-mechanical exchanges are suppressed in the solid phase, due to atomic localization and (b) that the bulk of the potential energy of interaction among the helium atoms can be captured by a pair-wise, central potential featuring a repulsive core at short interatomic separation; three-body interaction terms are neglected.
The development of an accurate, ab initio pair potential, describing the interaction of two helium atoms, has been the subject of an intense theoretical effort spanning nearly a century and still ongoing [21][22][23][24].The first pair interaction to have been extensively used in PIMC simulations whose results have been compared to available experimental data for the superfluid phase of 4 He [25] is the earliest version of the Aziz pair potential [22].The agreement between theory and experiment afforded by that potential is generally satisfactory, as confirmed by subsequent studies [26,27] overcoming the main limitations of the original calculation, mainly the relatively small size of the simulated system (64 atoms).
For the solid phase, on the other hand, the comparison of experimental and theoretical estimates has been considerably more limited in scope, for a number of reasons.The most important is the simple fact that the wealth of experimental data available for the superfluid phase [28] is presently unmatched in the crystal.The resurgence of interest in the solid phase of helium that took place two decades ago, following the report of possible "supersolid" behavior [29], prompted new investigations, significantly expanding the available comparison dataset.However, assessing different theoretical calculations is complicated by the plethora of helium pair potentials [24,[30][31][32][33][34][35][36] that have been proposed and utilized since Ceperley and Pollock's 1986 calculation.
Many of these potentials are either refinements of, or at any rate based upon the original Aziz model interaction; the differences among them are often relatively small and pertain to specific features, such as the depth of the attractive well, or the behavior of the repulsive core at short interatomic separation.To be sure, some of these aspects become quantitatively important as the system is compressed, and for pressures, e.g., in the megabar range, three-body interactions must be included in order to achieve a quantitative reproduction of the experimental EOS [37,38].It is not clear, however, how sensitive to the details of the pair interaction are the most important experimental quantities, e.g., the single-particle kinetic energy, especially at low pressure (less than ∼ 100 bars).
We present here the results of PIMC simulations of the hcp phase of 4 He at low temperature (typically around 1 K, but also as low as ≲ 0.1 K for specific densities), in a range of pressure between 25 and 150 bars.We consider a perfect helium crystal, namely one free of point (e.g., vacancies or interstitials) or extended (e.g., dislocations) defects.We make use of different versions of the Aziz pair potential and compare the results for the single-particle kinetic and total energy per atom and for the pressure to the most recent experimental estimates.For the physical quantities considered here, agreement with available experimental data seems satisfactory; the differences between the estimates yielded by the different potentials are typically (much) smaller than the experimental uncertainty.In other words, in the region of parameter space considered here, all of the pair potentials that have been utilized over the past two decades give essentially equivalent results.This study also provides additional, strong evidence that a perfect crystal of 4 He is not superfluid, at arbitrarily low temperatures.Indeed, quantum-mechanical exchanges of atoms remain strongly suppressed even at a temperature at low as 60 mK, at the melting density.Nevertheless, the effects of quantum statistics are detectable in the one-body density matrix (and therefore, experimentally, in the momentum distribution).We also compare our low-temperature energy estimates to recent ground-state Monte Carlo calculations; the comparison again points to the greater reliability of finite temperature techniques.
This paper is organized as follows: in section 2 we describe the model of interest and briefly review the computational technique utilized; in section 3 we present our results, offering a discussion in section 4.

Model and Methodology
The system is described as an ensemble of N point-like, identical particles with a mass m equal to that of a 4 He atom, and with spin zero, thus obeying Bose statistics.This is where a basic assumption is built into the model, namely the Born-Oppenheimer approximation, which allows us to integrate out the electrons and regard 4 He atoms as elementary particles.At ordinary conditions of temperature and pressure, this procedure is justified by the large decoupling between electronic and ionic energy scales.In different physical settings, e.g., in the interior of Jovian planets, such separation is no longer possible, and one has to take into account electronic and ionic degrees of freedom separately [39], but, in this study, we regard the individual 4 He as our elementary constituents.
The system is enclosed in a cell of volume Ω so that n = N/Ω is the nominal density.The cell is shaped like a cuboid, with periodic boundary conditions in all three directions.The sizes are adjusted to fit a perfect classical crystal, with no defects.
The quantum-mechanical many-body Hamiltonian reads as follows: where the first (second) sum runs over all particles (pairs of particles), λ ≡ h2 /2m = 6.0596415KÅ 2 , r ij ≡ |r i − r j | and v(r) is the pair potential which describes the interaction between two helium atoms.As mentioned in the Introduction, several different forms for v have been used in this study.Most of the results shown here were obtained with the original (1979) Aziz pair potential [22] (henceforth referred to as Aziz I), which is largely phenomenological, the result of a careful combination of all available theoretical and experimental data for the interaction of two helium atoms in the gas phase.The Aziz I pair potential has been shown to provide a rather accurate quantitative description of the superfluid phase of 4 He [19].Generally speaking, at high density (in both the liquid and solid phases), the pressure obtained from calculations based on the Aziz I potential is higher than the experimentally measured one, signaling that the Aziz I potential is too "stiff" at short distances and/or the presence of attractive forces, presumably arising from the interaction of clusters of atoms (e.g., triplets).Two different strategies have been pursued, in order to improve the agreement between theory and experiment.In the effective potential approach [40], one attempts to incorporate three-body effects into a single pair potential, which is parameterized to reproduce the experimental equation of state.Alternatively, and this is the approach that has been most extensively considered in recent times, the aim is that of moving away from empirical, parametrized forms, utilizing instead expressions arising from ab initio quantum chemistry calculations.In that case, three-body terms are included as separate contributions, and this is the spirit in which the original Aziz potential has been subsequently refined and other versions proposed.The result has been typically that of slightly worsening the agreement between theory and experiment, which would then be improved by the explicit inclusion of three-body terms [41,42].As mentioned above, however, no systematic comparison has been carried out of different versions of the Aziz pair potential when it comes to reproducing the physical properties of the solid phase of helium at low pressure, where the effect of three-body interactions is expected to be small.Here, for comparison purposes we are going to show results obtained with the Aziz I potential, as well as with the versions thereof proposed in 1987 [30] (Aziz II) and 1995 [33] (Aziz III).
We performed QMC simulations of the system described by Eq. ( 1), based on the canonical version [43,44] of the continuous-space Worm Algorithm [26] based on Feynman's space-time approach to quantum statistical mechanics [45].Since this methodology is extensively described in the literature, it will not be reviewed here; rather, all the relevant technical details will be provided, enabling others to attempt to reproduce the results presented here.The first aspect to discuss is that in all the simulations carried out in this study, 4 He atoms are treated as indistinguishable, i.e., the many-particle paths in imaginary time are allowed to "entangle", with particles trading places [19,45].Previous attempts to assess quantitatively, by computer simulation, the importance of quantum statistics in a perfect crystal of helium concluded that the effect was negligible [46,47]; in this work, we extend those studies to lower temperatures.We note that the WA has been shown particularly effective at sampling permutations of identical particles, and therefore we may be reasonably confident that failure to detect a significant presence of exchanges should constitute a true physical reflection of the nature of the system, as opposed to possible algorithmic inefficiency.Details of the simulation are standard; the short-time approximation to the imaginarytime propagator used here is accurate to fourth order in the time step τ (see, for instance, Ref. [48]).We have carried out numerical extrapolation of the estimates to the τ → 0 limit, and observed convergence of the thermal averages for a value of τ = 3.125 × 10 −3 K −1 for all quantities of interest here.Fig. 1 shows the extrapolation procedure for the kinetic energy per atom, e k .Standard estimators were used for all the energetic properties, including the well-known virial estimator for the pressure [19].
Physical quantities of interest calculated in this study include, besides pressure and energetics, structural correlations such as the pair correlation function, which is related to the experimentally accessible static structure factor.We have carried out simulations for systems of different sizes, the largest comprising N = 512 atoms.In general, we found that the numerical estimates of e K obtained on a simulated system comprising N = 216 atoms are indistinguishable, within statistical uncertainties, from those obtained on systems of larger sizes.The range of system sizes considered in this work gives us reasonable confidence in our ability to gauge finite-size effects, i.e., that the numerical values quoted herein are representative of the thermodynamic limit, within their statistical uncertainties.
For the calculation of the (potential) energy per particle and the pressure, we estimated the contribution of particles outside the main simulation cell by assuming a value of the pair correlation function g(r) = 1, for r greater than half the (shortest) cell side.Based on the observed quantitative consistency (within statistical uncertainties) of results for systems of different sizes, we contend that this procedure is numerically reliable.We come back to this point more quantitatively when discussing our results.

Effect of temperature
We begin by discussing the effect of temperature on the energetics and structure of solid 4 He.All of the calculations whose results are illustrated in this subsection were carried out using the Aziz I pair potential (1979).
The left (right) panel of Fig. 2 shows the computed kinetic energy per 4 He atom in a hcp 4 He crystal of density 0.0286 (0.0312) Å −3 , in the temperature range T < 1.5 K, the lowest temperature being T ≲ 0.1 K.The higher value of density corresponds to a pressure close to 56 bar; the lower value is very close to the T = 0 melting density and yields a computed pressure of approximately 25 bars (in both cases the pressure is nearly temperature-independent).Also shown for comparison are the most recent experimental estimates of the kinetic energy per atom in solid helium, namely those of Ref. [16].The two experimental works with which a comparison seems meaningful are Refs.[15] and [16].The comparison of theoretical and experimental estimates shows quantitative agreement, taking into account the relatively large uncertainty in the experimental determination of the kinetic energy, as well as the uncertainty on the quoted values of the density in Ref. [16] and the slight differences in temperature.
Within the statistical uncertainties of our calculation, no significant dependence on the temperature of the kinetic energy per atom can be detected in the temperature range considered.This is consistent with the observation made in previous theoretical [49] and experimental [16] work (see Fig. 2).We generally find this to be the case for all the basic energetic properties, namely kinetic and total energy per atom, as well as the pressure and for structural correlation as well, in the entire range of density considered in this work.Henceforth, therefore, unless otherwise stated, all of the results quoted here are at temperature T = 1 K, which can be for all practical purposes considered equivalent to the ground state for the helium crystal.We come back to a more comprehensive and quantitative discussion of our results for the energetics of the helium crystal in 3.3.Fig. 3 displays the spherically averaged pair correlation function g(r), computed by numerical simulations at temperature T = 1 K, for three different densities, namely n = 0.0286 Å −3 , which corresponds to a pressure of approximately 25 bars, n = 0.0312 Å −3 , for which the computed pressure is 56 bars, and the highest density considered here, namely n = 0.0353 Å −3 , for which the pressure is 141 bars.These are results obtained using the Aziz I potential, but as we show more quantitatively below, the differences in the values of the pressure computed with different versions of the Aziz potential are relatively small (of the order of a fraction of a bar), in the range of densities considered here.These values of the pressure (and all others reported in Table 1 below) are in excellent quantitative agreement with the experimental equation of state of Driessen, van der Polls and Silvera [50,51].
The calculations whose results are shown in Fig. 3 were performed on a system comprising N = 216 atoms.As shown in Fig. 3, the value of the pair correlation function g(r) is very close to unity, giving us confidence that the approximation utilized to compute the contribution to the potential energy and to the pressure, namely assuming g(r) = 1 outside the simulation cell, is quantitatively reliable.One might think of such a procedure as being somewhat "crude", but its use is justified, in practice, by both the exponential decay of the fluctuations around unity of g(r), as well as the rapid (∼ 1/r 6 ) decay with distance of the pair potential.In order to obtain a numerical check, we have carried out calculations on a system of N = 512 atoms.and computed the pressure based on the same approximation.The results are all consistent with those obtained for a system of N = 216 atoms; for example, at a density n = 0.0312 Å −3 we obtain a pressure of 56.3(5) bars at T = 1 K, on a system of N = 512 atoms.

Effects of quantum statistics
As mentioned in the introduction, the calculations carried out in this work fully include the effect of quantum statistics, which means that exchanges of indistinguishable helium atoms are allowed.However, the frequency with which they occur remains exceedingly low, even at the lowest temperature and density considered here.For example, at n = 0.0286 Å −3 , which is very close to the T = 0 melting density, and at a temperature as low as T = 0.0625 K, the measured likelihood of a single 4 He atom to be involved in an exchange is less than one percent, three-particle exchanges occurring in the basal plane being overwhelmingly the most frequent.Within the statistical errors of our calculation, the effects of Bose statistics on the energetic properties of the helium solid are essentially not detectable.We carried out a separate simulation in which 4 He atoms are treated as boltzmannons, i.e., distinguishable quantum particles, and estimate the kinetic energy per atom to remain unchanged, within our statistical uncertainties, as a result of the neglect of quantum statistics.
That is not to say, however, that quantitative effects of quantum statistics cannot be observed, chiefly in the momentum distribution, which is experimentally measurable [15,17].Specifically, effects of Bose statistics have been observed in solid 4 He near melting, at a temperature as high as 1.6 K, through an enhancement of the momentum distribution at low momenta, compared to what one would predict based on a theoretical approach neglecting quantum statistics.
Our calculations confirm such predictions.Fig. 4 displays the spherically averaged one-particle density matrix n(r), which is defined as where ⟨...⟩ stands for thermal average, Ψ and Ψ † are the usual Bose field operators and the integration dα is over the solid angle of r.The one-body density matrix is connected to the momentum distribution through a Fourier transformation.The comparison of the one-body density matrix computed with the full inclusion of quantum statistics with that which arises from an identical calculation in which helium atoms are considered distinguishable (dotted line in Fig. 4) shows how the former extends to considerably longer distances, roughly up to three lattice constants.This can be attributed to the increasing particle delocalization associated with quantum-mechanical exchanges.and is consistent with the observation that Table 1.Energetics of solid 4 He computed at T = 1 K by Monte Carlo simulations based on the continuous-space WA.Shown are the total energy per particle e (in K), as well as the kinetic energy e k (in K), and the pressure in bars.The results are mainly for the Aziz I [22] pair potential, but for comparison purposes, some results are included of calculations based on the Aziz II [30] and Aziz III [33] pair potentials.The results presented in this table pertain to calculations for a system comprising N = 216 atoms.Statistical errors (in parentheses) are on the last digit.They are not explicitly shown for the pressure, as they are consistently of the order of half a bar.Shown within square brackets are the ground state estimates from Ref. [60] n three-particle ring exchanges are the most frequent.The enhancement of the n(r) at large distances corresponds to a strengthening of the momentum distribution at low momenta.It is worth noting that a similar manifestation of Bose statistics can be observed in the momentum distribution of other Bose fluids and crystals, e.g., in liquid parahydrogen near freezing [52].The results shown in Fig. 4 also support the conclusions reached in early Monte Carlo studies [46,47], namely that in a perfect crystal of 4 He the one-body density matrix decays exponentially at large distances, displaying no noticeable dependence on the temperature, i.e., there is no Bose-Einstein condensation in the T → 0 limit [53].Indeed, all credible scenarios of possible "supertransport" in the solid phase of 4 He involve extended defects such as dislocations [54][55][56]; there is overwhelming evidence that supersolidity is not underlain by point defects such as vacancies or interstistials [57,58].Otherwise, theoretical work carried out mainly over the past two decades indicates that the existence of a supersolid phase requires that the pairwise interaction among particles not feature the kind of repulsive core that characterizes instead helium and other rare gas solids [59].

Energetics
All of our results for the energetics of solid 4 He are compiled in Table 1.We provides estimates for the total energy per atom e (in K), as well as for the kinetic energy e k (also in K) and for the pressure, obtained with different versions of the Aziz pair potential, specifically the Aziz I, Aziz II and Aziz III.The first general remark, regarding the dependence of results on the pair potential utilized in the calculation, is that the differences are very small, much smaller than the current experimental uncertainties.In general, newer versions of the Aziz pair potential yield an energy per atom some ∼ 0.2 K lower than that afforded by the original (Aziz I) potential.On the other hand, the pressure estimates are very close, whereas those for the kinetic energy are virtually identical for the Aziz I and Aziz II pair potentials, while the Aziz III yields values that are between 0.1 and 0.2 K higher.There is otherwise no detectable difference between the results obtained with the various versions of the Aziz potential for quantities such as the one-body density matrix, or the pair correlation functions.Altogether, given the typical uncertainties affecting the experimental determination of, e.g., the kinetic energy, there scarcely seem to be compelling reasons to pick any of the various refinements of the Aziz I potential (the only likely outcome being that of rendering more complicated the comparison between different theoretical calculations).
As mentioned above, the inclusion of three-body terms in the potential energy of interaction among helium atoms, becomes necessary in order to bring the computed equation of state in agreement with experiment, at pressures significantly higher than those considered in this work (e.g., in the MPa range).On the other hand, the effect on the kinetic energy is rather small [37,38].

Comparison with recent calculations
Also shown in Table 1 are the most recent estimates [60] published from solid 4 He in the same range of density considered here, all based on the Aziz III potential.These results were obtained by means of Diffusion Monte Carlo simulations, i.e., they are strictly speaking ground state estimates.However, given the weak dependence on the temperature displayed by our results (shown in Fig. 2) below 1 K, the comparison of our results, which are at finite temperature, to those of Ref. [60] is quite appropriate.Within DMC, Bose statistics can be incorporated by a proper choice of trial wave function [61], out of which the exact ground state is (in principle) projected.
Taken individually (i.e., considering just one specific density) the energy estimates obtained in this work may be considered in reasonable agreement with the DMC ones, in that the differences may not be necessarily regarded as statistically significant.However, the comparison across all density values consistently shows estimates at finite (T = 1 K) temperature to be slightly lower than the (supposedly "exact") T = 0 results.This situation has been observed repeatedly during the past decades, for various Bose systems [62][63][64].One may go back to the early GFMC [65] studies of solid 4 He based on the Aziz I potential [66], yielding a value of the ground state per 4 He atom of −5.175K in the solid phase at n = 0.02934 Å −3 , almost 1 K higher than that obtained here (see Table 1).For a long time, those estimates were held as reliable, largely due to the absence of experimental data and of any alternative calculation.
Given that the DMC calculations carried out in Ref. [60] have comparable numbers of atoms to those carried out here, and in any case finite-size effects on the energy are relatively small, the most plausible explanation for the (small) discrepancy between their results and ours can be ascribed to intrinsic limitations of the DMC projection method.Fundamental reasons have been proposed as to why finite temperature methods are a superior option to ground state ones, when it comes to studying the ground state of Bose systems [19], due to various sources of bias from which finite temperature techniques are unaffected.Examples of such sources of bias are the need for a trial wave function, which affects the estimation of all the quantities, as well as the finite population of random walkers, whose effect can be rather large [62,67].
There are also small differences between the results for the kinetic energy per particle obtained in this work and those reported in Ref. [60].This is less noteworthy than the disagreement between total energies yielded by the two calculations, however, first of all because the statistical errors on the kinetic energy quoted in Ref. [60] are fairly large, secondly because it is well-known that DMC does not allow for unbiased, reliable estimations of thermodynamic averages of observables that do not commute with the energy.

Discussion
In this paper we present results of extensive numerical studies based on the continuousspace WA, of the low-temperature properties of the solid phase of 4 He.We focused on a pressure interval ranging from approximately 25 bars (near the T = 0 melting line) all the way to approximately 140 bars.The calculations made use of several versions of the most popular interparticle pair potential traditionally utilized in computer simulation studies of solid helium.We have compared the results of our simulations to available experimental data, mainly for the EOS and for the kinetic energy per atom.The agreement between theoretical results and existing experimental estimates is satisfactory; within the experimental uncertainties, it is not possible to identify a specific version of the Aziz pair potential that affords a closer reproduction of available experimental data.In general, however, a theoretical model only including pairwise interactions appears sufficient to obtain a reliable quantitative account of the energetics, structure, and dynamics of the helium crystal, at least at moderate pressures.

Figure 1 .
Figure 1.Time step extrapolation of the kinetic energy per particle e k (K) in hcp 4 He at temperature T = 1 K, for a system comprising N = 216 atoms.The density is n = 0.0286 Å −3 .The results are obtained with the Aziz I pair potential.The solid line is a quartic fit to the data.Arrow points to the value of τ for which convergence is observed.

Figure 2 .
Figure 2. Kinetic energy per 4 He atom e k computed as a function of temperature, for a hcp crystal of density 0.0286 (0.0312) Å −3 , left (right) panel.The simulated system comprises N = 216 atoms, and the interatomic potential utilized is the Aziz I (1979).Dashed lines represent fits to the data, based on a constant value.Experimental data shown on the right panel are from Ref. [16].Not shown in the left panel is the experimental result from Ref. [15] at T = 1.6 K and n = 0.0288 Å −3 , namely 24.25(30) K.

Figure 3 .
Figure 3. Spherically averaged pair correlation function g(r) in hcp solid 4 He at three different values of the density, at T = 1 K.The calculations whose results are shown here were performed on a system comprising N = 216 atoms.

Figure 4 .
Figure 4. Spherically averaged one-body density matrix n(r) in hcp solid 4 He at three different temperatures.The density is n=0.0286Å −3 , and the simulated system comprises N = 216 atoms.Dotted line represents the T = 1 K estimate for a 4 He crystal in which atoms are regarded as distinguishable ( boltzmannons).