Phase behaviour of the quantum Lennard-Jones solid

The Lennard-Jones potential is perhaps one of the most widely-used models for the interaction of uncharged particles, such as noble gas solids. The phase diagram of the classical LJ solid is known to exhibit transitions between hcp and fcc phases. However, the phase behaviour of the quantum Lennard-Jones solid remains unknown. Thermodynamic integration based on path integral molecular dynamics and lattice dynamics calculations are used to study the phase stability of the hcp and fcc Lennard-Jones solids. The hcp phase is shown to be stabilized by quantum effects in PIMD while fcc is shown to be favoured by lattice dynamics, which suggests a possible re-entrant low pressure hcp phase for highly quantum systems. Implications for the phase stability of noble gas solids are discussed. For parameters equating to Helium, the expansion due to zero-point vibrations is associated with quantum melting: neither crystal structure is stable at zero pressure.


I. INTRODUCTION
Since its inception in 1924, the Lennard-Jones (LJ) potential, 1 has remained the canonical model for short-ranged particle interactions. The properties of this potential are uniquely defined due to the presence of only two parameters, σ and ε, which set the length and energy scales respectively. Despite this simplicity, the LJ system shows remarkably rich phase behavior, including transitions between hexagonal close packed (hcp) and face-centered cubic (fcc) solids. [2][3][4][5][6] For the classical system hcp is the most stable phase at low temperature and pressure conditions, with fcc becoming preferred upon heating and/or compression. However, the free energy differences between hcp and fcc are small, and careful calculation of long-range interactions as well as both harmonic and anharmonic thermal effects are important. 7,8 The attractive 1/r 6 part of the potential describes van der Waals interactions, making it suitable for studying the noble gas elements. Indeed, at low pressures the heavy noble gas elements such as Ar, Xe and Kr 9,10 , as well as small molecules like methane, 11 are well described as classical particles interacting via the LJ potential. However, as the density of the system is increased or the mass of the particles is decreased, quantum effects become increasingly important. He, the lightest of the noble gas elements, is dominated by quantum effects. This has motivated investigations into the quantum LJ system, which is far more complicated than the classical system. For one thing the uniqueness of the phase diagram is lost: in the quantum system there is an additional free parameter,h 2 /2m, which corresponds to the "quantumness" of the particles.
The solid phase of the quantum LJ system has been investigated using a number of computational methods. [12][13][14][15][16] These include quasi-harmonic lattice dynamics (QHLD), classical molecular dynamics, and path-integral methods. QHLD 17 is exact in the low temperature limit, and captures quantum effects such as the zero-point energy which are inaccessible to classical methods. However, it cannot capture anharmonic effects associated with interacting phonons, something well known to play an important role in quantum crystals at moderate temperatures 18 . Path integral methods 19,20 do not suffer from this shortcoming, capable in principle of providing exact properties of quantum systems. These methods exploit the path integral formulation of quantum mechanics, and sample the configuration space of the quantum system using molecular dynamics (PIMD) or Monte Carlo. However, the downside of these methods is that the computational effort required to obtain accurate results grows exponentially with the quantumness of the system, to the point that highly quantum systems are intractable without resorting to severe approximations. Thus QHLA and path integral methods are complimentary, the former accurately describing the low temperature limit, and the latter being suitable for moderate temperatures.
The focus of previous studies has been on properties such as the thermal expansion and heat capacity of the quantum LJ solid. These studies all assumed fcc to be the stable solid phase. However, as mentioned above, the classical LJ solid exhibits regions of both fcc and hcp stability, and thus the expectation is that the phase diagram of the quantum LJ solid should also exhibit both these phases. Interestingly, an exploration of the solid phase diagram in this system has not yet been undertaken. Experimentally, while there is some uncertainty with regards to the question of hcp vs. fcc stability in Ne, Ar, Kr and Xe at low pressures, it appears that the fcc phase is at least metastable in these systems. 21 By contrast, He is exceptional in that the hcp is the observed phase (notwithstanding a small region of bcc stability near melting at low pressures). 18 It is not known why He readily takes the hcp structure. One might speculate that quantum effects somehow act to stabilise the hcp structure, and that the effect is more pronounced in He than the heavier noble gas elements due to its high quantumness. An investigation into the phase behaviour of the quantum LJ solid will shine light on this.
Here, we use PIMD in combination with QHLD to determine the phase diagram of the quantum LJ solid, focusing on low pressure regime. We determine the relative stability of hcp and fcc for a range of quantumnesses, ranging from the classical limit to that comparable to He. In Section II we describe the methodology which underpins our calculations. Then in Section III we describe our model for the LJ solid, including approximations utilised in our calculations. Moreover in this section we provide computational details regarding our calculations. Results of our PIMD calculations are presented in Section IV, followed by results of our QHLD calculations in Section V. In Section VI we reconcile the PIMD and QHLD results and discuss their implications for the noble gas elements. Finally, in Section VII we restate our main conclusions.

