Quantum Monte Carlo study of thin parahydrogen films on graphite

The low-temperature properties of one and two layers of parahydrogen adsorbed on graphite are investigated theoretically through Quantum Monte Carlo simulations. We adopt a microscopic model that explicitly includes the corrugation of the substrate. We study the phase diagram of a monolayer up to second layer promotion, and the possible occurrence of superfluidity in the second layer. We obtain results down to a temperature as low as 8 mK. We find second-layer promotion to occur at a considerably greater coverage than obtained in previous calculations and estimated experimentally; moreover, we find no evidence of a possible finite superfluid response in the second layer, disproving recent theoretical predictions.


Introduction
The physics of a thin (few layers) film of parahydrogen molecules (p-H 2 ) adsorbed on a strongly attractive corrugated substrate like graphite, has been extensively studied experimentally [1,2,3,4,5] and theoretically [6,7,8,9,10,11].The basic features of the phase diagram of the system at low temperature are well understood, but there exist specific aspects which still require elucidation.
Adsorption of p-H 2 on essentially any substrate [12] occurs through the formation of successive crystalline layers.On a strongly attractive substrate such as graphite, the first layer is physically distinct, featuring an equilibrium phase registered with the underlying carbon lattice, which transitions at higher coverage to a compressible, incommensurate crystal through a series of domain-wall phases.It is interesting to compare p-H 2 with 4 He, i.e., the archetypal (Bose) quantum fluid, and it seems fair to state that there is little qualitative difference between the phase diagram of the first layer of either substance on graphite, as the energetics is in both cases dominated by the substrate attraction.On the other hand, the second layer displays very different behavior for the two systems, with 4 He forming a liquid (superfluid) film while p-H 2 remains solid.
An outstanding puzzle is represented by the value of coverage θ p of a p-H 2 monolayer at which molecular promotion to the second layer begins to occur.There exists an experimental estimate [3] for θ p ≈ 0.094 Å −2 , considerably lower (by about 15%) than the corresponding coverage for a 4 He monolayer on the same substrate [13].This seems counter-intuitive, as a) one would expect a parahydrogen monolayer to be more compressible than one of helium, which displays much more pronounced quantum effects [14], and b) the attractive well of the effective potential between a p-H 2 molecule and a graphite substrate [15] is about three times deeper than that for a 4 He atom and the same substrate [16].
Microscopic calculations for realistic models of the system have yielded conflicting results; specifically, an early study [10] making use of Quantum Monte Carlo (QMC) simulations predicted second-layer promotion at a coverage very close to that reported in Ref. [3], but a subsequent study [17] of the same system based on a slightly different QMC methodology showed no promotion to second layer up to a coverage ∼ 0.110 Å −2 , i.e., very similar to that observed for a 4 He monolayer.Both calculations treated the underlying graphite substrate as smooth, i.e., ignored its corrugation.While the cause of the disagreement between the two calculations is unclear, it is conceivable that corrugation could have a significant effect on θ p .However, in Ref. [17] the same calculation was carried out for a p-H 2 monolayer on graphene (which is some 10% less attractive than graphite), explicitly including substrate corrugation, and the second layer promotion coverage was found to be close to 0.11 Å −2 .One certainly expects a comparable or even greater value for the more attractive graphite.
A second point of contention has to do with the second (solid) layer, for which the claim of a possible finite superfluid response in the ground state has been recently made [11].This seems very similar to other predictions of superfluidity (SF) of either solid helium or hydrogen in various settings (usually in confinement) arrived at using the same ground state methodology, all of which have been shown to be incorrect.There is presently no experimental evidence of superfluid behavior of p-H 2 , except (possibly) in very small clusters (∼ 15 molecules or less) [18]; all reliable theoretical studies of parahydrogen, including in reduced dimensions [19], thin films [20] and in disordered and/or inhomogeneous environments [21,22,23] have shown that exchanges of indistinguishable molecules, which are known to underlie SF, are strongly suppressed even at low temperature.However, no microscopic calculation for the second layer of p-H 2 on graphite has yet been carried out, independent of that of Ref. [11], in which the corrugation of the graphite substrate is taken into account.Although a recent study [24] for 4 He provides strong evidence that the effect of substrate corrugation on the superfluid properties of the second layer is minimal, it seems worthwhile to investigate in greater depth the scenario of a possible superfluid response on the second layer of p-H 2 on graphite, taking the corrugation of the substrate into account.
To clarify the above outstanding issues, we have carried out extensive QMC simulations of a thin (1 and 2 layers) film of p-H 2 adsorbed on graphite at low temperature (as low as 8 mK), making use of the same ab initio potential proposed in Ref. [10], describing the interaction of a p-H 2 molecule with the substrate, specifically designed to capture the effects of substrate corrugation.We aim to compare directly our results to those obtained therein.For the 2-layer system, we focused on ascertaining whether there is a significant enhancement of quantum exchanges of indistinguishable p-H 2 molecules at subKelvin temperatures, and with that the possible onset of a weak superfluid response as predicted in Ref. [11], at the same thermodynamic conditions.
Our results yield, for the monolayer system, a second-layer promotion coverage close to 0.110 Å −2 , i.e., consistent with both that found on graphene as well as on a smooth graphite substrate [17], significantly above the existing experimental estimate from Ref. [3].While our results are generally qualitatively consistent with those of Ref. [10], there are significant numerical discrepancies, which we believe may be at the root of the difference between the promotion coverage pre-dicted therein and in this work.Still, the disaccord of our calculation with experiment remains, and at this point can only be solved by new, independent experimental and theoretical studies.
No evidence of any "supersolid" phase is observed for the second layer.Our simulations, carried out down to a temperature as low as 8 mK, yield the same paucity of quantum-mechanical exchanges of indistinguishable molecules observed at temperatures of the order of 1 K, down to temperatures more than two orders of magnitude lower; no discernible difference can be seen in the one-body density matrix computed at temperature T = 4 K and T =0.03 K, once again confirming that this system forms an insulating crystal which does not undergo a superfluid transition in the T → 0, which is the conclusion reached in many previous theoretical studies, consistent with all existing experimental observations.It is therefore our submission that the prediction of Ref. [11] is incorrect.
The remainder of this paper is organized as follows: in section 2 we describe the microscopic model of the system; in Sec. 3 we briefly describe our methodology; we present and discuss our results in Sec. 4 and finally outline our conclusions in Sec. 5.

