Theoretical study of H2 in a two-dimensional crystalline matrix

We have carried out extensive path integral Monte Carlo simulations of two-dimensional para-hydrogen (p-H2) embedded in a crystalline matrix of alkali atoms. Our results show that, at low temperatures (⪅5 K), the thermodynamically stable phase of p-H2 is a solid, commensurate with the underlying lattice. A nonzero superfluid signal for p-H2 is observed in simulated systems of very small size (e.g. 13 molecules); however, results for systems of larger sizes are altogether consistent with absence of superfluidity in the thermodynamic limit.

temperature (T ∼14 K) significantly higher than that (∼6 K) at which such phenomena as Bose condensation and superfluidity are expected to occur. This is due to the depth of the attractive well of the potential between two hydrogen molecules, significantly greater than that between two helium atoms.
The question of whether p-H 2 would display a SF transition at low temperature, should a sufficiently clever procedure be imagined to suppress crystallization, is clearly of great fundamental interest. Currently, helium is the only known substance displaying superfluidity in the condensed (i.e. non-gaseous) phase. Establishing the presence of superfluidity in at least another condensed matter system might provide further understanding of this phenomenon, possibly suggesting new ways of searching for it in other systems.
The currently accepted theoretical picture of superfluidity relates it to long permutation cycles among identical particles [2,3]; these cycles, in turn, arise from elementary exchanges involving pairs, triplets, etc of particles. In order for these elementary exchanges to occur, the particle thermal wavelength λ T ∝ T −1/2 must be comparable to the average interparticle distance. Thus, superfluidity can only occur at sufficiently low temperature. However, if at low T the system crystallizes, quantum exchanges are suppressed, as particles enjoy a high degree of localization in the crystal phase. Indeed, for this reason, superfluidity and crystalline order have long been thought to be incompatible [4,5]. Although this notion has been challenged by some authors [6]- [8], who have maintained that superflow could occur in a solid as well, and the claim has been recently made of the actual observation of a supersolid phase of Helium [9], the conventional wisdom is that the liquid phase of matter is more likely to display superfluidity.
A fair amount of experimental and theoretical effort is thus aimed at identifying ways of stabilizing a liquid phase of p-H 2 at low temperature. Several attempts have been made [10]- [13] to supercool bulk liquid p-H 2 , but the search for a SF phase has so far not met with success. Confinement and reduction of dimensionality are widely regarded as plausible avenues to the stabilization of a liquid phase of p-H 2 at temperatures sufficiently low, such that a SF transition may be observed. Indeed, computer simulations yielded evidence of SF behaviour in very small (less than 20 molecules) p-H 2 clusters [14], and its experimental observation has been reported [15].
A different, intriguing suggestion was made by Gordillo and Ceperley [16] (GC), namely that superfluidity ought to occur in a strictly two-dimensional (2D) p-H 2 fluid embedded in a crystalline matrix of alkali atoms. Their contention, based on path integral Monte Carlo (PIMC) simulations, is that the presence of the underlying lattice of foreign atoms, incommensurate with the equilibrium crystal structure of pure p-H 2 , should cause a substantial reduction of the equilibrium density of the 2D fluid of p-H 2 molecules, stabilizing a liquid phase, turning SF at T ∼ 1 K.
While undoubtedly intriguing, the above predictions are based on simulations of very small systems (of the order of 10-15 p-H 2 molecules). While in principle it is possible to perform finitesize scaling of the results, attempting extrapolation to the thermodynamic limit, this procedure is not always reliable, particularly with such small system sizes and given the typical statistical uncertainties that occur in Monte Carlo simulations. Thus, it seemed worthwhile to us, to try and assess the robustness of GC's contention by performing PIMC simulations of the same system, but with up to ten times more particles. Seven years after GC's original work, such simulations can be performed rather comfortably on commonly available workstations. We modelled our system of interest exactly as in [16], namely, as a collection of N pointlike quantum particles (p-H 2 molecules) confined to moving in 2D, and in the presence of an underlying lattice of impurities, represented by fixed particles, interacting with the p-H 2 molecules via a chosen potential. Although it has been suggested [17] that the potential selected in [16] may not accurately describe any known interaction, we regard it as an interesting problem to identify conditions (if any) under which p-H 2 superfluidity is possible, even though they may be regarded as somewhat 'artificial', and, at least presently, difficult to realize experimentally.
Our results suggest a different interpretation than that offered by GC, regarding the structure in the equilibrium phase. We contend that the equilibrium phase of the system is not a liquid, as argued in [16], but rather crystalline, commensurate with the underlying lattice of impurities. Thus, if this system should indeed display superfluidity, it would be a supersolid, based on the commonly accepted definition of a supersolid as a phase of matter simultaneously featuring long-range crystalline order and dissipationless flow. However, in our study we only obtained a finite SF signal when simulating the same, small-size system as in GC's original calculation. On systems of larger size, not only does the SF signal disappear, but we also observe a substantial suppression of permutations among p-H 2 molecules. All of this suggests that the commensurate crystal is not a superfluid in the thermodynamic limit.
The remainder of this paper is organized as follows: in the next two sections we introduce the model and provide computational details; in the following section, we illustrate our results and provide theoretical justifications for our analysis. Finally, we outline our conclusions, and motivate our overall pessimism, regarding the possibility of observing p-H 2 superfluidity in systems such as the one considered here.