A. Path integral molecular dynamics
The path integral formalism exploits the isomorphism between a system of N quantum particles and a set of P interacting replicas of the system, each consisting of N classical particles. The exact quantum partition function Z = Tr[e −βĤ ] is mapped onto a classical one 22 Z P , such that In this equation m is the mass of the particles, β ≡ 1/(k B T ) is the inverse temperature,r i,j is the position vector of particle j in replica i, and is the Hamiltonian for the set of replicas. Note that Z P is only strictly equal to the quantum partition function Z in the limit P → ∞. However, in practice calculations are necessarily limited to finite P. Ideally P is large enough that results are indistinguishable from the P → ∞ limit. It can be seen from Eqn. 3 that each particle interacts with its corresponding particles in adjacent replicas via a harmonic potential with spring constant The strength of the inter-replica interactions is therefore directly proportional to both the mass of the particles and the temperature. Moreover, the condition r P+1 = r 1 in Eqn. 3 signifies that r P+1,j = r 1,j for all j. Therefore the replicas form a closed loop, and so the resulting system is often referred to as a ring polymer, with each replica representing a bead in the polymer chain. From Eqn. 3 it can also be seen that the particles additionally interact within each replica according to the given interatomic potential U (|r i,j − r i,l |), which in this case is the LJ potential.
The quantum system described by H P can be sampled using molecular dynamics techniques. To do so, conjugate momenta p i,j are added to the Hamiltonian such that