Model
We consider an ensemble of N p-H 2 molecules, regarded as point-like spin-zero bosons, moving in the presence of a graphite substrate, modeled as described below.The system is enclosed in a simulation cell of sizes L x × L y × L z , with periodic boundary conditions in all directions (but L z is taken large enough to make boundary conditions in the z direction irrelevant).The graphite substrate lies on the z = 0 plane; in our monolayer simulations, we set L x = 34.08Å, L y = 36.8927Å; on the other hand, those of a two-layer system were carried out for a significantly smaller simulation cell, one with L x = 21.3Å and L y = 19.6761Å.The reason is that our goal in this case is exclusively that of ascertaining whether the conditions for SF exist, and the superfluid response is typically enhanced in a finite system with periodic boundary conditions.Thus, if no superfluid signal is observed in a system of small size, a fortiori that will be the case in the thermodynamic limit.
The quantum-mechanical many-body Hamiltonian reads as follows: The first and third sums run over all the N p-H 2 molecules, λ = 12.031 KÅ 2 ; the second sum runs over all pairs of molecules, r i j ≡ |r i − r j |, r i ≡ (x i , y i , z i ) being the position of the ith molecule, and v(r) is the accepted Silvera-Goldman pair potential [25], which describes the interaction between two p-H 2 molecules.V is the potential describing the interaction of a p-H 2 molecule with the graphite substrate; for consistency with the calculation of Ref. [10], we use the same potential utilized therein, which can be expressed as follows: V 0 (z) is a term that only depends on the distance of the molecule from the basal plane; it corresponds to the laterally averaged potential normally utilized to represent a flat substrate [15].The sum in the second term of the right-hand side of ( 2) runs over all reciprocal lattice vectors G ≡ (G x , G y , 0) of the graphite substrate.The functions V 0 (z), V G (z) are provided in Ref. [10].We come back to the details of the evaluation of (2) in our simulations in Sec. 3.

