Muonium hydride: the lowest density crystal

A muonium hydride molecule is a bound state of muonium and hydrogen atoms. It has half the mass of a parahydrogen molecule and very similar electronic properties in its ground state. The phase diagram of an assembly of such particles is investigated by first principle quantum simulations. In the bulk limit, the low-temperature equilibrium phase is a crystal of extraordinarily low density, lower than that of any other known atomic or molecular crystal. Despite the low density and particle mass, the melting temperature is surprisingly high (close to 9 K). No (metastable) supersolid phase is observed. We investigated the physical properties of nanoscale clusters (up to 200 particles) of muonium hydride and found the superfluid response to be greatly enhanced compared to that of parahydrogen clusters. The possible experimental realization of these systems is discussed.


I. INTRODUCTION
An intriguing open question in condensed matter physics, one with potential practical significance, is whether there exists a lower bound for the density of a crystal. In other words, how large can the average distance between nearest neighbor atoms be, before crystalline long-range order is suppressed by thermal and quantum fluctuations? It remains unclear whether a lower limit exists, or how to approach it experimentally. We begin our discussion by rehashing a few basic facts about the crystalline phase of matter.
Crystallization occurs at low temperature (T ) in almost all known substances, as the state of lowest energy (ground state) is approached. Classically, the ground state is one in which the potential energy of interaction among the constituent particles is minimized, a condition that corresponds to an orderly arrangement of particles in regular, periodic lattices. In most cases, quantum mechanics affects this fundamental conclusion only quantitatively, as zero-point motion of particles results in lower equilibrium densities and melting temperatures, with respect to what one would predict classically; typically, these corrections are relatively small. Only in helium is the classical picture upended by quantum mechanics, as the fluid resists crystallization all the way down to temperature T = 0,K, under the pressure of its own vapor.
In many condensed matter systems, the interaction among atoms is dominated by pairwise contributions, whose ubiquitous features are a strong repulsion at interparticle distances less than a characteristic length σ, as well as an attractive part, which can be described by an effective energy well depth . The dimensionless param- * m.boninsegni@ualberta.ca eter Λ = 2 /[m σ 2 ], where m is the particle mass, quantifies the relative importance of quantum-mechanical effects in the ground state of the system.
It has been established for a specific model pair potential incorporating the above basic features, namely the Lennard-Jones potential, and for Bose statistics, that the equilibrium phase is a crystal if Λ < Λ c ≈ 0.15; it is a (super)fluid if 0. 15 Λ 0.46, while, for Λ > 0.46, the system only exists in the gas phase [1]. The value of Λ c is estimated [2] to be slightly higher (≈ 20 %) for systems obeying Fermi statistics. For substances whose elementary constituents are relatively simple atoms or molecules, this result provides a useful, general criterion to assess their propensity to crystallize, as measured by their proximity in parameter space to a quantum phase transition to a fluid phase. It can be used to infer, at least semi-quantitatively, macroscopic properties of the crystal, such as its density and melting temperature, both of which generally decrease [1] as Λ → Λ c .
The two stable isotopes of helium have the highest known values of Λ, namely 0.24 (0.18) for 3 He ( 4 He); for all other naturally occurring substances, Λ is considerably lower. The next highest value is that of molecular hydrogen (H 2 ), namely Λ = 0.08, quickly decreasing for heavier elements and compounds. At low temperature, H 2 forms one of the least dense crystals known, of density ρ = 0.0261Å −3 (mass density 0.086 gr/cc). The low mass of a H 2 molecule (half of that of a 4 He atom) and its bosonic character (its total spin is S = 0 in its ground state) led to the speculation that liquid H 2 may turn superfluid at low temperature [3], just like 4 He. In practice, no superfluid phase is observed (not even a metastable one), as molecular interactions, and specifically [4] the relatively large value of σ (∼ 3Å) cause H 2 to crystallize at a temperature T = 13.8 K.
Present consensus is that H 2 is a non-superfluid, insulating crystal at low temperature, including in reduced dimensions [5,6]; only small clusters of parahydrogen (∼ 30 molecules or less) are predicted [7][8][9][10] to turn superfluid at T ∼ 1 K, for which some experimental evidence has been reported [11]. If, hypothetically, the mass of the molecules could be progressively reduced, while leaving the interaction unchanged, thus increasing Λ from its value for H 2 all the way to Λ c , several intriguing scenarios might arise, including a low temperature superfluid liquid phase, freezing into a low-density crystal at T = 0, and even a supersolid phase, namely one enjoying at the same time crystalline order and superfluidity [12]. Obviously, the value of Λ can also be modified by changing one or both of the interaction parameters ( and σ) independently of the mass; in this work, however, we focus for simplicity on the effect of mass on the physics of the system.
One potential way to tune the mass is via substitution of protons or electrons with muons. For example, an assembly of molecules of muonium hydride (HMu) differ from H 2 by the replacement of one of the two protons with an antimuon (µ + ), whose mass is approximately 11% of that of a proton. Quantum chemistry calculations have shown that it is very similar in size to H 2 , and has the same quantum numbers in the ground state [13]. It is therefore not inconceivable that the interaction between two HMu might be quantitatively close to that between two H 2 (we further discuss this aspect below). In this case the value of the parameter Λ is ∼ 0.14, i.e., very close to Λ c for a Bose system. This leads to the speculation that this substance may crystallize into a highly quantal solid, displaying a strikingly unique behavior, compared to ordinary crystals. In order to investigate such a scenario, and more specifically to gain insight into the effect of mass reduction, we studied theoretically the low temperature phase diagram of a hypothetical condensed phase of HMu, by means of first principle computer simulations of a microscopic model derived from that of H 2 .
The main result is that the equilibrium phase of the system at low temperature is a crystal of very low density, some ∼ 5% lower than the T = 0 equilibrium density of liquid 4 He. Despite the low density, however, such a crystal melts at a fairly high temperature, close to 9 K, i.e., only a few K lower than that of H 2 . No superfluid phases are observed, either fluid or crystalline, as exchanges of indistinguishable particles, known to underlie superfluidity [14], are strongly suppressed in this system, much as they are in H 2 , by the relatively large size of the hard core repulsion at short distances (i.e., σ). As a result, the behavior of this hypothetical system can be largely understood along classical lines. This underscores once again the crucial role of exchanges of identical particles in destabilizing the classical picture, which is not qualitatively altered by zero-point motion of particles alone [15]; it also reinforces the conclusion [4] reached elsewhere that it is the size of the hard core diameter of the intermolecular interaction that prevents superfluidity in H 2 , not the depth of its attractive part.
To complement our investigation, we also studied nanoscale HMu clusters of varying sizes, comprising up to few hundred molecules. We find the behavior of these clusters to be much closer to that of 4 He (rather than H 2 ) clusters. For example, at T = 1 K the structure of HMu clusters is liquid-like, and their superfluid response approaches 100%, even for the largest cluster considered (200 HMu molecules). Thus, while mass reduction does not bring about substantial physical differences between the behavior of bulk HMu and that of H 2 , it significantly differentiate the physics of nanoscale size clusters.