B. Thermodynamic integration in PIMD
To evaluate the relative stability of the fcc and hcp LJ solids, we require a comparison of the free energies of both phases. We use thermodynamic integration, a robust and widelyused technique for calculating free energies from MD simulations, 27 to obtain the free energies of both phases.
In thermodynamic integration the free energy difference between two states A and B is obtained by introducing a coupling parameter λ to the partition function Z of the system such that Z(λ = 0) = Z A and Z(λ = 1) = Z B . Since F = −k B T ln Z, the free energy difference can then be accessed as If A is a reference state with a known free energy then the above equation can be used to determine the free energy of state B. First ∆F BA is determined by integrating ∂ ln Z/∂λ over λ, and then this ∆F BA is added to F A to obtain F B . This approach is routinely used to obtain the free energy of a given classical crystal. In this case state A is chosen to be an Einstein crystal whose free energy is known analytically, while state B is chosen to be the true crystal. The methodology for performing such calculations is well-documented. [28][29][30] The generalization of thermodynamic integration to a quantum system is straightforward: we must include contributions to the free energy from nuclear quantum effects and so we add a second thermodynamic path between the quantum crystal (state B) and the classical crystal (state A). Thus the overall free energy F of the quantum crystal is separated into two terms: the classical free energy F c plus the excess quantum free energy ∆F q : Nuclear quantum effects in PIMD are controlled by the bead-bead interaction term and so ∆F q can be calculated via thermodynamic integration by tuning the strength of this interaction. The most straightforward choice is to use the particle mass µ as the coupling parameter and slowly vary it from the true atomic mass m 0 to infinite mass in the classical limit: Evaluation of this derivative yields where T prim µ is the primitive estimator for the quantum kinetic energy. This estimator has poor convergence properties with large P, however, and so without loss of generality it can be replaced with the more well-behaved centroid-virial estimator T vir , [31][32][33] where Finally, for ease of numerical integration a change of variables to g = m 0 /µ is done [34][35][36] to allow for integration in the range [0,1]: As this procedure is quite computationally expensive, it is not feasible to use mass thermodynamic integration to fully explore the phase diagram of the quantum LJ solid. Instead, we chose to perform a full calculation of the quantum-corrected free energy for a single reference point F 0 (V 0 , T 0 ) and then use Gibbs-Helmholtz integration of the free energy to generate the rest of the phase diagram. For a given (V, T ) point the free energy can be calculated from the reference point using the thermodynamic relationships for constant volume and for constant temperature. PIMD trajectories must still be run to obtain values for U (T ) and P (V ), but now only one trajectory is required for each point instead of a full mass thermodynamic integration.
C. Quasi-harmonic lattice dynamics QHLD 17 entails calculating the phonon density of states for a range of densities, and using these data, in conjunction with equations such as those given below, to calculate physical quantities such as the pressure and free energy. While this approach is less accurate than PIMD at moderate temperatures, it is insightful. In particular, it allows the hcpfcc free energy difference to be understood in terms of thermal and zero-point vibrational contributions. Such a decomposition is not possible in PIMD.
Consider a crystal phase at density ρ and temperature T . The Helmholtz free energy for the crystal can be decomposed as follows: where U GS is the ground state energy of the crystal and F vib is the vibrational contribution to the free energy. (Here, both U GS and F vib , are intensive quantities, and similarly for all other energies below). In the quasi-harmonic approximation the interatomic forces at the ground states for all ρ are assumed to be harmonic, with ρ-dependent force constants. In this case F vib can be expressed as follows: where N is the number of particles in the system, ω i is the angular frequency of the ith phonon (of which there are 3N − 3, excluding the translational modes) for the considered ρ, is the zero-point energy, and ω denotes the mean phonon frequency. In the classical limit, (up to an inconsequential temperature-dependent constant). Another important limit is Here the second term in Eqn. 16 vanishes, leaving F vib = F zp .
The free energy difference between the hcp and fcc phases at ρ, ∆F ≡ (F hcp − F fcc ) (and similarly for ∆U GS , ∆F vib , etc.), can be decomposed similarly to F above. Using the above equations it can be shown that with in the classical limit, and in the zero-temperature limit, where ω hcp(fcc) i denotes the angular frequency of the ith phonon for the hcp(fcc) crystal.
Note that, as can be seen from Eqn. 20, ∆F depends on the hcp and fcc phonon spectra only through ∆ ω in the zero-temperature limit. Here, quantum vibrational effects favour the structure with lowest zero-point energy, i.e. lowest ω . Similarly, it can be seen from Eqn. 20 that in the classical limit, vibrational effects act to stabilise the structure with the lowest mean log-frequency ∆ ln ω . Moreover, since the phonon frequencies ω i are proportional to 1/ √ m, ∆F is independent of the mass of the particles in the classical limit because changing m leaves ∆ ln ω unchanged. By contrast, for the general case (Eqn. 19) ∆F depends on the masses throughh/ √ m, which is a measure of the quantumness of the system.
The above equations can be used to determine F hcp , F fcc and ∆F at a given ρ and T from the hcp and fcc phonon densities of states at ρ. By considering many ρ and T , the regions of the ρ-T phase diagram where hcp is stable (∆F < 0) and where fcc is stable (∆F > 0) can be deduced. Moreover, it is also possible to use the hcp and fcc densities of states over a range of ρ to determine the P -T phase diagram. This is achieved by first calculating the hcp and fcc Gibbs free energies as functions of P and T via where ρ in this expression is the ρ such that P (ρ , T ) = P , and is the pressure at a given ρ and T . Then, the Gibbs free energy difference ∆G = (G hcp −G fcc ) is evaluated as a function of P and T . Finally, ∆G(P, T ) is used to deduce the hcp and fcc regions of the P -T phase diagram similarly to above for the ρ-T phase diagram: ∆G < 0 indicates that hcp stable; ∆G > 0 that fcc is stable.