Methodology
The QMC methodology adopted here is the canonical [26,27] version of the continuous-space Worm Algorithm [28,29], a well-established finite temperature QMC technique.Details of the simulations carried out in this work are standard, and therefore the reader is referred to the original references.However, because as mentioned above there are numerical discrepancies between our results and those reported in Ref. [10], where a similar computational technique was utilized, we are going to provide here enough technical information to enable others to repeat the calculation, in order to investigate the possible causes of the disagreement.
The calculation of the potential energy of interaction between a p-H 2 molecule and the graphite substrate (Eq.2) has been carried out in the simulation by tabulating and interpolating the functions V 0 (z) and V G (z), and by performing a direct, on-the-fly summation over terms associated to a subset of reciprocal lattice vectors; we found that the inclusion of the twelve shortest reciprocal lattice vectors is sufficient to achieve a numerically accurate representation of V(r).
The key physical quantities computed in this work are the energetics and the superfluid fraction ρ S (T ) of the top layer as a function of temperature, for which we use the well-known winding number estimator [30].We also evaluate structural properties, such as the averaged molecular density profile along the direction z, perpendicular to the substrate, in order to investigate the occurrence of promotion to the second layer.Crystalline order in both layers may also be monitored through the visual inspection of the imaginary time paths.
Simulations of a two-layer system have been carried out by a) considering molecules in the top and bottom layer as two distinct species and b) regarding molecules in the bottom layer as distinguishable quantum particles (i.e., "Boltzmannons").Both approximations are justified by the de facto absence of quantum-mechanical exchanges involving molecules in the bottom layer; on the other hand, molecules in the top layers are considered indistinguishable and their quantum (Bose) statistics is fully accounted for (as we shall see, even changes of molecules in the top layers turn out to be exceedingly infrequent even at the lowest temperature considered here).
We used the primitive approximation for the short imaginary-time (τ) propagator and carried out extrapolation of all the computed physical observable in the τ → 0 limit.This point warrants a quantitative discussion because this is where the first major difference with Ref. [10] is noted.Fig. 1 show the computed energy per p-H 2 molecule at temperature T = 1 K as a function of the time step τ, for coverage θ = 0.0636 Å −2 , namely the coverage at which a commensurate monolayer crystal forms, as will be shown below.Because the primitive approximation is used, one expects quadratic behavior as a function of τ, in the τ → 0 limit.As shown in the figure, the estimate for τ = 1.5625 × 10 −3 K −1 is indistinguishable, within statistical uncertainties, from the ex-trapolated one, namely −480.3(4)K.This is only a few K lower than the value reported in Ref. [10], namely −476.9K, but what is puzzling is that the value of the time step used therein was 0.005 K −1 , which in our calculation yields an energy estimate some ∼ 17 K lower than the value reported in Ref. [10].The reason for this disagreement is unknown to us, at this time, as the calculations are based on the same potentials; it cannot, in our view, be attributed to the different system sizes (the calculation of Ref. [10] used N = 60 molecules, whereas N = 80 was the number in this work, at this coverage), nor to the difference in temperature between the two calculations (we come back to this point below), nor to a possible, different choice of wave vectors included in the expansion of the potential (2).
The energy estimates quoted below for the energy were all obtained either by carrying out explicitly the τ → 0 extrapolation, as shown above, or using a time step τ = 1.5625 × 10 −3 K −1 .On the other hand, it has been observed in this work, as well as in 4 He studies based on the same methodology [24], that convergence of the results for structural and superfluid properties is achieved with a considerably greater (up to four times) time step than that needed for the energy.