II. MODEL
We model the HMu molecules as point-like, identical particles of mass m and spin S = 0, thus obeying Bose statistics. For the bulk studies, the system is enclosed in a cubic cell of volume V = L 3 with periodic boundary conditions in the three directions, giving a density of ρ = N/V . For the cluster studies, the N particles are placed in a supercell of large enough size to remove boundary effects. The quantum-mechanical many-body Hamiltonian reads as follows: where the first (second) sum runs over all particles (pairs of particles), λ ≡ 2 /2m = 21.63 KÅ 2 (reflecting the replacement of a proton with a µ + in a H 2 molecule), r ij ≡ |r i − r j | and v(r) denotes the pairwise interaction between two HMu molecules, which is assumed spherically symmetric.
FIG. 1. Geometry utilized in the calculation of the effective interaction between two HMu molecules. Shown is the line connecting the centers of mass of the two molecules, as well the three angles describing their relative orientation, i.e., the two polar angles θ1, θ2, as well as the azimuthal angle φ.
In order to decide on an adequate model potential to adopt in our calculation, we use as our starting point the H 2 intermolecular potential, for which a considerable amount of theoretical work has been carried out [16][17][18]. We consider here for definiteness the ab initio pair potential proposed in Ref. 18. The most important effect of the replacement of one of the protons of the H 2 molecule with a µ + is the shift of the center of mass of the molecule away from the midpoint, and toward the proton. We assume that this effect provides the leading order correction with respect to the H 2 intermolecular potential, and we set out to obtain a corrected version of the interaction for the new geometry, as illustrated in Fig. 1. We use the program provided in Ref. [18] to generate a potential as a function of the distance between the midpoints and the angular configurations, and then transform that to a potential as a function of the distance between the centers of mass and the angular configurations. Finally, we average over the angular configurations to obtain a onedimensional isotropic potential. We take the distance between the proton and the µ + to be that computed in Ref. [19], which differs only slightly from the distance between the two protons in the H 2 molecule.
In Fig. 2, we compare the potential energy of interaction between two HMu molecules resulting from this procedure with that of obtained for two H 2 molecules, i.e., with the center of mass at the midpoint. The comparison suggests that the displacement of the center of mass of the molecule results in a slight stiffening of the potential at short distance. Indeed, in the range of average interparticle separations explored in this work (namely the 3.5-3.7Å range), the differences between the interactions are minimal. Also shown in Fig. 2 is the Silvera-Goldman potential [16], which is arguably the most widely adopted in theoretical studies of the condensed phase of H 2 , and has proven to afford a quantitatively accurate [20] description of structure and energetics of the crystal. As one can see, it has a significantly smaller diameter and is considerably "softer" than the potential of Ref. 18; the reason is that it incorporates, in an effective way, non-additive contributions (chiefly triplets), whose overall effect is to soften the repulsive core of the pair-wise part computed ab initio.
It therefore seems reasonable to utilize the Silvera-Goldman potential [16] to carry out our study, as the use of a quantitatively more accurate potential is not likely to affect the conclusions of our study in a significant way. Furthermore, because the majority of the theoretical studies of condensed H 2 have been carried out using the Silvera-Goldman potential, its use in this study allows us to assess the effect of mass alone.

