{'fraction_non_alphanumeric': 0.03637427878585166, 'fraction_numerical': 0.01726733004431809, 'mean_word_length': 4.658623136976579, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular hydrogen isotopes adsorbed on krypton-preplated graphite: Quantum Monte Carlo simulations

Adsorption of ortho-deuterium and para-hydrogen films on a graphite substrate, pre-plated with a single atomic layer of krypton, is studied theoretically by means of quantum Monte Carlo simulations at low temperature. Our model explicitly includes substrate corrugation. Energetic and structural properties of these adsorbed films are computed for a range of hydrogen coverages. Thermodynamically stable adsorbed films are solid, with no clear evidence of any liquid-like phase. Quantum exchanges of ortho-deuterium and para-hydrogen are essentially absent in this system, down to zero temperature; consequently, this system displays no superfluidity in this limit. Our simulations provide evidence of a stable domain wall fluid at low temperature, consistently with recent experimental observations.


I. INTRODUCTION
Studies of low temperature adsorption of highly quantal fluids, such as helium or molecular hydrogen are motivated by the fascinating properties that such adsorbed quantum films display, often considerably different than those of the bulk. Fluids of ortho-deuterium (o-D 2 ) and para-hydrogen (p-H 2 ) molecules are interesting physical systems, because these molecules have mass equal to (in the case of o-D 2 ) or half (in the case of p-H 2 ) the mass of a helium atom, rendering zero-point motion significant at low temperature. In addition, each molecule is a composite boson, and therefore it is conceivable that, at low enough temperature, a hydrogenic fluid might display physical behavior similar to that of liquid helium, including superfluidity (SF). 1 However, unlike helium, bulk hydrogen solidifies at low temperature (T c ≈ 14 K); this solidification has so far prevented the observation of phenomena such as Bose Einstein condensation (BEC) and SF which, in the bulk, are speculated to occur in the liquid phase below T ≈ 6 K. Solidification occurs because of the depth of the attractive well of the potential between two hydrogen molecules, which is significantly greater than that between two helium atoms. Several attempts have been made 2,3,4,5 to supercool bulk liquid p-H 2 , but the search for SF (in the bulk) has so far not met with success.
Toward the aim of maintaining a stable liquid phase of molecular hydrogen down to low enough temperatures for SF to be observed, several avenues can be explored, including reduction of dimensionality. The predicted 6 theoretical value of the melting temperature T m for a solid film of p-H 2 in two dimensions (2D), is slightly less than 7 K, i.e., approximately half that of three-dimensional bulk p-H 2 , but still significantly higher than that at which the system, if it remained a liquid, would undergo Bose Condensation and turn superfluid, estimated at ∼ 2 K in 2D. However, in an actual experimental realization of a quasi-2D system, i.e., a film of hydrogen molecules adsorbed on a substrate, quantum zero-point motion of adsorbed particles in the direction perpendicular to the substrate can be significant, as calculations for adsorbed helium films on alkali metal substrates have shown. 7 One might speculate that such zero-point motion may result in an effective screening of the hydrogen intermolecular interactions, possibly leading to a reduction of T m , with respect to the purely two-dimensional case.
With this aim in mind, extensive experimental studies of hydrogen films adsorbed on a different substrates have been carried out, using various techniques. 8,9,10,11,12,13,14 One of the most remarkable findings, is that the melting temperature T m of a solid p-H 2 monolayer 11 can indeed be significantly less than that of purely 2D p-H 2 . The interesting question is whether by "fine tuning" of the substrate, T m may be lowered to the point where interesting, collective quantum many-body phenomena could become observable, in some thermodynamically stable liquid phase.
In a recent neutron scattering investigation of o-D 2 films adsorbed on a krypton pre-plated graphite substrate, 15 evidence of a stable "liquid-like" phase of o-D 2 down to T ∼ 1.5 K was reported (termed a "domainwall" fluid by the authors). This result is obviously of great interest, as it appears that the substrate considered in this work may indeed provide an avenue to the observation of the long sought superfluid phase of molecular hydrogen, given the low temperature down to which a liquid-like phase may exist. Motivated by this experiment, we have undertaken a theoretical study of the low temperature phase diagram of p-H 2 and o-D 2 films adsorbed on such a substrate.
In a previous study, 16 we determined that p-H 2 , adsorbed upon a substrate consisting of graphite preplated with a commensurate monolayer of krypton forms two thermodynamically stable monolayer phases at low T , both solid; one is commensurate with the krypton layer, the other incommensurate. No evidence was observed of a liquid phase at T = 0 for these or any intermediate coverages, and quantum exchanges of p-H 2 were found to be greatly suppressed. Quantum zero-point motion in the direction perpendicular to the substrate was found to be less significant than that observed for p-H 2 adsorbed on lithium, 17 from which one can conclude that this mechanism ought not yield a significant reduction in the melting temperature. In that study, while the Kr layer was modeled explicitly, the underlying graphite substrate was assumed smooth.
In this paper, we again explore this problem, but with two main modifications with respect to Ref. 16, namely a) we assume this time a compressed krypton layer preplating the graphite, i.e., incommensurate with the underlying substrate lattice structure, and b) we explicitly model corrugation of the underlying graphite substrate, in some of our simulations. Furthermore, in order to allow for a direct comparison with Ref. 15, we study adsorption of both p-H 2 and o-D 2 . Energetic, structural and superfluid properties are investigated using Path Integral ground state (PIGS) Monte Carlo simulations, as well as the continuous-space worm algorithm 18,19 for simulations at finite temperature.
The main results of this study are the following: 1. No evidence is observed of a thermodynamically stable liquid phase in o-D 2 or p-H 2 . In all cases, equilibrium phases are solid monolayers, commensurate with the underlying krypton plating, with coverage (i.e., 2D density) θ • =0.07253Å −2 for both isotopes.
2. An o-D 2 monolayer composed of commensurate domains separated by domain walls is found to be thermodynamically stable, compressible in a relatively extended range of coverages.
3. Quantum exchanges of p-H 2 or o-D 2 molecules are essentially absent in the T → 0 limit; consistently, no evidence of a finite superfluid response of either isotope is ever observed.
4. Our results are consistent with the presence of a "domain-wall fluid", as suggested in Ref. 15. Though no dynamical information can be reliably extracted using the QMC techniques employed here, visual inspection of instantaneous configurations generated by our sampling procedure suggests that such a domain-wall fluid indeed form; however, such a fluid is classical in character, and does not support SF of either hydrogen isotope even at T =0.
The remainder of this manuscript is organized as follows: Sec. II offers a description of the model used for our system of interest, including a discussion of the potentials and the justifications for the main underlying assumptions. Sec. III involves a brief discussion of the computational techniques and specific details of implementation, in addition to details of calibration and optimization. The results are presented in Sec. IV; finally, Sec. V is a summary of the findings and our concluding remarks.