A. Reduced units
We consider a system of N distinguishable quantum particles interacting via the LJ potential (Eqn. 1). As the length and energy scales are set by the σ and ε parameters, it is convenient to define all physical properties of the system in terms of dimensionless reduced units, as in Table I. The inclusion of quantum effects adds a second lengthscale to the system: the de Boer parameter 37 Λ * . This is a dimensionless quantity which describes the relationship between the particle diameter σ and the de Broglie wavelength of particles with energy ε. Large values of Λ * indicate a more delocalized quantum system, while Λ * = 0 corresponds to the classical limit. For the noble gases Λ * ranges from ≈ 0.01 (for Xe) to ≈ 0.4 (for He). 18 Thus we consider the range Λ * = 0 to 0.4 in this work.

B. Truncation scheme
In MD simulations it is necessary to truncate the potential interactions in order to avoid artefacts due to self-interaction through the periodic boundary. Previous work 3,7,8

Quantity Expression
Length has demonstrated that for the classical LJ solid one must take great care with regards to the truncation scheme and the treatment of the long-range interactions. The most commonlyused scheme is to shift the potential so that it is continuous at some cutoff radius r * c : This treatment, referred to as the spherically truncated and shifted (STS) model, avoids errors from discontinuous jumps in the potential but it fails to account for the interactions occurring beyond r * c . It therefore displays differing phase behaviour from the "true" LJ potential. 3 Conventional tail corrections, 27 which assume that the radial distribution function g(r * ) is uniform and equal to 1 at r * > r * c , are not useful for our purposes because they are independent of the crystal structure.
Instead, the contributions of these long-range interactions to the total energy of the system can be accounted for using what we will refer to here as the ground state perturbation (GSP) model. 3 Here the ground state of the LJ system, where all particles reside on their lattice sites, is treated exactly and only the excitations of the system are subject to truncation.
To do so, we introduce a correction term to the potential: where U LRC is defined by U GS (ρ * ) is the ground state energy of the untruncated "true" LJ system at density ρ * and R ij is the inter-particle distance in said ground state. The U GS (ρ * ) term is found from lattice-summation 2,38-40 as The A 12 and A 6 terms have been tabulated for different phases of the LJ solid, and in this work parameters for the fcc and hcp phases were taken from ref. 41. For simulations in the N V T ensemble where density is constant, the U LRC term simply amounts to a constant shift in the relative fcc and hcp energies and will therefore not affect the dynamics of either system.

C. Computational details
Classical and PIMD simulations were performed with LAMMPS 42 using the i-PI wrapper. 43 Simulations were run with σ = 2.96 Å and ε = 0.00295 eV, parameters which have been shown to give good results for PIMC calculations of noble gas solids. 15,16 Systems of 256 Lennard-Jones particles with either the fcc or hcp crystal structures were initialized at a specified density ρ * . Trajectories were then initiated in the NVT ensemble with orthorhombic periodic boundary conditions. Temperature was kept constant using the stochastic PILE-G thermostat 24 with a relaxation time of 0.01 t * . Simulations were run at temperatures ranging from T * = 0.10 to 0. 50