METHODOLOGY
The low temperature phase diagram of the thermodynamic system described by Eq. (1) as a function of FIG. 2. The intermolecular interaction energy E (in K) as a function of distance (Å) between the centers of mass of the molecules, for the H2 (red, dashed) and HMu (black, solid line) cases, obtained with the aid of the programs provided in Ref. [18]. Also shown (blue, solid line) is the Silvera-Goldman potential.
density and temperature has been studied in this work by means of first principles numerical simulations, based on the continuous-space Worm Algorithm [21,22]. Since this technique is by now fairly well-established, and extensively described in the literature, we shall not review it here; we used a variant of the algorithm in which the number of particles N is fixed [8,9]. Details of the simulation are standard; we made use of the fourth-order approximation for the short imaginary time (τ ) propagator (see, for instance, Ref. 23), and all of the results presented here are extrapolated to the τ → 0 limit. We generally found numerical estimates for structural and energetic properties of interest here, obtained with a value of the time step τ ∼ 3.0 × 10 −3 K −1 to be indistinguishable from the extrapolated ones, within the statistical uncertainties of the calculation. We carried out simulations of systems typically comprising a number N of particles. In the cluster calculations, N can be chosen arbitrarily. For the bulk, the precise value of N depends on the type of crystalline structure assumed; typically, N was set to 128 for simulations of body-centered cubic (bcc) and face-centered cubic (fcc) structures, 216 for hexagonal close-packed (hcp). However, we also performed a few targeted simulations with twice as many particles, in order to gauge the quantitative importance of finite size effects.
Physical quantities of interest for the bulk calculations include the energy per particle and pressure as a function of density and temperature, i.e., the thermodynamic equation of state in the low temperature limit. We estimated the contribution to the energy and the pressure arising from pairs of particles at distances greater than the largest distance allowed by the size of the simulation cell (i.e., L/2), by approximating the pair correlation function g(r) with 1, for r > L/2. We have also computed the pair correlation function and the related static structure factor, in order to assess the presence of crystalline order, which can also be detected through visual inspection of the imaginary-time paths.
We probed for possible superfluid order through the direct calculation of the superfluid fraction using the wellestablished winding number estimator [24]. In order to assess the propensity of the system to develop a superfluid response, and its proximity to a superfluid transition, we also rely on a more indirect criterion, namely we monitor the frequency of cycles of permutations of identical particles involving a significant fraction of the particles in the system. While there is no quantitative connection between permutation cycles and the superfluid fraction [25], a global superfluid phase requires exchanges of macroscopic numbers of particles (see, for instance, Ref. 14).