Model
As mentioned in section 1, we model our system of interest as a collection of N point-like particles (p-H 2 molecules) of mass m, moving in two dimensions in the presence of an external potential arising from a lattice of static, identical impurities. The system is enclosed in a rectangular simulation cell of sides L x , L y (and area A = L x × L y ), with periodic boundary conditions in all directions.
The quantum-mechanical Hamiltonian of the system is the following: Here, V is the interaction potential between any two p-H 2 molecules, only depending on their relative distance r ij ≡ |r i − r j |; the accepted Silvera-Goldman [18] potential is used to describe this interaction. As there appears to be some ambiguity in the literature, regarding the values of the parameters of this potential, for the sake of clarity we list in table 1 the values used in this calculation (see also [19]). The potential takes the form: where x = r/r m and F(x) = exp{−[(D/x) − 1] 2 }, if x < D, and 1 otherwise. Parameters of the Silvera-Goldman potential describing the interaction between two p-H 2 molecules (equation (2)). All parameters are dimensionless, except for (given in K) and r m (given in Å).
0.07468938 C 10 0.38990868 31.5763295 The system also includes M impurities, positioned at regular lattice sites R σ , with σ = 1, 2, . . . , M. In this work, we only considered a triangular lattice arrangement of impurity sites; 1 the concentration of impurities is given by M/A, and is set at 0.01155 Å −2 . This results in a lattice constant of 10 Å for the triangular lattice formed by the impurities. The nominal density The interaction between a p-H 2 molecule and an impurity (i.e. the U term in (1)) is described by a Lennard-Jones potential with parameters = 9.54 K and σ = 3.75 Å, just like in [16]. This potential, with these parameters, was chosen in [16] with the aim of providing a reasonably realistic description of the interaction of a p-H 2 molecule with a potassium atom.
Different values of the parameters and σ have been proposed elsewhere for Lennard-Jones potentials, describing the interaction of a p-H 2 molecule with alkali atoms [20]. In particular, a more realistic value of for the interaction between a p-H 2 molecule and a potassium atom may be as much as three times greater than that used here, whereas σ may be underestimated in this work by about 20%. In principle, there is also the question of whether a 2D lattice of potassium atoms such as the one considered here may arise in nature.
In the context of this work, these issues (namely, the potential and the impurity lattice) do not seem crucial, as we are primarily interested in exploring the theoretical problem described by (1), rather than in aiming at providing results to be compared directly to experiment. We note, finally, that, as suggested in [16], it may be feasible to produce a lattice of alkali atoms such as the one described here, by adsorbing a fraction of a monolayer of alkali metal atoms (Rb, Cs and K) on a Ag(111) or on a graphite substrate [21,22].
There is also, of course, the question of whether a fully 3D treatment of the problem may result in different physics. It is conceivable, for example, that zero-point motion in the direction perpendicular to the substrate might strengthen the disordered (liquid) phase. The degree of confinement of hydrogen molecules in the direction perpendicular to the substrate depends, of course, on the strength of the underlying substrate. A recent study of p-H 2 films adsorbed on lithium (a relatively weakly attractive substrate) has yielded evidence that such a zero-point motion, while significant, is not nearly as large as for helium films, and does not dramatically alter the properties of the system with respect to a fully 2D treatment [25].