A. Classical free energies
As a starting point for the PIMD calculations, we first determined the classical free energy term F * c at the chosen reference point of T * 0 = 0.10 and ρ * 0 = 1.07255 using the standard Einstein crystal method [28][29][30] . The interaction cutoff r * c was initially chosen to be 2.5, as it is the most popular choice in the literature. The sensitivity of the F * c term to the choice of cutoff length r * c was investigated by running simulations with r * c = 2.2 and 2.8 as well. The results are listed in Table II, with the phase stability expressed relative to the hcp phase as Predictably, the phase stability of the STS model is highly dependent on the chosen cutoff, 7,8 with r * c = 2.2 and 2.5 favouring hcp while 2.8 favours fcc. To reduce this cutoff effect the GSP correction to ∆F * c was calculated from the ground state fcc and hcp structures using Eqn. 26 as ∆U (ρ . The resulting correction term and corrected free energy difference ∆F * c,corr = ∆F * c + ∆U (ρ * ) LRC is also listed in Table II The resulting phase diagrams are shown in Figure 1 for both the STS and GSP models.
The STS model shows that hcp is stable at low densities and is in agreement with previous assessments of the phase behaviour of the uncorrected r * c = 2.5 LJ solid. 41 Inclusion of longrange GSP correction shifts the phase boundary by stabilizing fcc over hcp, and as such the hcp phase is only stable at low density and temperature. Since the fcc phase has been shown to have higher entropy than the hcp phase, 47-51 it is then perhaps unsurprising that inclusion of long-range order works to stabilize fcc. All in all, these results demonstrate the delicate balance between solid phases and the surprising complexity in this simple model.
Further information about the phase stability may be gleaned from a consideration of finite size effects, but due to the large computational cost for PIMD simulations a study of large systems is beyond the scope of this work.

B. PIMD convergence
The number of beads used in a PIMD simulation is a very important choice; P must be large enough to accurately probe the quantum limit, but small enough that the computational cost is still affordable. The number of beads required depends on the relative strength of the quantum harmonic energy levels versus the thermal energy. A typical rule of thumb given for the minimum number of beads is P min = 4βhω max , 20 where ω max is the highest vibrational frequency in the system. However, in practice the required number of beads is often much higher than this minimum limit. To choose an appropriate number for the quantum LJ solid, PIMD simulations were run with increasing P and the system was deemed to be converged when the average internal energy E * = U * LJ + T * vir was within 0.34 ε, which corresponds to ≈ 1 meV/atom in real units. The convergence of the structural properties of the system was also monitored via the radial distribution function of the ring polymer beads. Convergence was reached at P = 144 for T * = 0.10. Convergence plots are shown in Figure 2 for both energetic and structural properties.

C. Quantum free energies
The excess quantum free energies ∆F * q for each phase were obtained from PIMD using a 13-point mass thermodynamic integration from m = m 0 to m = ∞ at the chosen reference point T * 0 = 0.10, ρ * 0 = 1.07255. The free energy difference between the two phases was monitored and trajectories 150-200 t * in length were required for adequate convergence of this quantity. An example is shown in Figure 3. The phase stability of the quantum LJ solid was then calculated using Eqns. 13 and 14 with F * 0 = F * c + ∆F * q at the reference point of T * 0 = 0.10, ρ * 0 = 1.07255. For each Λ * value, calculations were run from ρ * = 1.30 to whichever density gave P * ≈ 0. Note that only densities where the system remained solid across the whole temperature range were considered here, thus in some cases it was not possible to reach P * ≈ 0 due to melting at higher temperature. The resulting phase diagrams are shown in Figure 4 for both the STS and GSP models. Phase stability is represented as ∆F * = F * hcp − F * fcc . The quantum contribution to the pressure is immediately apparent in these plots, since the highly quantum systems require significantly lower densities to reach P * ≈ 0 than in the classical case. Even though the GSP model favours fcc more than the STS model, the inclusion of nuclear quantum effects stabilizes hcp over fcc in both cases. This is evidenced by the increase in size of the hcp region with increasing quantumness, which is observed even with a relatively heavy particle at Λ * = 0.1. Since these phase diagrams are plotted against ρ * , this growth of the hcp region means that the fcc region is being pushed to higher and higher pressures with increasing quantumness and is therefore being destabilized relative to hcp. Interestingly, the slope of the phase boundary changes sign at Λ * = 0.4. This may be due to fluctuations in this highly quantum system, or it may indicate a change in the nature of the phase boundary at high quantumness.

A. Initial investigations
The hcp and fcc phonon density of states (DoS) for r * c = 10.0 at ρ * = 1.07255, the density for the classical LJ solid at T * = P * = 0, are shown in Figure 5. The figure shows that the hcp DoS has both low and high frequency peaks, while fcc has more intermediate frequency modes. Moreover, both structures share the same peak at high frequency, though the highfrequency peak is larger in fcc than hcp. However, despite the qualitatively different shapes of the densities of states of the two crystals, the mean frequencies and mean log-frequencies, which, as discussed in Section II C, play an important role in determining which of the structures is stable, are indistinguishable on the scale of this figure. Hence the fine structure of the DoS must be considered to determine which of hcp and fcc is stable. We return to this point in a moment. Preliminary calculations revealed that the hcp and fcc crystals were mechanically unstable for densities less than ρ * = 0.8: at these densities the system exhibited phonons with imaginary frequencies for all cutoffs considered. Hence, keeping in mind that we are interested in the low pressure region of the phase diagram, we focused on densities ranging from ρ * = 0.8 to 1.3.
The hcp pressure is shown as a function of T * over this density range for various Λ * in Figure 6. The fcc pressure is indistinguishable from that of hcp on the scale of this figure.
As expected, increasing the quantumness while fixing the density results in an increase in the pressure of the system. This is primarily due the zero-point vibrations. To elaborate, from Eqns. 17 and 23 the contribution to the pressure from this energy is Noting that ∂ hω /∂ρ is positive and proportional toh/ √ m, and hence also Λ * , it follows that P zp is also proportional to Λ * . Figure 6 also reveals that for Λ * = 0.2, 0.3 and 0.4 there is no mechanically stable hcp or fcc density corresponding to P = 0 in the quasi-harmonic approximation: the figure implies that the hcp and fcc densities for P = 0 would be achieved at ρ < 0.8, which, as just mentioned, are mechanically unstable within the approximation.
Since the quasi-harmonic approximation is valid in the low temperature limit, for Λ * ≥ 0.2 some phase other than hcp or fcc must therefore be stable at P = T = 0. To validate our implementation of QHLD we also considered Λ * = 0.0103, 0.0166, 0.0296 and 0.0896, which correspond to Xe, Kr, Ar and Ne respectively, 14 and compared results for r * c = 10 to those of ref. 14. Plots of P * vs. ρ * for these Λ * at T * = 0 (not shown) were in good agreement with those given in ref. 14, as were plots of ρ * vs. T * at P * = 0 (not shown).

B. Sensitivity to cutoff
To investigate the sensitivity of the phase behaviour to r * c , we focused on ρ * = 1.07255. We considered r * c = 2.2, 2.5 and 2.8, supplementing our calculations described above for this density using r * c = 10. The hcp and fcc DoS for r * c = 2.2, 2.5 and 2.8 are almost identical to those shown in Figure 5 (which recall are for r * c = 10) on the scale of the figure. However, differences in the fine structure of the DoS for different r * c have important implications for the stability of hcp vs. fcc.
Recall that in the zero-temperature limit ∆F depends on the zero-point energy through the difference in the mean phonon frequencies ∆ ω ≡ ( ω hcp − ω fcc ) (Eqn. 21), while in the classical limit ∆F depends on the difference in the mean log-frequencies ∆ ln ω (Eqn. 20). In Table III ∆ ω * and ∆ ln ω * are compared for various values of r * c . It can be seen that ∆ ln ω * is positive for all r * c . This implies that vibrational effects act to stabilise fcc in the classical limit for all considered r * c . This follows from Eqn. 20: if ∆ ln ω * > 0, ∆F * increases with T , which corresponds to fcc stabilisation.
By contrast, the sign of ∆ ω * depends on r * c . ∆ ω * is negative at r * c = 2.2, implying that the zero-point energy of hcp is lower than that of fcc. Since increasing quantumness increases the size of the zero-point contribution to ∆F * (c.f. Eqn. 20), this means that increasing quantumness stabilises hcp for r * c = 2.2 (in the zero-temperature limit). On the other hand, for r * c = 2.5, ∆ ω * is positive, which implies the opposite, i.e. that fcc is stabilised by quantum effects. The same is true for r * c = 2.8 and r * c = 10.0, though for these cutoffs ∆ ω * is of a smaller magnitude and hence the stabilisation of fcc by quantum effects is less pronounced than at r * c = 2.5. The above discussion is borne out in Figure 7, which shows ∆F * vs. Λ * at this density for various T * and r * c . Note that for T * = 0 increasing the quantumness increases ∆F * for r * c = 2.5, 2.8 and 10.0, with a more pronounced increase for r * c = 2.5; while increasing quantumness decreases ∆F * for r * c = 2.2. These trends are in accordance with the values of ∆ ω * for each cutoff described above. Note also that, at Λ * = 0, increasing T * increases ∆F * , in accordance with the above discussion that ∆ ln ω * is positive for all r * c , i.e. that vibrational effects always stabilise fcc in the classical limit. This is also the case away from the classical limit: thermal effects always act to stabilise fcc.
A key result to be drawn from the above discussion is that the difference in the zero-point TABLE III. Differences in the phonon mean frequencies and mean log-frequencies between the hcp and fcc structures, i.e. ∆ ω * and ∆ ln ω * , for the LJ solid at ρ * = 1.07255 and various cutoffs r * c .
A positive value indicates that ∆ ω * or ∆ ln ω * is higher for the hcp structure than for the fcc structure. energies of the hcp and fcc structures in the LJ solid is highly sensitive to r * c , to the point that changing r * c can reverse whether hcp or fcc is stablised by increasing the quantumness of the system. There is no analogous problem in the classical case; the discussion above, and previous studies, have revealed that changing r * c does not change the fact that the fcc structure is stabilised by thermal effects.

C. Phase behaviour
We now consider the phase behaviour for r * c = 10 over the density range ρ * = 0.8 to 1.3. ∆F * vs. ρ * is shown in Figure 8 for various Λ * . The figure shows that increasing quantumness generally has the effect of stabilising fcc: as Λ * is increased, ∆F * also increases. The effect on the phase diagram is shown in Figure 9, which shows the hcp-fcc phase boundaries for various Λ * in the ρ * -T * plane. Note that increasing quantumness has the effect of reducing the size of the hcp region; the hcp region of the phase diagram is "compressed" towards ρ * = 0.8 as Λ * is increased.
This behaviour is largely due to the zero-point energy. Recall that for r * c = 10 we found for ρ * = 1.07255 that the zero-point energy favours fcc: ∆E * ZP > 0. The same is true at all other densities we considered. As the quantumness is increased, ∆E * ZP becomes a larger contribution to ∆F * , and therefore fcc becomes increasingly favoured. This is especially true at high densities, where, as shown in Figure 8, ∆E * ZP is larger than at lower densities, leading to a more pronounced effect. Of course, temperature also affects ∆F * through the difference in the vibrational energies.
At low pressures ∆E * ZP is smaller and temperature plays a larger role. At very low densities the phase behaviour is non-trivial: at ρ * = 0.84 as Λ * is increased there is an increase in the hcp-fcc transition temperature until Λ * = 0.2, followed by a decrease as Λ * is increased further.
As can be seen from Figure 6, much of the range covered in Figures 8 and 9 pertain to densities which correspond to negative pressures. ∆G * vs. P * is shown in Figure 10, and the hcp-fcc phase boundaries in the P * -T * plane are shown in Figure 11. Recall that neither crystal structure is mechanically stable at P * = 0 for Λ * > ∼ 0.2 in the quasi-harmonic approximation, which is why not all the curves in these figures extend to P * = 0. Figure 11 reveals that as the quantumness is increased the location of the hcp-fcc transition is moved to lower pressures. Moreover, at P * = 0, increasing Λ * from 0 to 0.1 moves the P * = 0 transition and 9). Finally, the hcp-fcc free energy differences obtained from PIMD are typically about an order of magnitude larger than those obtained from QHLD (see Figures 4 and 9).
There are a number of possible causes for the discrepancies between the two methods.
Firstly, the models used in the QHLD and PIMD calculations were different: the QHLD calculations employed r * c = 10, while the PIMD calculations used a value of r * c which scaled commensurately with the density (Eqn. 29). To investigate this further, we repeated our QHLD calculations using the same cutoff scheme used in our PIMD calculations. However, we found that this worsened the agreement between the two methods: using the scaled cutoff scheme makes the fcc more stable in the quasi-harmonic approximation than it is for r * c = 10. This can be seen by comparing Figures 8 and 12, the latter of which shows ∆F vs. ρ * using the scaled cutoff scheme. Another key difference between our QHLD and PIMD calculations is that they pertain to different system sizes, namely N = ∞ and N = 256 particles, respectively. It is well known that the quantitative details of the phase diagram for the classical LJ solid are sensitive to the system size, and surely the same is true for the quantum solid. However, we do not believe that finite size effects are the main cause of the discrepancies. Rather, the main cause is the approximations which underpin the two methods.
To elaborate, QHLD is functionally exact in the zero-temperature limit but will break down at some finite temperature, and it is well known that this break down occurs very quickly in quantum molecular crystals 18 . PIMD, on the other hand, is more accurate than QHLD at finite temperature but calculations cannot be performed at temperatures near zero due to the prohibitive computational cost. The two methods therefore provide complimentary data: QHLD informs us about the zero-temperature limit, and PIMD speaks to the high-temperature limit. The true phase diagram of the quantum system may be a combination of the results from these two methods, where at low temperatures fcc is increasingly favoured as the quantumness is increased, and at high temperatures hcp is increasingly favoured. This suggests a potential re-entrant fcc-hcp-fcc transition at low pressures for a highly-quantum LJ system, as illustrated schematically in Figure 13.