A. Bulk
We simulated crystalline phases of HMu assuming different structures, namely bcc, fcc and hcp. At low temperature, all of these crystals remain stable in the simulation. The energy difference between different crystal structures is typically small, of the order of the statistical errors of our calculation, i.e., few tenths of a K. This is similar to what is observed in H 2 [20]). Consequently, we did not attempt to establish what the actual equilibrium structure is, as this is not central to our investigation.
As a general remark, we note that the estimates of the physical quantities in which we are interested remain unchanged below a temperature T = 2 K. Thus, results for temperatures lower than 2 K can be considered ground state estimates.
We begin by discussing the equation of state of the system in the T → 0 limit. Fig. 3 shows the energy per HMu molecule computed as a function of density for a temperature T = 1 K. The solid line is a quadratic fit to the data, whose fitting parameters yield the equilibrium density ρ 0 = 0.02090(5)Å −3 , corresponding to an average intermolecular distance 3.63Å, as well as the ground state energy e 0 = −45.75(2) K. The ground state energy is almost exactly [20] one half of that of H 2 , with a kinetic energy of 70.4(1) K, virtually identical to that of H 2 at equilibrium, i.e., at a 30% higher density.
The equilibrium density ρ 0 is lower than that of superfluid 4 He, which is 0.02183Å −3 . The system is in the crystalline phase, however, as we can ascertain through the calculation of the static structure factor S(q), shown in Fig. 4. The result shown in the figure pertains to a case in which particles are initially arranged into a hcp crystal. The sharp peaks occurring at values of q corresponding to wave vectors in the reciprocal lattice signals long-range crystalline order.
It is interesting to compare the structure of a HMu crystal to that of two other reference quantum solids, namely H 2 (more precisely, para-H 2 ) at its equilibrium density, namely ρ = 0.0261Å −3 , and solid 4 He near melting, at density ρ = 0.0287Å −3 . The pair correlation functions shown in Fig. 5 were computed for hcp crystals at temperature T = 1 K (the results for H 2 and 4 He were obtained in Ref. 26). While the fact that the peaks appear at different distances reflects the difference in density, 4 He being the most and HMu the least dense crystals, the most noticeable feature is that the height of the peaks is much more pronounced in H 2 than in HMu and 4 He, whose peaks height and widths are similar. Fig. 6 shows the pressure of a HMu crystal computed at T = 1 K, as a function of density. In the relatively narrow density range considered, a linear fit is satisfactory, and allows us to compute the speed of sound v = (mρκ) −1/2 , where κ = ρ −1 (∂ρ/∂P ) is the compressibility. We obtained the speed of sound 1300 ± 100 m/s at the equilibrium density. Next, we discuss the possible superfluid properties of the crystal, as well as of the fluid into which the crystal melts upon raising the temperature. There is no evidence that the crystalline phase of HMu may display a finite superfluid response in the T → 0 limit. Indeed, the observation is that the behavior of this crystal is virtu-ally identical to that of solid H 2 , as far as superfluidity is concerned. The main observation is that exchanges of identical particles, which underlie superfluidity, are strongly suppressed in HMu, much like they are in H 2 , mainly due to the relatively large diameter of the repulsive core of the pairwise interaction. Indeed, exchanges are essentially non-existent in solid HMu at the lowest temperature considered here (T = 1 K), much like in solid H 2 ; the reduction of particle mass by a factor two does not quantitatively alter the physics of the system in this regard.
As temperature is raised, we estimate the melting temperature to be around T ∼ 9 K. We arrive to that conclusion by computing the pressure as a function of temperature for different densities (at and below ρ 0 ), and by verifying that the system retains solid order even at negative pressure, all the way up to T = 8 K. No evidence is seen of any metastable, under-pressurized fluid phase; on lowering the density, the system eventually breaks down in solid clusters, just like in H 2 [5]. Above 8 K, the pressure of the solid phase jumps and stable fluid phases appear at lower densities, signaling the occurrence of melting. These fluid phases do not display any superfluid properties; exchanges of two or three particles occur with a frequency of approximately 0.1%. It is known that quantum mechanical effects in H 2 contribute about one half of the Lindemann ratio at melting [26]; we did not perform the same calculation for this system, but in this system quantum mechanics should be more important, on account of the lighter mass. However, qualitatively melting appears to occur very similarly as in H 2 . Also, our results suggest that thermal expansion, which is negligible [27] in solid H 2 , is likely very small in this crystal.