II. MODEL
We consider a system of N hydrogen molecules, regarded as point particles of nuclear spin zero (i.e., obeying Bose statistics; this discussion is entirely equivalent for both isotopes). Molecules sit above a substrate consisting of a single atomic layer of krypton, below which is a graphite substrate. The Kr layer consists of (L) pointlike atoms, pinned at fixed positions R k (k = 1, 2, ..., L) (an assumption justified by their relatively large mass). They are arranged in a perfectly two-dimensional triangular lattice, with a spacing 15 of 3.99Å and sit at a height of 3.46Å over the top layer of C atoms of the graphite substrate; this latter distance corresponds to the energetic minimum of the most accurate Kr-graphite potential available, 20,21 if graphite is regarded as smooth.
If a smooth graphite substrate is assumed, the model quantum many-body Hamiltonian is therefore as follows: Here, m is the mass of a hydrogen molecule, {r i } (with i=1,2,...,N ) are the positions of the hydrogen molecules, r ij ≡ |r i − r j |; z i is the height of the ith hydrogen molecule above the graphite surface. V is the potential describing the interaction between any two hydrogen molecules, and U represents the interaction of a hydrogen molecule with a Kr atom. Finally,Ũ represents the interaction of a hydrogen molecule with the smooth graphite substrate.
As we are investigating p-H 2 and o-D 2 , which have nearly spherically symmetric electron densities, all pair potentials are assumed to depend only on relative distances. The interaction V is described by the Silvera-Goldman potential, 22 which provides an accurate description of energetic and structural properties of condensed p-H 2 and o-D 2 . 23, 24 The interaction of a hydrogen molecule and a Kr atom is modeled using a standard 6-12 Lennard-Jones (LJ) potential; we make use of the Lorentz-Berthelot mixing rule, 25,26 yielding ǫ = 75.6 K and σ = 3.3Å, for our purposes consistent with current state-of-the-art potential energy surfaces (for example, numerical agreement is within ≈ 5% for the well depth). 27 The Kr atomic monolayer is uniform, with each Kr atom kept fixed in its lattice position for the duration of the simulation.
We use a simple "3-9" potential to describe the interaction of hydrogen molecules with the smooth graphite substrate, 28 i.e.,Ũ where C=7913.24Å 3 K and D=259.39 K are parameters derived from the original p-H 2 -C Lennard-Jones parameters 29 (σ = 3.18Å, ǫ = 32.05 K) and the density of carbon atoms in graphite (ρ = 0.114Å −3 ).
The assumption built in (1) of a smooth graphite substrate seems justified, on account of the relatively large (∼ 6Å) distance between the hydrogen molecules and the substrate, due to the presence of the Kr layer. Nevertheless, corrugation may be important in this problem, given that the Kr layer is not registered with the underlying graphite substrate. Therefore, in order to gauge the effect of substrate corrugation, we performed a number of simulations with an explicitly modeled top layer of C atoms, placing a smooth graphite slab underneath. For these simulations, the pair-wise interaction between a hydrogen molecule and a C atom was taken to be the same LJ potential utilized to infer the "3-9" interaction of hydrogen molecules with a smooth substrate.
The results obtained using this composite substrate indicate no qualitative change in the physical character of the system as a result of introducing the smooth plane approximation for graphite. Quantitatively, the potential energy per particle, for example, was found to be ∼2 K lower when assuming explicit corrugation of the graphite substrate, a difference of less than 1%. Henceforth, therefore, we shall confine our discussion to the case of a smooth graphite substrate, i.e., model (1).
A slice of the effective potential energy surface for molecular hydrogen resting at its average height above the substrate is shown in Figure 1; included are the interactions with both the Kr monolayer and the graphite slab. There are two energetically degenerate sublattices of preferential adsorption sites for this geometry -these are denoted sublattice A and sublattice B .
The model (1) clearly contains important physical simplifications, such as the neglect of zero-point motion of Kr atoms, as well as the restriction to additive pairwise interactions (to the exclusion of, for example, three-body terms), all taken to be central, and the use of the highly simplified LJ and "3-9" potentials. Nonetheless, it seems a reasonable starting point, and even quantitatively we expect it to capture the bulk of the physical picture. The two energetically degenerate sublattices of preferential adsorption sites, A and B, are shown. A commensurate coverage is defined as the occupation of one entire sublattice (also termed " 1 2 filling"). This is a low resolution image; contact the authors for high resolution image.