Monolayer
We begin by discussing the energetics of a monolayer.In general, we observe that the estimates for the energy, as well as for all structural properties, remain essentially unchanged (the energy within less than 1%) for temperatures below 4 K.In the zero-coverage limit (i.e., one particle), we obtain an energy per molecule of −454.7(1)K, which is very close to the estimate obtained by Crowell and Brown for the laterally averaged potential [15].
In our monolayer calculations, we cut off the pair potential at 17 Å, i.e., we set it to zero beyond that distance.We estimate the energy contribution arising from pairs of molecules at distances greater than that to amount to less than 0.01% of the total energy.Fig. 2 shows the computed energy per p-H 2 molecule as a function of coverage, at the two temperatures T = 1 K (filled circles) and T = 4 K (filled boxes).Also shown (open circles) are the estimates from Ref. [10].There is an obvious numerical discrepancy between the results obtained here, which on the scale of the figure are temperature independent, and those offered in Ref. [10].The difference is just a few K at the equilibrium coverage θ = 0.0636 Å −2 , which is the same on both works and corresponds to the formation of a solid monolayer commensurate with the underlying substrate, but increases at higher coverage and becomes as large as ∼ 20 K in the proximity of the coverage estimated in Ref. [10] to correspond to second layer promotion.
The origin of such a large disagreement between these two calculations is unclear, especially given that, as shown above, the choice of time step made by the authors of Ref. [10] should have resulted in an underestimation of the energy.While the overall trend is similar in both calculations, the specific physical predictions made in Ref. [10] are, in our view, likely to be affected by the quantitative inaccuracy of their energy estimates.
Altogether, the physics of the monolayer is very similar to that observed on graphene; we did not attempt to characterize it in detail in this work, focusing instead on the issue of second-layer promotion, attempting to elucidate the outstanding disagreement between some of the theoretical calculations and the existing experimental estimates.
Fig. 3 shows the p-H 2 density profile n(z) in the direction perpendicular to the substrate, computed at T = 1 K for a coverage θ = 0.1098 Å −2 , i.e., well above the experimentally estimated second-layer promotion coverage θ p [3].In Ref. [10] it was argued that simulations based on the many-body Hamiltonian (1) yield a value of the second-layer promotion coverage in agreement with experiment, based on an analysis of density pro-  files such as that shown in Fig. 3, displaying a "bump" at distances ∼ 6.5 Å from the substrate for coverages close to, or above the experimentally estimated θ p .Our simulations, on the other hand, give no evidence of that until a much higher value of coverage is reached, slightly above 0.11 Å −2 , which is consistent both with the value found on graphene [17] as well as for 4 He on graphite [31].
Fig. 4 shows an instantaneous configurational snapshot (particle world lines) from a simulation of a p-H 2 monolayer at temperature T = 1 K for coverage θ = 0.1018 Å −2 .Although the world lines extend far out, up to distances considerably greater than the average width of the layer, nevertheless all particles remain localized within it, as shown by the positions of the path centroids.There is no evidence of the formation of an actual second layer.We therefore conclude that promotion to the second layer does not take place for θ ∼ 0.094 Å −2 as contended in Ref. [10], but at considerably higher coverage, ∼ 0.11 Å −2 , value is very close to that predicted for 4 He as well, on the same substrate [31].