B. Clusters
It is also interesting to study the physics of nanoscale size clusters of HMu, and compare their behavior to that of parahydrogen clusters, for which, as mentioned above, superfluid behavior is predicted at temperature of the order of 1 K, if their size is approximately 30 molecules or less. Crystalline behavior emerges rather rapidly for parahydrogen clusters of greater size [28], with "supersolid" behavior occurring for specific clusters [29]. This calculation is carried out with the same methodology adopted for the bulk phase of the system, the only difference being that the simulation cell is now taken large enough that a single cluster forms. In this respect, the behavior of HMu clusters is closer to that of parahydrogen than 4 He clusters, in that no external potential is required to keep them together (as in the case of 4 He), at least at sufficiently low temperature ( 4 K).
However, the reduced molecular mass makes the physics of HMu clusters both quantitatively and qualitatively different from that of parahydrogen clusters. The first observation is that the superfluid response is greatly enhanced; specifically, it is found that clusters with as many as N = 200 molecules are close to 100% superfluid at T = 1 K, and remain at least 50% superfluid up to T 4 K, at which point they begin to evaporate, i.e., they do not appear to stay together as normal clusters. This is in stark contrast with parahydrogen clusters, where exchanges are rare in clusters of more than 30 molecules at this temperature, and they would involve at most ∼ 10 particles. Consistently, the superfluid response is insignificant in parahydrogen clusters, and they stay together almost exclusively due to the potential energy of interaction. Our results highlight the importance of the energetic contribution of quantum-mechanical exchanges in the stabilization of HMu clusters at low temperature. Indeed, at T = 1 K we observe cycles of exchanges involving all of the particles in HMu clusters comprising as many as 200 molecules. Also shown for comparison is the same quantity for a cluster of N = 30 parahydrogen molecules, at the same temperature; this cluster has no significant superfluid response. Fig. 7 shows density profiles computed with respect to the center of mass of the system for three different clusters of HMu, comprising N = 30, 100 and 200 molecules. As mentioned above, all these clusters are essentially fully superfluid at this temperature. These density profiles are qualitatively similar, nearly featureless and reminiscent of those computed for 4 He droplets [30]. The comparison with the density profile computed for a cluster of 30 parahydrogen molecules, also shown in Fig. 7, illustrates how the latter is much more compact and displays pronounced oscillations, which are indicative of a welldefined, solid-like shell structure.
As the number of molecules increases, clusters of HMu ought to evolve into solid-like objects, their structure approaching that of the bulk. The calculation of the radial superfluid density [25,31] suggests that crystallization begins to occur at the center of the cluster; for example, at T = 1 K the largest HMu cluster studied here (N = 200) displays a suppressed superfluid response inside a central core of approximately 5Å radius, in which crystalline order slowly begins to emerge, while the rest of the cluster is essentially entirely superfluid. We expect the non-superfluid core to grow with the size of the cluster. In other words, large clusters consist of a rigid, insulating core and a superfluid surface layer, with a rather clear demarcation between the two. This is qualitatively different from the behavior observed in parahydrogen clusters, in which crystallization occurs for much smaller sizes, and "supersolid" clusters simultaneously displaying solid-like and superfluid properties can be identified [29].