III. COMPUTATIONAL METHOD
Accurate ground state expectation values for quantum many-body systems described by a Hamiltonian such as (1) can be computed by means of Quantum Monte Carlo (QMC) simulations. In this work, the method utilized is Path Integral ground state (PIGS), which is an extension to zero temperature of the standard, Path Integral Monte Carlo method. 30 PIGS is a projection technique, which filters the exact ground state wave function out of an initial trial state. It is therefore closely related to other ground state projection methods, such as Diffusion Monte Carlo (DMC), but has a few distinct advantages (for a discussion, see, for instance, Ref. 31). Because this method is described in a number of publications, it will not be reviewed here. Some of the technical details of the calculation performed in this work (mainly, the short imaginary time propagator) are the same as in Ref. 32.
The trial wave function utilized is of the Jastrow type: The values of the parameters α = 30Å 3 , γ = 250Å 5 and µ = 750Å 5 were obtained empirically, by minimizing the energy expectation value computed in separate variational calculations. Using the trial wavefunction as defined above, we observe convergence of the ground state energy estimates with a projection time 0.250 K −1 , using a time step τ = (ν/2560) K −1 , where ν = 1 for p-H 2 and ν = 2 for o-D 2 .
Finite-temperature results were obtained making use of a reformulation of standard path integral Monte Carlo (PIMC) called the worm algorithm (WA), a scheme which allows for efficient computation of thermodynamic properties, including winding numbers and off-diagonal correlations, for systems of much greater size than that accessible to conventional PIMC simulations. 18,19 The only input parameter to the WA, as for any finite temperature method, is the microscopic model (1). Just as in ground state simulations, an optimal value τ w of the imaginary time step must be determined, such as accurate and unbiased estimates of the observables are obtained; for this system, satisfactory results are yielded by the choice τ w = (ν/320) K −1 .
Calculations for a range of p-H 2 and o-D 2 coverages were carried out, starting from an initial configuration of molecules sitting atop the Kr layer. The simulation cell consists of a 12×14 triangular lattice of Kr atoms with 3.99Å nearest neighbor spacing (resulting in a simulation box of dimension 47.880Å × 48.3762Å). 15 Periodic boundary conditions are used in the three directions, but the simulation cell is chosen sufficiently large in the z direction that they have no effect vertically, because of the strongly attractive character of the composite substrate (i.e., p-H 2 and o-D 2 molecules do not evaporate). It is important to note that to capture the features observed experimentally in Ref. 15 requires that a large number of particles be simulated (the simulation cell must be large enough to support at least two crystalline domains).
The systematic errors of our calculation are attributable to finite projection time (for the T =0 PIGS calculations) and the finite time step τ (for both PIGS and the WA). Based on comparisons of results obtained from simulations with different values of projection time and/or time step, we estimate our combined systematic error on the total energy per p-H 2 molecule to be of the order of 0.7 K or less (corresponding to less 0.5%). The simultaneous use of a ground state and a finite temperature method allows us to obtain an independent check of our calculations in the T → 0 limit, where the two methods must yield the same results, within statistical uncertainties.