Bilayer
We now discuss the results of our simulation of a two-layer p-H 2 system, aimed at investigating the possible occurrence of a finite superfluid response in the top layer, in the low-temperature limit.The possible SF of a film of p-H 2 has been investigated and ruled out in numerous theoretical calculations, characterized by a wide variety of physical settings [20,21,22,23].Ev- ery time the conclusion was arrived at that exchanges among identical molecules, which underlie any superfluid response, are too strongly suppressed in this system, mainly due to the large value of the diameter of the repulsive core of the potential at short distances.Nonetheless, it was recently contended that the top layer of a two-layer film of coverage 0.165 Å −2 may display a small but finite superfluid response in the ground state [11].If confirmed, this would be a remarkable finding, as the top layer still displays crystalline longrange order, which would render this the first example of a supersolid system (there is no indication that the top solid layer would be commensurate with the underlying one).
In order to provide an independent theoretical assess-ment of this intriguing prediction, we carried out simulations at finite temperature of the same system, based on model (1), which is physically equivalent to that utilized in Ref. [11], as explicitly verified in Ref. [24].As mentioned above, we purposefully utilized a relatively small simulation cell and overall number of particles, enabling us to explore as wide a range of temperature as possible, enabling us to detect any sign of build-up of superfluid order.We set the coverage to be the same as that of Ref. [11], namely 0.165 Å −2 , with a twodimensional density equal to 0.100 Å −2 in the bottom layer [32], and 0.065 in the top one, just as in Ref. [11].
It is worth restating that the absence of SF in such a small system which comprises 42 p-H 2 molecules on the bottom layer and 27 on the top one, necessarily implies no SF in the thermodynamic limit.Fig. 5 shows two-dimensional density maps for the bottom and top layers of the simulated two-layer system, at a temperature T =8 mK.The first obvious remark is that the system forms crystal phases on both layers, as expected; p-H 2 molecules are initially given the positions shown in the figure for the bottom layer, namely a triangular lattice incommensurate with the graphite substrate is assumed, consistently with our present experimental and theoretical understanding of the system.Molecules in the top layer, on the other hand, are initially placed at random initial position at a distance z = 6 Å from the graphite basal plane, i.e., the regular crystalline arrangement shown in the right part of Fig. 5 forms spontaneously, despite the fact that the cell is not designed to accommodate a triangular lattice with that number particles; just as observed in other studies, if the system cannot form the preferred structure, due to geometrical constraints, it will do the "next best thing", namely form whatever triangular crystal can be formed inside the cell.
Another thing that can be noticed is that p-H 2 molecules in both layers are highly localized, that there is no visible overlap between quantum delocalization clouds of adjacent molecules to indicate that exchanges among molecules are essentially non-existent, and this is precisely what is observed in the simulations, i.e., even those carried out at a temperature as low as 8 mK show the nearly complete absence of exchanges of p-H 2 molecules.Consequently, no evidence of any finite superfluid response is observed; the superfluid fraction of the top layer is zero (within the precision of the computing machine) even at a temperature as low as T = 8 mK.This observation is consistent with all other studies of p-H 2 performed by us, though it is worth mentioning that, at least to our knowledge 8 mK is the lowest temperature attained so far in QMC simulations of either 4 He or p-H 2 .Fig. 6 compares the laterally averaged one-body density matrix computed for the top layer, at the two temperatures T = 4 K and T = 0.03 K.The most remarkable aspect of this result, to our knowledge, is the virtual indistinguishability of the two results, at least within the scale of the figure.The absence of any noticeable temperature dependence within a temperature range that spans over two orders of magnitude constitutes, in our view, strong evidence that the system is not superfluid, even in the T = 0 limit.One could argue, of course, that the temperature that we considered is not low enough to observe the onset of SF.We do not believe that to be the case, however, as the hypothetical superfluid transition should conform to the Bereszinskii-Kosterlitz-Thouless paradigm, and in particular, fulfill the well-known "universal jump" condition [33], based on which one can obtain a fairly reliable quantitative estimate [34] of the superfluid transition temperature T c .In this case, we can expect T c ≈ 10 mK, i.e., simulations of a system of size as small as that considered here at T = 8 mK should yield clear evidence of quantum-mechanical exchanges of top-layer molecules and/or a finite in-plane superfluid response, if the predictions of Ref. [11] were correct.