B. Origin of hcp stabilisation in PIMD
The origin of the fcc stabilisation with increased quantumness in QHLD was shown to be due to its lower mean vibrational frequency and thus lower zero point energy. The origin of hcp stabilization in PIMD is less clear, due to the opaque nature of PIMD simulations. To investigate this further we partition the ∆F * values from Figure 4 into contributions from internal energy and entropy. The average difference in internal energy ∆E * can be extracted directly from the PIMD trajectories and so −T * ∆S * is easily accessible as: The results are shown in Figure 14 for the cases of Λ * = 0, 0.1 and 0.3 at the reference temperature T * = 0.1. Here we only consider the GSP-corrected model, as we do not wish to include fluctuations in the energy due to cutoff effects. The ∆F * curve is reminiscent of the QHLD result in Figures 8 and 12, but larger in magnitude and shifted to more dramatically favour hcp at low densities. From these plots we can clearly see that ∆E * , which contains the zero point energy, is practically negligible and that the free energy difference between the two phases is due to the difference in entropy −T * ∆S * . Therefore the stabilization of quantum hcp in PIMD must be due to anharmonic effects which are not accounted for in QHLD.

C. Implications for noble gases
The LJ potential is often used to model interactions in noble gas solids and it is therefore interesting to consider our results in that context. Of the noble gases, only He has been observed in the hcp phase; all others adopt the fcc structure. 52,53 Explanations for this phenomenon have been offered in the literature 54-56 but a consensus has not been reached.
In Figure 15 we compare the values of the de Boer parameter Λ * for the noble gases with our PIMD results for the GSP-corrected model. By and large these results agree with the experimental observations. The heavy noble gases (Ar, Kr and Xe) are clustered near the fcc phase, while He has a strong preference for hcp. The only discrepancy is Ne, which is shown to prefer hcp with this treatment but has an fcc structure in reality. However, the LJ potential is a somewhat simplistic treatment and a more realistic potential may show more accurate phase behaviour for Ne. Overall, our results show that the hcp structure of He is due to quantum effects, while the heavy noble gases prefer the classically-favoured fcc structure. Moreover, there are very large contributions of zero-point vibrations to the pressure, in the case of helium these are large enough to drive the density below the region in which the Lennard-Jones potential can stabilize a crystal structure.