IV. CONCLUSIONS
We have investigated the low temperature phase diagram of a bulk assembly of muonium hydride molecules, by means of first principle quantum simulations. Our model assumes a pairwise, central interaction among HMu molecules which is identical to that of H 2 molecules. By a mapping of the H 2 inter-atomic potentials derived from ab initio calculations, we showed that this model provides a reasonable description of the HMu-HMu interation. It is certainly possible to carry out a more accurate determination of the pair potential, but the main effect of the lower µ + mass should be that of increasing the diameter of the repulsive core of the interaction, in turn suppressing exchanges even more. As illustrated in Ref. 15, exchanges of identical particles play a crucial role in destabilizing the classical picture in a many-body system; when exchanges are suppressed, quantum zero point motion can only alter such a picture quantitatively, not qualitatively. Alternatively, this can be understood by the fact that Λ is decreased as the hard core radius is increased, making the system more classical. Obviously, regarding the interaction as spherically symmetric is also an approximation, but one that affords quantitatively accurate results for parahydrogen [20], for which the potential energy of interaction among molecules plays a quantitatively more important role in shaping the phase diagram of the condensed system.
Perhaps the most significant observation of this study is the different physics of bulk and nanoscale size clusters of HMu. Bulk HMu is very similar to H 2 ; despite its very low density (lower than that of superfluid 4 He), the equilibrium crystalline phase is stable below a temperature of about 9 K, much closer to the melting temperature of H 2 (13.8 K) than the mass difference may have led one to expect. No evidence is observed of any superfluid phase, neither liquid nor crystalline, underscoring once again how, in order for a supersolid phase to be possible in continuous space, it requires some physical mechanism to cause a "softening" of the repulsive part of the pair potential at short distances [32], even if only along one direction, like in the case of a dipolar interaction [33].
On the other hand, clusters of HMu including up to a few hundred molecules display superfluid behavior similar to that of 4 He clusters. This suggests that, as the value of the quantumness parameter Λ approaches Λ c from below, one may observe nanoscale superfluidity in clusters of rather large size.
We conclude by discussing the possible experimental realization of the system described in this work. The replacement of elementary constituents of matter, typically electrons, with other subatomic particles of the same charge, such as muons [34], has been discussed for a long time, and some experimental success has been reported. Even a bold scenario consisting of replacing all electrons [35] in atoms with muons (the so-called "muonic matter") has been considered; recently a long-lived "pionic helium" has been created [36]. Thus, it also seems plausible to replace a proton in a H 2 molecule with an antimuon; indeed, muonium chemistry has been an active area of research for several decades [37]. In order for "muonium condensed matter" to be feasible, a main challenge to overcome is the very short lifetime of the µ + , of the order of a µs.