Discussion and Conclusions
We have carried out extensive QMC simulations of monolayer and bilayer p-H 2 films adsorbed on a graphite substrate, making use of a microscopic model in which the corrugation of the substrate is fully taken into account.For a monolayer, we report important numerical differences with respect to prior theoretical work [10], the most important being the estimated second-layer promotion coverage, which we find to be some ∼ 15% higher.The resolution of the discrepancy requires additional independent theoretical studies carried out with the same methodology and microscopic model.For the moment, however, we note that significantly greater binding energy of the p-H 2 molecules obtained in this work, at coverages above the equilibrium one.Our estimate is also in disagreement with the only (at least to our knowledge) existing experimental determination of the second-layer promotion coverage for this system; while we cannot comment on experimental details, it seems counterintuitive that second-layer promotion should occur at a lower coverage for a system such as p-H 2 , which experiences a significantly stronger attraction to the graphite substrate and whose behavior is considerably less quantum than 4 He.Nonetheless, it should be mentioned that our result is also at variance with that of a recent ground state calculation [11], based on a different model of the substrate, which predicts a value of the second layer promotion coverage in agreement with the experiment.
We have also investigated the possible occurrence of a finite superfluid response in the top layer of a bilayer system, at the same condition of coverage in which other authors have maintained that SF could occur in the low-temperature limit.Our simulations, reaching a temperature as low as 8 mK, confirm the same observations consistently and repeatedly made for p-H 2 in all physical settings explored so far, namely that this system fails to display SF due to the strong suppression of quantum-mechanical exchanges.This parallels the conclusion reached for 4 He monolayers on graphene, for which a very similar contention of SF was made, based on the same approach [35].In our submission, our calculation disproves the contention of Ref. [11], the latest of a series of incorrect predictions [19,22,23,36,37] of a finite SF response based on a computational technology (Diffusion Monte Carlo) known to be affected by serious limitations, chiefly the bias due to the use of a finite population of random walkers [38,39,40,41,42].
In general, all bulk supersolid phases that have been established theoretically hinge on some "softening" of the repulsive core of the pairwise interaction at short distances [43,44].

Figure 1 :
Figure 1: Energy per molecule at T = 1 K, for coverage θ = 0.0636 Å −2 , computed as a function of the time step τ.Dotted line is a quadratic fit to the data for τ → 0. It is not possible to include in the quadratic fit data points obtained for values of the time step above 1.5625 × 10 −3 K −1 .

Figure 2 :
Figure 2: Energy per p-H 2 molecule on a corrugated graphite substrate computed in this work as a function of coverage θ, at the two temperatures T = 1 K (filled circles) and T = 4 K (filled boxes).Dashed line is the ground state energy at the equilibrium density, attained for coverage θ = 0.0636 Å −2 .Open circles refer to data from Ref. [10] at T = 2 K. Statistical errors are smaller than symbol sizes. o

Figure 3 :
Figure 3: Density profile n(z) (in Å −3 ) in the direction perpendicular to the graphite substrate for a monolayer p-H 2 film of coverage θ = 0.1098 Å −2 , computed by simulation at temperature T = 1 K. Dashed line shows the corresponding profile for coverage θ = 0.0636 Å −2 .

Figure 4 :
Figure 4: Snapshot of instantaneous many-particle world lines for a p-H 2 monolayer of coverage θ = 0.1016 Å −2 at temperature T = 1 K, projected on the x − z plane.

Figure 5 :
Figure 5: Density map of p-H 2 molecules in the bottom (left) and top (layer), at temperature T = 8 mK.The total coverage is 0.165 Å −2 , the two-dimensional density of the bottom layer is 0.100 Å −2 .

Figure 6 :
Figure 6: Laterally averaged one-particle density matrix for the top layer of p-H 2 in a two-layer system of coverage 0.165 Å −2 .Results shown are for the two temperature T = 4 K (circles) and T = 0.03 K (diamonds).Statistical errors are smaller of symbol sizes.