IV. RESULTS
Physical quantities of primary interest, for both p-H 2 and o-D 2 , include the energy per molecule, e(θ), and the superfluid fraction, ρ s . PIGS (T =0) was used to compute e(θ) as a function of the coverage (two-dimensional density), θ, from θ=0.06735Å −2 to θ=0.09498Å −2 . For a subset of these coverages, the worm algorithm (WA) was employed to yield finite temperature estimates for all of the above physical quantities, for a set of temperatures between T =5 K and T =0.25 K. The results for e(θ) are shown in Figure 2; note that e(θ) for p-H 2 has been shifted by −50 K for clarity. In the case of p-H 2 , the main feature is an energetic minimum at θ=0.07253Å −2 , corresponding to the formation of a thermodynamically stable solid layer, that is commensurate with the underlying krypton lattice. This coverage is approximately 8% greater than the equilibrium density of p-H 2 in strictly 2D. 6 Within the precision of the calculation, e(θ) is numerically consistent with the absence of other thermodynamically stable coverages, in the range explored in this work.
For o-D 2 , we find that there is an energetic minimum at θ=0.07253Å −2 , again corresponding to a commensurate solid film. In addition, one finds two coverage ranges, separated by a cusp, where θe(θ), shown in Figure 3, has a positive second derivative (meaning that these configurations satisfy the condition of positive compressibility, dµ/dθ=d 2 [θe(θ)]/dθ 2 > 0, and are thus thermodynamically stable); we will call the first of these IC-I, and the second IC-II. IC-I extends from θ=0.07598Å −2 to θ=0.08030Å −2 , while IC-II extends from θ=0.08030Å −2 to θ=0.8635Å −2 . The much richer behavior of e(θ) in the case of o-D 2 , is a direct consequence of the greater mass, and reduced zero-point motion, of this isotope; as a result, the potential energy plays a greater role in shaping the physical properties. Phases of intermediate coverages found for o-D 2 are apparently "washed out" by zero-point motion in the case of p-H 2 . While this result is surprising, similar anomalous behavior is noted in  related studies. 33,34 Of particular interest for the present study is IC-II, where the o-D 2 is found to form a structure best characterized as being composed of several crystalline domains commensurately occupying one of the two triangular sublattices available. Between these commensurate domains are regions of high o-D 2 density -domain walls. Molecules that are part of these domain walls have a decreased average interparticle distance. One may expect that such o-D 2 , having the same zero-point motion as o- In order to provide accurate permutation statistics and superfluid density estimates, additional simulations were carried out for several coverages at finite temperature using the worm algorithm. Figure 2, inset, shows the energy values obtained in the range of temperature 0.25 K ≤ T ≤ 5 K, for θ=0.07253Å −2 and θ=0.08389Å −2 .
Shown in Figure 4 are vertical density profiles, ρ(z), for θ=0.07253Å −2 and θ=0.08389Å −2 at T =1 K. One can see that the vertical extent of the o-D 2 film is slightly enhanced for increased coverages.
While dynamical information cannot be extracted using the computational techniques employed in this study, one does have access to real-space instantaneous configurations of the system, from which one can extract qualitative insight, and potentially use to preclude particular system behaviors. Such representative many-particle configurations are shown in Figure 5 for o-D 2 and p-H 2 with θ=0.07253Å −2 at T = 0, and for o-D 2 with θ=0.08389Å −2 at T =0.5 K. In each of (a) to (c), we find a high degree of particle localization, even along the domain walls observed in (c). No single configuration observed was found to be inconsistent with the instantaneous configurations one would necessarily observe in a domain-wall fluid, while the spatial configuration of domain-walls was not fixed as the system evolved in Monte Carlo time.
Finally, consistent with the high degree of o-D 2 localization observed, permutations are found to not take place at any temperature for any coverage examined, and, accordingly, the finite temperature estimator for the superfluid density consistently yields a value of ρ s =0. In light of the absence of structural and energetic differences found between T =0 and T =0.25 K, it is unlikely that novel phases of o-D 2 (namely, phases with a finite superfluid density) remain unexplored at intermediate temperatures.