D. Prospect of other phases
One key conclusion of our QHLD calculations was that P = 0 was inaccessible for Λ > ∼ 0.3 on account of the mechanical instability of the fcc and hcp crystals at very low pressures.
Keeping in mind that QHLD is a perturbation method, the question remains as to whether any crystal structure is stable for Λ > ∼ 0.3 at such pressures. In this work we have only considered the hcp and fcc phases, which are the only two phases known to be stable in the classical LJ solid. Of course, there is the prospect that phases other than hcp and fcc are stable at certain temperatures and low pressures for the quantum LJ solid.

VII. CONCLUSION
Through a combination of PIMD simulations and lattice dynamics calculations, the inclusion of nuclear quantum effects in the LJ solid has been shown to stabilize the fcc phase at low temperature and the hcp phase at high temperature. In addition, the quantum effects on the phase behaviour of the quantum LJ solid in PIMD is relatively insensitive to the truncation scheme or cutoff length used, which makes it somewhat easier to draw definitive conclusions than their classical counterparts -despite the increased computational cost. We also offer an explanation for the experimentally observed phase behaviour of the noble gas solids. Helium, which is highly quantum, has a sufficiently large de Boer parameter that it falls squarely in the quantum-favoured hcp region. Thus the contrast between hcp helium and fcc structures for the other noble gases solids is purely due to quantum effects.