Calculations
We studied physical properties of the system described by (1), in the range of temperature 1-4 K, using PIMC simulations based on the standard procedure [24]; technical details of this study are the same as in [25]. In particular, we have utilized a high-temperature density matrix [26], accurate up to fourth-order terms in the expansion of the exact density matrix in powers of the imaginary time step τ. Also, the algorithm utilized here includes the all-important sampling of permutations of p-H 2 molecules, necessary to study superfluidity. We extensively tested our code by reproducing the results of Pollock and Ceperley for liquid 4 He in two and three dimensions [27,28]. We utilized three simulation cells in our study, all with periodic boundary conditions; the smallest one is identical with that used by GC, namely 20 × 17.32 Å 2 , which accommodates four impurities; we also used a simulation cell measuring 40 × 34.64 Å 2 , which accommodates 16 impurities, and one measuring 50 × 51.96 Å 2 , with 30 impurities.
We observed convergence of the physical estimates with a value of τ = τ 0 = 1/320 K −1 ; figure 1 illustrates this point with results for the energy per p-H 2 molecule obtained at T = 2 K, with different values of the time step τ (in K −1 ). The density is θ = 0.0381 Å −2 , and the system utilized for this test calculation is the one with the smallest size (i.e. four impurities). Figure 2 shows computed energetics of 2D p-H 2 films at different densities. Our results are for a temperature T = 2 K; we find that the results are little changed below T ∼ 4 K. In all of our calculations, we estimated the contribution to the potential energy attributed to particles outside the main simulation cell by assuming that the pair correlation function g(r) equals one outside the cell. Inspection of the results for g(r) shows that this is not a very good approximation, but we estimate the resulting error to be at the most of the order of our statistical uncertainties (for the largest system size studied). Comparison of our results with those of GC shows broad agreement, as far as the general trend and the location of the equilibrium density (θ e = 0.0381 Å −2 ) are concerned. There are, however, significant numerical discrepancies between our results and GC's, ours being consistently lower by a few K. This is puzzling, particularly considering that our estimates are obtained at a higher temperature; in fact, even our estimates at 4 K are still considerably lower than those of GC's at 1 K. Although our results pertain to systems of size 7-8 times greater than those of GC's calculation, a numerical discrepancy of this magnitude cannot be attributed to finite-size corrections. In order to illustrate this point, in figure 2 we show our energy estimate (star) at the equilibrium density, for a system of 13 p-H 2 molecules, i.e. identical in size with that of GC. Such an estimate is in close agreement with the one that we obtain for a larger (99 particles) system.