V. DISCUSSION
Using numerically exact ground state Quantum Monte Carlo (PIGS) and finite-temperature path integral Monte Carlo (worm algorithm) methods, we studied p-H 2 and o-D 2 adsorption onto a graphite substrate preplated with a compressed krypton monolayer. We performed calculations based a simple model, in which graphite corrugation is ignored, the Kr atoms in the spacer layer are assumed static and point-like, and p-H 2 (o-D 2 )-substrate interactions are given by Lennard-Jones type potentials. Using  an exact ground state technique, for both p-H 2 and o-D 2 , we find that a solid monolayer commensurate with the Kr layer is thermodynamically stable, and, in both cases, is the coverage corresponding to the minimum energy per particle. For o-D 2 , we also find that there are two distinct compressible incommensurate solid regions, with the denser region, IC-II, corresponding to regions of commensuration separated by domain walls. Investigation of this coverage region at finite temperature, using the worm algorithm, down to T =0.25 K, yields no evidence of superfluidity, consistent with our observation of a high degree of localization of o-D 2 molecules. No quantum exchanges, and by extension permutations, are found to take place, though they are sampled efficiently. Though our simulations support the notion of a "domain-wall fluid", as all instantaneous configurations generated by our QMC algorithm do indeed show such domain walls, this mechanism is found not to support superfluidity of either hydrogen isotope.
There are obviously several sources of uncertainty in this calculation which need to be discussed. The potentials used to describe the interactions between the p-H 2 and the substrate are very rough; this does not seem too important an issue, as far as the interaction of p-H 2 molecules with graphite is concerned, given the relatively large average distance at which molecules sit, supported by calculations using explicit corrugation of the graphite substrate. On the other hand, a more realistic interaction potential between p-H 2 and krypton may quantitatively alter the energetics shown here. Despite these issues, and other simplifications, it does not seem likely that the qualitative structural information will change. Thus, we conclude this system is not a good candidate for further consideration in the search for superfluid o-D 2 or p-H 2 .
In general, the results shown here confirm those obtained in a other studies, 35 to the effect that frustration arising from incommensuration does not lead to the formation of a stable superfluid phase of molecular hydrogen at low T . A more complex physical mechanism, causing a weakening of the effective interaction between two hydrogen molecules may have to be investigated, if some progress is to be achieved in this direction. Possibly, the interaction of the hydrogen molecules with the surface electrons of a metal substrate, or of a nanostructure, might have such a "renormalizing" effect.