Results
Without knowing all the details of the calculations performed in [16], we may only speculate as to what could cause such a significant difference in the estimates of a basic quantity as the energy. One possible source of error in PIMC simulations is the choice of the time step. While in principle one should always extrapolate all estimates to the τ → 0 limit, this procedure is quite time-consuming, and seldom carried out in practice. Normally, one performs all simulations using a fixed value of time step τ 0 , chosen sufficiently small that estimates obtained with this time step may be expected to coincide with extrapolated ones, within the statistical uncertainties. In our study, we ascertained (figure 1) that the time step τ 0 utilized would not result in a greater systematic error on the energy estimates than ∼0.1 K, which is much less than the difference between GC's results and ours. It is possible than GC's results may be affected by a more significant time step error. An alternative, perhaps simpler, explanation would attribute the discrepancy to slightly different values for the parameters of the p-H 2 Silvera-Goldman potential, utilized in the two calculations. We finally mention that, using our PIMC technology, and the Silvera-Goldman potential with parameters given in table 1, we have obtained results in the T → 0 limit, both for pure 2D p-H 2 [23] as well as in the presence of impurities, in quantitative agreement with ground state quantum Monte Carlo calculations for the same systems, performed by other authors [17]. Thus, we are quite confident in the accuracy of the results presented here. Figure 3 shows an instantaneous configuration of 99 p-H 2 molecules at the equilibrium density θ e . The arrangement of p-H 2 molecules on a regular (essentially kagomé) crystal is clear, and can be confirmed by a calculation of the static structure factor for a lattice vector associated to such crystalline arrangement or by the pair correlation function g(r) (shown in figure 4). These results are not qualitatively or quantitatively different from those reported in [16]; however, the conclusion of the authors (GC) therein, is that p-H 2 is a liquid at the equilibrium density θ e .
Their (GC's) argument is based on a comparison of the static structure factor calculated for the system in the presence of impurities and for pure p-H 2 at the same density (figure 2 of [16]). The observed similarity of the structure factor for these two cases (except for those wave vectors associated to the underlying impurity lattice), and the presumption that pure p-H 2 should be a (metastable) liquid at this low density (which is well below the equilibrium density of 2D p-H 2 , namely 0.067 Å −2 , at which the system is crystalline; see, for instance, [23]), suggest that p-H 2 is a liquid, even in the presence of the impurity lattice.
We propose here a different assessment than that of GC. It has been recently shown that no low-temperature metastable liquid phase of pure 2D p-H 2 exists. The spinodal density of 2D p-H 2 in the T → 0 limit is approximately 0.059 Å −2 ; the system is a crystal at equilibrium, remains a solid down to the spinodal density and breaks down into solid clusters at lower densities [23]. Moreover, little or no overlap is seen among paths, suggesting that permutations are suppressed, which is typical of quantum crystals, and has been consistently found to be the case in PIMC studies of low-dimensional p-H 2 at low temperature, even in low density metastable phases Pair correlation function of p-H 2 molecules, computed by PIMC at T= 1 K and at the equilibrium density (θ e = 0.0381 Å −2 ) on a 99-particle system. Distances (r) are given in Å. [23,25,29]. Indeed, p-H 2 molecules remain spatially localized, as shown in figure 3. All of this, in our view, suggests that p-H 2 molecules, which are prevented by the impurities from forming the preferred crystal arrangement at the equilibrium density, will do the 'next best thing', namely form a crystal commensurate with the underlying impurity lattice. Such a crystal is shown in figure 3. Therefore, if p-H 2 were indeed a superfluid in the physical arrangement considered here, it ought to be regarded as a 'supersolid', based on the commonly accepted definition of the supersolid phase as one featuring, simultaneously, crystalline long-range order and dissipationless flow. The question is, of course, whether the system of our interest is, in fact, SF.
PIMC simulations allow one to search for superfluidity in a given many-body physical system, at least in principle, by performing a direct computation of the SF fraction ρ s [24]. Using the 'winding number' estimator [25], we have computed ρ s at the equilibrium density for systems with 13, 53 and 99 p-H 2 molecules. The smallest system is enclosed inside the rectangular cell shown in figure 3 (lower left corner). For the 13-molecule system, we find a finite value of ρ s at T = 1 K, namely ρ s = 0.22 ± 0.02, consistent with that found in [16]. However, we obtain ρ s = 0 as soon as the system size is increased to 53 molecules. Specifically, we observe the complete disappearance of permutations among p-H 2 molecules on systems of size greater than 15 particles.
In and of itself, the fact that ρ s → 0 in PIMC simulations of finite systems of increasing size, is not sufficient to conclude that ρ s = 0 in the thermodynamic limit. The SF signal in PIMC is related to long cycles of permutation of p-H 2 molecules, which are sampled less and less efficiently as the system size is increased. Thus, one may simply be facing an ergodicity problem. However, a careful analysis of the results leads us to surmise that, in this case, the SF density in the thermodynamic limit is indeed zero at the temperatures considered here (down to 1 K), and quite likely all the way to T = 0.
In our PIMC simulation, p-H 2 permutation cycles of arbitrary length can occur as (i) the result of direct sampling, as well as (ii) the global result of a number of smaller permutation cycles, e.g. involving two or three particles only (see, for instance, [24]). On histogramming both the probability of acceptance P acc (n) of permutation cycles of length n, as well as the probability of occurrence P(n) of those cycles in the simulation, we find that, in a PIMC simulation with 13 p-H 2 molecules (the only one yielding a finite value of ρ s ), both probability distributions feature a sharp peak at n = 5 (see figure 5). Permutations and cycles of two and three particles, instead, occur exceedingly rarely. The SF signal obtained for a system of 13 p-H 2 molecules is clearly related to the cycles of length 5, i.e. the only ones that can result in a nonzero winding, as they span the entire simulation box. These cycles exclusively occur through direct sampling, rather than as a result of smaller permutation cycles, involving, for example, two or three molecules. They can be obtained, both horizontally as well as vertically, simply by having each particle 'hop' into the position of one that is immediately adjacent, which in turn does the same, in the same direction. The last particle in the cycle hops into the position of the first one, across the periodic boundary. An example is shown in figure 6. In our simulations with greater number of molecules, the equivalent cycles never occur; with 53 particles, the corresponding cycles would consist of 10 particles, whereas with 99 particles a cycle of the above type would involve 12-15 particles. Although our algorithm can construct such long cycles, acceptance is found to be zero in all of our simulations.
There are, in our view, two reasons to believe that this mechanism will not contribute a finite SF signal at any finite temperature, in the thermodynamic limit. Firstly, because these permutation cycles only occur as a result of direct sampling, rather than through many, more frequently sampled smaller cycles (as is the case, for example, in liquid 4 He), their contribution to the physical observables should vanish in the thermodynamic limit. Secondly, because they are essentially 1D in character, i.e. the same physical mechanism yields a finite SF density in a finite 1D system, for which the critical temperature is known to be zero.

Conclusions
Based on our extensive PIMC study of 2D p-H 2 embedded in a crystalline matrix of alkali atoms, modelled in exactly the same way as in a previous study, we conclude that this system is not a candidate for observing superfluidity in p-H 2 . At low temperature, the system forms instead a 2D crystal, commensurate with the underlying lattice of impurities. We have presented results for triangular impurity lattices, but the same results were seen with other lattices as well (e.g. rectangular). The finite SF signal obtained by other authors can be attributed, in our view, to the finite size of the simulation cell used in [16]. It is unlikely that such a null outcome could be significantly changed by a 'fine tuning' of the potential parameters describing the interaction of p-H 2 molecules with the impurities. For example, similar calculations carried out on the same physical system considered here, but with very different potentials, have essentially confirmed our findings [17].
What are, therefore, the prospects of observing superfluidity of p-H 2 in 2D films? The absence of a metastable liquid phase of p-H 2 appears to render the investigation more complex than it may have been originally anticipated. It is not obvious how incommensuration, or even disorder, may reduce the tendency of the system to solidify. It seems that one will always end up with a commensurate crystal or, at best, a 'glassy' solid. And, given the scarce propensity of the pure system to display permutations among molecules, it is unclear how disorder or incommensuration may enhance particle exchange.
Nonetheless, it remains possible that a different substrate may indeed stabilize a liquid phase of molecular hydrogen. In fact, recent experiments [30] on D 2 films, co-adsorbed on graphite and preplated by a monolayer of Kr, have provided some evidence that a liquid phase of D 2 may be stable, down to a temperature as low as 1.5 K. At such a low temperature, one would certainly expect p-H 2 to turn SF, and because p-H 2 molecules are lighter than D 2 , it is conceivable that a liquid phase of p-H 2 should be observable as well, in the same physical conditions. This system seems therefore quite worthy of further experimental and theoretical study.