Two-Dimensional Bose-Hubbard Model for Helium on Graphene
Jiangyong Yu, Ethan Lauricella, Mohamed Elsayed, Kenneth Shepherd Jr., Nathan S. Nichols, Todd Lombardi, Sang Wook Kim, Carlos Wexler, Juan M. Vanegas, Taras Lakoba, Valeri N. Kotov, Adrian Del Maestro
TTwo-Dimensional Bose–Hubbard Model for Helium on Graphene
Jiangyong Yu, Ethan Lauricella, Mohamed Elsayed, Kenneth Shepherd Jr., Nathan S. Nichols,
1, 2
Todd Lombardi, Sang Wook Kim, Carlos Wexler, JuanM. Vanegas, Taras Lakoba, Valeri N. Kotov, and Adrian Del Maestro
5, 6, 1 Department of Physics, University of Vermont, Burlington, VT 05405, USA Materials Science Program, University of Vermont, Burlington, VT 05405,USA Department of Physics and Astronomy, University of Missouri, Columbia, MO 65211, USA Department of Mathematics & Statistics, University of Vermont, Burlington, VT 05405, USA Department of Physics and Astronomy, University of Tennessee, Knoxville, TN 37996, USA Min H. Kao Department of Electrical Engineering and Computer Science,University of Tennessee, Knoxville, TN 37996, USA
An exciting development in the field of correlated systems is the possibility of realizing two-dimensional (2D) phases of quantum matter. For a systems of bosons, an example of strong corre-lations manifesting themselves in a 2D environment is provided by helium adsorbed on graphene.We construct the effective Bose-Hubbard model for this system which involves hard-core bosons( U ≈ ∞ ), repulsive nearest-neighbor ( V >
0) and small attractive ( V (cid:48) <
0) next-nearest neighborinteractions. The mapping onto the Bose-Hubbard model is accomplished by a variety of many-bodytechniques which take into account the strong He-He correlations on the scale of the graphene latticespacing. Unlike the case of dilute ultracold atoms where interactions are effectively point-like, thedetailed microscopic form of the short range electrostatic and long range dispersion interactionsin the helium-graphene system are crucial for the emergent Bose-Hubbard description. The resultplaces the ground state of the first layer of He adsorbed on graphene deep in the commensuratesolid phase with 1 / I. INTRODUCTIONA. Helium on Two-Dimensional Materials: AMany-Body Paradigm
The problem of He atoms deposited on solid sub-strates has been identified for many decades as a bosonicmany-body problem that could exhibit a rich phase dia-gram including the possibility of dimensional crossover[1–8]. Graphite was first recognized as an ideal two-dimensional substrate due to its exceptional homogene-ity, [9] and extensive experimental [10, 11] and theoreticalstudies [12–14] have demonstrated that under the rightcircumstances a superfluid He film can develop on thegraphite surface. Because He atoms are neutral, themany-body interactions that determine the behavior ofthis system are the van der Waals (VDW) interactionsbetween He atom pairs and between He and graphite.Since VDW interactions are typically fairly weak, butlong range, the possibility of superfluidity, and at whichdensity (and film coverage) it can exist depends on theinterplay between the two-body He–He interactions andthe interaction of He with the substrate (in this case car-bon) atoms.Since the discovery of the two-dimensional (2D) ver-sion of graphite, namely graphene [15], the problem ofHe–substrate interactions has been revisited with greatenthusiasm [5, 16–18]. As graphene is a purely 2D sys-tem, the VDW adsorption potential that tends to local-ize helium-4 atoms is 10% weaker (compared to graphitewhich is a bulk material) and therefore there is the ex- otic possibility of purely 2D He superfluidity (atomicwidth film). While graphite’s properties are set by itsbulk structure, graphene’s 2D lattice and (related) elec-tronic structure can be manipulated in a variety of ways.This is the reason why graphene and 2D materials moregenerally have become an attractive area of theoreticaland applied electronics research [19]. For example, dop-ing (addition of electrons or holes into the layer) canbe easily done, or the hexagonal structure can be dis-torted, or hydrogenation agents can be introduced (mak-ing graphene effectively an insulator) [15, 20]. All of theseaffect the graphene lattice and electronic state and, byextension, the VDW potential between He and graphene[16]. Finally, graphene’s dielectric environment can beeasily changed. For example, putting graphene on dif-ferent dielectric substrates immediately affects (screens)the electronic charge resulting in a modified strength ofthe VDW force.For all of the above reasons the problem of He ongraphene, and its extensions, has become a pressing prob-lem due to its potential to produce purely 2D collectivebosonic states. The first question to answer is the behav-ior of helium-4 on pristine graphene in vacuum. So far,theoretical studies [17, 18, 21–23] have concluded thatthe first adsorbed He layer on graphene forms an insu-lating state where He atoms occupy 1/3 of all graphenehexagon centers (energetically preferred location), in atriangular lattice pattern, the so-called C1/3 commensu-rate solid. There is still some limited controversy on thepossible existence of a competing classical or quantumliquid at zero temperature [22] based on the exact form a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b FIG. 1. A schematic of the √ ×√ He atoms (blue) adsorbed on graphene (gold) showing 1 / N s = 72 sites filled where the He atoms are localizedon the sites of a triangular lattice (see shadow). The size ofthe indicated region of graphene is L x × L y (cid:39) (cid:6) A × (cid:6) A. of the helium–graphene potential utilized in simulations,but energy differences are on the order of statistical un-certainties. This pinning of He atoms in this insulatingstate (depicted schematically in Fig. 1) is due to a com-bination of the He–graphene attractive VDW potentialand He–He repulsion, as we discuss below. The secondHe layer can become superfluid [18, 24, 25], as it is far-ther away from the attractive graphene potential, eventhough studies show that a number of other states, in-cluding incommensurate solid phases, are very close to itin energy. Overall, the emergence of superfluidity turnsout to be a very complex many-body problem due to afine balance between fairly weak VDW forces.The aim of this work is to conclusively develop an ef-fective 2D Bose–Hubbard (BH) model for the first layerof He on graphene. The reasons why such a model ishighly desirable are as follows. (1) The results men-tioned above about the existence of the 1/3 insulatingstate are obtained by different zero and finite tempera-ture quantum Monte Carlo (QMC) techniques. In fact,we will complement those with our own version of groundstate continuum QMC. However, to gain intuition aboutthe stability of the C1/3 phase and its proximate phases,it is advantageous to develop an effective lattice Bose–Hubbard model where the most important interactionsare identified. Of course, the phases predicted by the ef-fective BH model must agree with the QMC results. Wewill see that this is indeed the case. (2) It is clear from theoutset that the resulting BH model is highly non-trivialto develop, compared, for example, with BH models usedin cold atom physics (where optical lattice potentials arethe equivalent of the graphene potential here). The rea-son is that in cold atom physics the atom density is verylow (a billion times lower), while in our case of He ongraphene the coverage is high, and atoms are separatedfrom each other on the scale of the graphene lattice (sev-eral Angstroms), which is smaller than the range of theVDW potential. Thus, while interaction effects in coldatom physics are generally easy to incorporate by assum-ing s -wave scattering between atoms [26–29], this is notthe case in our solid-state context where there is a finiterange over which interactions are important. It is not a priori clear that a consistent 2D effective BH modeldescription even exists since the QMC techniques previ-ously mentioned are fully 3D. Thus a careful comparisonbetween “2D restricted” QMC and several other tech-niques has to be made. (3) Finally, armed with such aneffective 2D BH model, one can use it as a first step in theanalysis of a variety of other systems, including situationswhere graphene’s properties are modified (as previouslydescribed), or generalizing to other 2D materials. B. Approach and Summary of Main Results
The main result of this work is that the behavior ofthe first layer of He atoms adsorbed on graphene can becaptured via a single-band “hard-core” Bose–Hubbardmodel with strong (effectively infinite) on-site Hubbardrepulsion ( U ≈ ∞ ). The resulting low energy Hamilto-nian has the form: H BH = − t (cid:88) (cid:104) i,j (cid:105) ( b † i b j + h.c. ) + V (cid:88) (cid:104) i,j (cid:105) n i n j + V (cid:48) (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) n i n j + . . . , (1)where t is the hopping strength, b † i ( b i ) creates (destroys)a bosonic He atom on site i with [ b i , b † j ] = δ i,j , V is thenearest neighbor interaction, and V (cid:48) is the next-nearestneighbor interaction. The ellipsis indicates higher orderinteractions that are neglected here. The sites i, j corre-spond to the vertices of the triangular lattice formed bythe centers of graphene’s hexagons as seen in Fig. 2. Wefind by a variety of methods that: t ∼ V ∼
50 Kand V (cid:48) ∼ − t we employ maximally local-ized Wannier function. First, the VDW potential dueto graphene, acting on a single He atom is calculated bytechniques described in our previous work [16]. Then theone-particle Schrodinger equation in the external VDWpotential is solved numerically and the Wannier func-tions are constructed. The electronic dispersion followsthe symmetry of the triangular lattice and the effectivehopping t from this analysis is ∼ t by independently computing theeffective 2D adsorption potential experienced by a Heatom above the graphene sheet using path integralground state quantum Monte Carlo (QMC) and twotypes of ab initio methods: Density Functional Theory(DFT) [30, 31] improved by including VDW energies inthe appropriate DFT functional [32–34], and 2nd orderMøller–Plesset (MP2) [35–37]. In all cases, the non-interacting band structure resulting from the periodicadsorption potential is determined, where the hoppingcan then be extracted from the bandwidth or overlap in-tegrals. The results are all consistent with the simple ~b ~b t, VV a √ a a ~a ~a FIG. 2. The triangular lattice defined by hexagon centersof the graphene lattice (shown in grey) with the lattice andbasis vectors in Eq. (2) shown in the upper right. The near-est neighbor hopping t and interaction V are experienced be-tween sites separated by √ a (blue) where a (cid:39) . (cid:6) A isthe length of a hexagon side. The dashed red lines indicatethe next-nearest neighbor interaction V (cid:48) between sites sep-arated by 3 a which form a triangular superlattice at 1 / a and a are sharedbetween the graphene and triangular lattice. Wannier theory and provide a powerful confirmation ofthe accuracy of the extracted value of t .Next, we turn to the He–He interaction inducedterms ( U, V, V (cid:48) ) in Eq. (1). The consistent incorpora-tion of interactions proves rather formidable and intro-duces unique technical challenges, not usually encoun-tered when building and estimating parameters in aneffective Bose–Hubbard model applicable to ultra-coldbosonic lattice gases [26, 27, 29, 38, 39]. The He–He po-tential itself between two isolated atoms in vacuum is wellunderstood and can be very accurately parametrized, asa result of decades of study [40–44]. The first observationis that He is special because the s -wave scattering length a s ∼ (cid:6) A [45, 46] is “by chance” (i.e. without any fine-tuning) very large, of the order >
10 times larger than therange of the potential (the effective VDW length). Forvery dilute lattice systems of trappable heavy atoms, theaverage atomic separation (cid:96) (cid:29) a s , and interactions in aneffective BH description can be computed by convolvinghighly localized spatial atomic wavefunctions with nar-row δ -function interactions with a strength proportionalto a s . For the adsorption geometry considered here – heliumatoms confined to move on a triangular lattice with spac-ing a ∼ (cid:6) A due to the proximate graphene structure– we are in the opposite limit where the calculation ofeffective interaction parameters is very sensitive to theshort-range part of the He–He potential. At small scales,below ∼ (cid:6) A, this potential rises rapidly to a very largestrength ∼ K yielding a hard-core description with U ≈ ∞ and promoting the nearest neighbor interaction V to play a dominant role. In addition, the fact thatHe–He is very strong at the lattice scale suggests that aself-consistent formulation has to be employed for the cal-culation of V , which operates on this scale. Calculating V by using the “bare” Wannier functions, i.e. the single-particle localized wavefunction in the field of graphene,produces an unphysically large value V ∼ K. Thisproblem suggests a strategy where a self-consistent ad-justment of the Wannier functions to accommodate thestrong repulsion is employed, for example in the spirit ofthe Jastrow factor commonly introduced in such situa-tions [47–49]. We have determined that instead of work-ing with Jastrow factors, it is more convenient to use theself-consistent Hartree-Fock equations [50, 51]. These areexpected to provide a very accurate description of two-body interactions due to the strongly localized nature ofthe Wannier functions around a given site. We find thatthe Hartree-Fock equations converge to the same result(“fixed point”) which is independent of the details of thepotential at ultra-small distances, producing V ∼
70 K(see § IV D).The value of V can also be calculated within three ad-ditional and complementary approaches. The continuumQMC method mentioned previously, provides a very ac-curate estimate for the adsorbed He wavefunctions andthe total interaction energy at unit filling can be con-verted into an effective V ∼
50 K (see § IV E). This is insatisfactory agreement with the Hartree–Fock method.Van der Waals corrected DFT provides a third, indepen-dent check of the above results which yields V ∼
20 K( § IV F) and ab initio 2nd order Møller–Plesset perturba-tion theory for two adsorbed He atoms on a variety of aro-matic carbons (up to circumcoronene) yields V ∼
50 K.The combination of all the aforementioned techniques(each subject to very different approximations) leads toan effective Bose–Hubbard model (Eq. (1)) with param-eters summarized in Table I. All energies are reportedin K, the natural scale in the adsorption system underconsideration. With the exception of the simple Wan-nier theory (which as discussed above does not properlytake into account the effects of interactions on the latticescale), all results are in good agreement, allowing us todefinitively place helium on graphene within the contextof the extended hard-core Bose–Hubbard model on thetriangular lattice.
Method t / (K)
V / (K) V (cid:48) / (K) t/V Wannier 1.45 7540 638 0.0002HF 1.45 69.7 -2.08 0.021QMC 1.38 54.3(1) -2.76(2) 0.025DFT 1.10 21.4 -1.36 0.051MP2 0.59 51.5 -1.97 0.011TABLE I. The hopping parameter t , nearest and nextnearest-neighbor interaction V and V (cid:48) , and the ratio of t/V of the effective Bose–Hubbard model defined in Eq. (1) ascalculated by the five different methods: Wannier functions,Hartree–Fock (HF), quantum Monte Carlo (QMC), DensityFunctional Theory (DFT), and Møller–Plesset perturbationtheory (MP2). In all cases, t is calculated via the bandstructure of single helium atom subject to a periodic two-dimensional adsorption potential V He − s . Note that t is thesame for Wannier and Hartree–Fock as they use the sameempirical potential. C. Implications for the Quantum Phase Diagram
The phase diagram of Eq. (1) in the limit of infinite U and considering only nearest neighbor interactions ( t − V model) can be analyzed within the mean-field theory[52, 53], as shown in Fig. 3. This result is known to be inqualitative agreement with lattice quantum Monte Carlofor hard-core bosons with extended interactions [54–56].For small values of the chemical potential (low fillingfraction) three phases are identified: the C1/3 phase dis-cussed previously, a supersolid phase, and uniform su-perfluid phase. The phase boundary between the solid f = 1 / f = 1 / t/V ∼ /
50 (Table I), places Heon pristine graphene firmly in the C1/3 phase, (as shownby the symbols) consistent with previous simulations ofthe full three dimensional system [18].
D. Paper Outline
In the remainder of the paper, we provide a discus-sion of the microscopic models we employ to characterizea three-dimensional system of helium atoms interactingwith a two-dimensional graphene membrane. We thendiscuss in what context or limits this system can be un-derstood within an effective 2D theory. Working withinthe 2D limit, we provide details of the approaches brieflydiscussed in the introduction to estimate the parametersof a hard-core extended Bose–Hubbard model. This in- .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . t/V µ / V superfluidsupersolidsolid f = 1 / f = 2 / FIG. 3. The mean-field phase diagram for hard-core bosonson the triangular lattice with nearest-neighbor interactions V and hopping t with density controlled by the chemical po-tential µ . Identified phases include commensurate solids (atfillings f = 1 / , / He on graphene resides deep in the commensuratesolid phase at 1 / µ/V has the same valueas the tip of the first lobe. cludes studying the band structure of a single He atomadsorbed on graphene, and determining wavefunctionsfor the many-particle system at different levels of sophis-tication. Finally, we conclude by comparing all our re-sults and describe the exciting future directions this workopens up for studying hard-core Bose–Hubbard modelsin a solid state setting.All data and code needed to generate the results in thispaper are available online [57].
II. MODEL: HELIUM ON GRAPHENE
We consider a system of He atoms in proximity to agraphene substrate frozen at z = 0 with lattice ( a ) andbasis ( b ) vectors: a = a (cid:16) √ , (cid:17) , b = a (cid:16) √ , (cid:17) a = a (cid:16) −√ , (cid:17) , b = a (0 ,
1) (2)where a (cid:39) . (cid:6) A is the carbon–carbon distance as de-picted in Figs. 1 and 2. In this paper we consider twoclasses of methods distinguished by how interactions arehandled. (A) In tight binding, Hartree–Fock, and quan-tum Monte Carlo calculations we employ empirical inter-action potentials, while (B) van der Waals corrected den-sity functional theory [30, 31] and Møller–Plesset [35] per-turbation theory utilizes an ab initio estimate for the in-teraction energy within the Born–Oppenheimer approxi-mation. The combination of these two classes of methodsensures a broader regime of applicability and improvedconfidence in our final results for the mapping of the mi-croscopic system to an effective extended Bose–Hubbardmodel.
A. Empirical
A system of N He atoms of mass m interacting withthe graphene membrane can be described in first quanti-zation via the Hamiltonian: H = − (cid:126) m N (cid:88) i =1 ∇ i + N (cid:88) i =1 V He − s ( r i ) + (cid:88) i 961 K and σ = 2 . (cid:6) A,which are different from previous studies that employedparameters determined for graphite [17, 18, 25, 66]. Theresulting empirical potential is shown in Fig. 4(b) which . 25 2 . 50 2 . 75 3 . 00 3 . 25 3 . 50 3 . r / (cid:6) A − − V H e − H e / K (a) r z / (cid:6) A − − V H e − s / K (b) z empiricalDFTMP2 FIG. 4. The interaction (a) and adsorption (b) potentialsin Eq. (3). The three curves in both panels indicate threedifferent approaches to the potentials. Empirical indicatesthe formulas discussed in the text where V He − He is taken fromRef. [43] while V He − s is determined from Eq. (4) at r = (0 , § IV F and IV G. has a minimum at the center of a graphene hexagon adistance z min (cid:39) . (cid:6) A above the membrane with depth ∼ − 190 K. Since the smallest reciprocal lattice vector, | g | (cid:39) (cid:6) A − , | gz | ≈ (cid:29) g converges rapidly such that onlya few sets with equal | g | need to be retained. B. Ab Initio Here we briefly discuss the main conceptual differ-ences between the “empirical” approach outlined above,where the van der Waals interactions are used to calcu-late physical quantities via many-body techniques (suchas the Hartree–Fock method and quantum Monte Carlo),and ab initio methods. In the latter, calculations areperformed using the usual Born–Oppenheimer approxi-mation where all atomic nuclei are considered classicaland only the electrons receive a full quantum-mechanicaltreatment. Using van der Waals corrected Density Func-tional Theory and 2nd order Møller–Plesset perturbationtheory, we have computed the effective interactions be-tween He–He and He–graphene with the results shown inFig. 4 with the details included in § IV F-IV G.For interactions between two helium atoms, MP2agrees very well with the empirical potential [43, 44] (asseen in Fig. 4(a)) motivating the choice of basis set em-ployed. The dispersion corrected DFT predicts He–Hehard-core interactions that are weaker than V He − He atshort distances but correctly captures the location of theminimum r min (cid:39) . (cid:6) A. The ab initio computation ofthe height ( z ) dependence of the adsorption potential ata fixed position in the xy plane corresponding to the cen-ter of a graphene hexagon, (as seen in Fig. 4(b)) yield avalue of z min (cid:39) . − (cid:6) A where the minima is observedwith a depth varying between -400 and − 300 K.The agreement between the He–He and He–graphenepotentials is remarkable, in light of the drastically differ-ent approximations at play (e.g. frozen nuclei vs. dis-persion) and the large variance in V He − s reported in theliterature for various ab initio approaches [67].All adsorption potentials lead to the existence of a well-defined monolayer of He on graphene. Even with thedifferences between the interactions on display in Fig. 4,the resulting effective 2D low energy model that describesthe system will turn out to be remarkably similar. III. DIMENSIONALITY OF THE FIRSTADSORBED LAYER Regardless of the form of the employed interaction po-tential in the microscopic model, the goal of this workis to obtain access to properties of the ground stateof the N -particle three-dimensional time-independentSchr¨odinger equation: H Ψ ( R ) = E Ψ ( R ) (6)in order to determine the parameters of an effectivetwo-dimensional Bose–Hubbard Hamiltonian describedby Eq. (1) where R ≡ { r , . . . , r N } are the spatial lo-cations of helium atoms.The basic physical picture of adsorption of helium ongraphene is clear. At low temperature and densities,atoms preferentially adsorb to the strong binding siteslocated at the center of graphene hexagons due to theattractive interaction seen in Fig. 4(b). If the density islow enough that interactions between helium atoms arenot relevant, Eq. (3) can be numerically integrated toobtain the z -dependence of the wavefunction in the ap-proximation where the corrugation is neglected and theadatoms experience an average smooth potential over the xy -plane ( i.e. taking only the g = 0 term in Eq. (4), seeAppendix A for details). The resulting single particledensity in the z -direction [17, 68, 69] is shown in Fig. 5along with values corresponding to the adsorption poten-tials computed via ab initio methods. Thus, single atomsare strongly localized around z ≈ (cid:6) A, regardless of the z / (cid:6) A0 . . . . . . ρ ( z ) / (cid:6) A − max z / ˚A ∆ z rms / ˚A ∆E / Kempirical 2.69 0.161 21.2DFT 3.07 0.149 23.5MP2 2.93 0.142 25.5 empiricalDFTMP2 FIG. 5. The particle density, ρ ( z ) ∝ | φ ( z ) | in Eq. (A5)obtained via the shooting method for V He − s computed forthe three different approaches to the potentials described in § II A-II B. The table indicates the value of z at which themaximum density occurs, the root mean squared value of z ∆ z rms = (cid:113) (cid:104) z (cid:105) − (cid:104) z (cid:105) , and the difference between the zeropoint energy and the minimum of V He − s for each poten-tial. All three methods yield density profiles for the adsorbedlayer that are effectively two dimensional with sub-angstromwidths. way the adsorption potential is calculated, with an rootmean squared width of ≈ . (cid:6) A and a zero point energythat lifts the ground state ≈ 20 K above the classical po-tential minimum.As the density of adatoms is increased, there is nowa competition between the energy gained due to attrac-tion of the graphene sheet, and the interaction poten-tial between helium atoms, V He − He , which has an at-tractive minimum at r min ≈ (cid:6) A and eventually becomesrepulsive at smaller distances (see Fig. 4). The lengthscales defining V He − He should be compared with thoseimposed by the graphene corrugation potential where thenearest neighbor distance between two hexagon centersis r NN = √ a (cid:39) . (cid:6) A while the next-nearest neigh-bor distance, corresponding to one out of every threehexagons occupied has r NNN = 3 a (cid:39) . (cid:6) A as seenin Fig. 2. Thus at low densities, the system stabilizesat a single well-defined 2D monolayer, that can exist atboth commensurate and incommensurate filling fractions f = N/N s (where N s is the number of triangular latticesites) in a regime of coverage where both the adsorptionand interaction energies are attractive. As the densitycontinues to increase, eventually the cost of repulsive in-teractions between helium atoms overcomes the reducedattraction felt further from the sheet and layer comple-tion is reached near f ≈ . 6. At this point, a secondlayer begins to form and the system can no longer beconsidered as effectively two dimensional (see Fig. 13 in § IV E 2).This simple picture has been validated by 50 yearsof experiments [1, 3, 70–74] and numerical simulations[12–14, 21, 22, 65, 69, 75, 76] on helium adsorbed ongraphite, where the adsorption potential is 10% strongerthan graphene. While no experiments yet exist in thegraphene system considered here, quantum Monte Carlosimulations [17, 18, 22, 24, 25, 66, 77] both at zero andfinite temperature show analogous behavior. As alreadydiscussed in the introduction, in the first layer, a com-mensurate √ × √ ◦ incompressible C1/3 solidphase (helium atoms occupy 1/3 of the strong bindingsites on a triangular lattice (hexagon centers) with con-stant √ a and axes rotated by 30 ° with respect to theoriginal graphene triangular lattice) is thermodynami-cally stable over a large range of chemical potentials [18](see Figs. 1 and 2) and may compete with a lower densityliquid [66] depending on simulation details and the em-ployed form of V He − s . All observed phases in the firstlayer are incompressible, with no systematic evidence offinite superfluid density surviving extrapolation to thethermodynamic limit.Thus, the ground state of the first adsorbed mono-layer of He on graphene can be described by an ef-fective two-dimensional system. We now discuss how itcan be mapped at low energies onto a single-band Bose–Hubbard model, which requires moving beyond the sim-ple continuum one-body model described here and un-derstanding the role of interactions in Eq. (3). IV. EFFECTIVE 2D BOSE–HUBBARDDESCRIPTION We attack this problem at various levels of sophistica-tion starting from the non-interacting band structure andWannier theory (where we analyze the corrugation of theadsorption potential) and systematically explore the ef-fects of interactions in different approximation schemes:Hartree, Hartree–Fock, quantum Monte Carlo, Møller–Plesset and dispersion corrected density functional the-ory.In this section we introduce an effective 2D V He − s ( r )potential where r = ( x, y ) is the in-plane coordinate.The way V He − s ( r ) is determined depends on the spe-cific method used and will be discussed on a case by casebasis. A. Mapping onto a Bose–Hubbard Model First, we briefly outline the well-known general proce-dure for mapping the interacting problem in Eq. (3), ontothe effective Bose–Hubbard model Eq. (1). This mappingis valid at low energies and therefore the two representa-tions lead to the same ground state properties. A similarmapping has been used to analyze the properties of di-lute Bose gases confined on optical lattices [26, 28, 29];however the physics in our case turns out to be funda-mentally different due to the importance of short-range correlations for a (fairly dense) collection of helium atomsconfined to the graphene lattice.We begin by expressing the first-quantized microscopicHamiltonian in Eq. (3) in second quantization for asingle 2D monolayer, via the introduction of bosonicfield operators, ˆΨ( r ) , ˆΨ † ( r ) such that the local densityis n ( r ) = ˆΨ † ( r ) ˆΨ( r ). In this notation, the effective 2DHamiltonian can be written as a sum of a one-particleterm, which includes the kinetic energy and the helium–graphene potential, and a two-body term (that originatesfrom the helium–helium interaction) [78, 79]: H = (cid:90) d r ˆΨ † ( r ) (cid:18) − (cid:126) m ∇ r + V He − s ( r ) (cid:19) ˆΨ( r ) +12 (cid:90) (cid:90) d r d r (cid:48) ˆΨ † ( r ) ˆΨ † ( r (cid:48) ) V He − He ( r − r (cid:48) ) ˆΨ( r (cid:48) ) ˆΨ( r ) , (7)where a discussion of V He − He is included in Appendix A.For helium atoms strongly confined near 2D triangularlattice locations r i defined by the centers of graphenehexagons (see Fig. 2), the field operators can be expandedover a complete orthonormal set of localized Wannierfunctions ψ ( r − r i ) and the bosonic annihilation and cre-ation operators b i , b † i [78, 79]:ˆΨ( r ) = (cid:88) r i ψ ( r − r i ) b i , ˆΨ † ( r ) = (cid:88) r i ψ ∗ ( r − r i ) b † i . (8)We use the shorthand notation b † i ≡ b † r i for an operatorthat creates a boson at r i , and ψ i ( r ) = ψ ( r − r i ) for theWannier function localized around the site i on a trian-gular lattice. The Wannier functions will be constructedin the next section, and we assume that they correspondto the lowest energy band.Substituting Eq. (8) into Eq. (7) one obtains the effec-tive lattice Bose–Hubbard Hamiltonian in Eq. (1), wherethe boson density is n i = b † i b i . The one-particle hop-ping ( t ) and the density–density interactions on-site ( U ),at nearest-neighbor sites ( V ), and next-nearest-neighborsites ( V (cid:48) ) on a triangular lattice are then given by thestandard expressions [79]: t = − (cid:90) d r ψ ∗ ( r ) (cid:20) − (cid:126) m ∇ r + V He − s ( r ) (cid:21) ψ ( r − a ) (9) U = (cid:90) (cid:90) d r d r (cid:48) | ψ ( r ) | V He − He ( r − r (cid:48) ) | ψ ( r (cid:48) ) | (10) V = (cid:90) (cid:90) d r d r (cid:48) | ψ ( r ) | V He − He ( r − r (cid:48) ) | ψ ( r (cid:48) − a ) | (11) V (cid:48) = (cid:90) (cid:90) d r d r (cid:48) | ψ ( r ) | V He − He ( r − r (cid:48) ) | ψ ( r (cid:48) − a − a ) | . (12)Here, the choice of lattice site for the computation isarbitrary due to translational invariance, and one canreplace a ↔ a . B. Band Structure and Effective Hopping t In order to calculate the overlap integrals in Eqs. (9)–(12), we start by evaluating the general band structureand specifically the hopping parameter t which are deter-mined by the purely one-particle Hamiltonian, − (cid:126) m ∇ r + V He − s ( r ). For a given effective 2D potential V He − s ( r ),the procedure is straightforward [50]. Bloch’s theoremstates that the solutions to the Schr¨odinger equation in aperiodic potential are the product of a periodic function u ( n ) k ( r ) and a plane wave, Ψ ( n ) k ( r ) = e i k · r u ( n ) k ( r ), where k is the 2D lattice quasi-momentum. The index n labels thedifferent bands with corresponding energies ε ( n ) ( k ). Weseek wave functions for the lowest energy band ( n = 1)and hence omit the band index for simplicity. Once theBloch wave-functions are found, the localized Wannierfunctions are constructed via ψ ( r − r i ) = 1 √ N s (cid:88) k ∈ BZ e − i k · r i Ψ k ( r ) , (13)where the summation is over the first Brillouin zone, and N s is the number of (triangular) lattice sites. Eq. (13)can now be used in the overlap integral for t defined inEq. (9) for a given value of V He − s computed within theempirical or ab initio approach. 1. Empirical Here the bare potential is given by Eq. (4) and weuse two approaches to construct an effective 2D potential V He − s ( r ).Following the discussion in § III, we can integrate thefull 3D helium–graphene interaction potential over theprobability density in the z -direction presented in Fig. 5as described in detail in Appendix A. This leads to a 2Dpotential (cid:101) V He − s ( r ) as defined in Eq. (A4c). The corre-sponding band structure is presented in Fig. 6 and theresulting Wannier function is plotted in Fig. 7(a). Basedon these results, Eq. (9) is evaluated in the lowest bandresulting in: t W = 1 . 45 K . (14)In addition, it is straightforward to show that in thetight-binding approximation the lowest band is describedby the explicit formula: ε ( k ) − ε = − t (cid:34) cos( k x a ) + 2 cos (cid:18) k x a (cid:19) cos (cid:32) √ k y a (cid:33)(cid:35) , (15)where a = √ a and ε is an energy offset. This meansthat the bandwidth, defined as the energy difference be-tween the K point (located at momentum (4 π/ a, , 0) in Fig. 6, is equal to 9 t . Equation (15)is plotted as the dashed line in Fig. 6, and the consider-able agreement provides further validation for mappingfrom the continuum to a lattice model. FIG. 6. The band structure obtained using (cid:101) V He − s ( r ) as in-troduced Eq. (A4c) along a high symmetry path in the firstBrillouin zone as shown in the inset. The first band (with n = 1) is well separated from the higher excited bands andthus the low energy properties of the system are determinedby the lowest band. The dashed line shows the tight bind-ing dispersion from Eq. (15), in excellent agreement with thecontinuum model supporting the use of an effective 2D lat-tice model. The shaded region corresponds to the bandwidthequal to 9 t by symmetry. An alternative approach to obtaining a 2D effectivepotential is to exactly simulate a single He atom sub-ject to the full 3D potential via quantum Monte Carlo asdescribed in detail in § IV E and obtain the adsorptionpotential as a function of the 2D coordinate in the plane, r , via Eq. (26): V He − s ( r ) = (cid:104)V He − s ( x, y ) (cid:105) . The corre-sponding hopping parameter calculated from this poten-tial is t QMC = 1 . . (16)where the parenthesis indicates the statistical uncertaintyin the last digit. We note that this value agrees withthat computed using the adsorption potential determinedfrom the 1D wavefunction in Eq. (14) at the order of 10%. 2. Ab Initio The hopping parameter t can also be estimated foran effective 2D potential computed within the ab initioapproximation. While it is computationally difficult toperform a DFT and MP2 calculation for every position r , these numerical methods can readily determine theadsorption potential at the high symmetry points corre-sponding to the minima, maxima, and saddle point (asshown in Fig. 8). Since the summation over | g | is dom-inated by the terms with the smallest magnitudes andconverges rapidly, the full 2D potential can be approxi-mated as V He − s ( r ) = V + c g (cid:88) | g | = g e i g · r + c g (cid:88) | g | = g e i g · r , (17) FIG. 7. Spatial dependence of the density ρ ( x, y ) = | ψ ( x, y ) | /N of an adsorbed He atom on graphene for a y = 0cut in the xy -plane corresponding to the lattice path shownin the top inset. (a) A comparison of the localized Wannierfunction defined in Eq. (13) to that computed via quantumMonte Carlo (QMC) for a single particle N = 1, f = 1 /N s byslightly biasing a single site at the level of the trial wavefunc-tion as discussed in § IV E. The wavefunction strongly pene-trates into neighboring lattice sites, leading to the breakdownof Wannier theory for the computation of interaction parame-ters. (b) The density at unit filling N = N s as computed viaQMC, and within the Hartree-Fock and Hartree approxima-tions showing the tendency towards exponential localizationon a single site. In both panels, the 2D normalization is com-puted over the full graphene sheet. where g = 4 π/ (3 a ) and g = √ g are the lengths ofthe two smallest set of g vectors . The coefficients c g , can be uniquely determined from the minimum, maxi-mum, and saddle point values of the potential as c g = − 19 ( V maxHe − s − V minHe − s ) ,c g = 18 (cid:0) V spHe − s − V minHe − s (cid:1) − (cid:0) V maxHe − s − V minHe − s (cid:1) . (18)A summary of the relevant parameters calculated withdifferent methods is presented in Table II.The calculation of the hopping parameter t then pro-ceeds as in the previous section, where the Wannier func-tions are determined using the 2D potential in Eq. (17),leading to: t DFT = 1 . 10 K , t MP2 = 0 . 59 K . (19) yx V minHe s V maxHe s V spHe s FIG. 8. The effective 2D potential V He − s ( x, y ) used to calcu-late the hopping t can be reconstructed from MP2 and DFTcalculations by determining three values corresponding to theminimum, maximum, and saddle-point values as indicated.The resulting scale along the black line can be seen in Fig. 18.Method (cid:0) V maxHe − s − V minHe − s (cid:1) / (K) (cid:0) V spHe − s − V minHe − s (cid:1) / (K)HF 21.2 17.5QMC 24.7 21.7DFT 39.2 36.1MP2 72.2 66.0TABLE II. The parameters taken from the adsorption poten-tial for the four different methods at the high symmetry pointscorresponding to the minima, maxima, and saddle point re-quired to calculate the coefficients c g , in Eq. (17) (HF andWanner use the same potential). The results for t from the QMC, and DFT adsorptionpotentials are remarkably similar given the variation intheir underlying approximations to the full 3D system.MP2, on the other hand, predicts a smaller value of t , asa result of the significantly stronger adsorption potential V He − s ( r ) from this method.As a final check on the physical realism of these re-sults, the WKB method can be used to estimate t asdiscussed in Appendix C, leading to results which are invery reasonable agreement with those presented above. C. Interaction Effects: Breakdown of the WannierTheory So far, we have been considering the mapping of the2D adsorbed He layer within the single particle approx-0imation. Now, we proceed with an evaluation of the in-teraction parameters U, V, V (cid:48) . In any parameterizationof the He–He interaction potential, the existence of astrong hard-core will preclude the double occupation ofa single site on the triangular lattice and thus effectively U = ∞ as it is the dominant scale with U (cid:29) t, V, V (cid:48) . Forthe potentials in Fig. 15(a) we find that Eq. (10) yields U > K. Therefore, the effective Bose–Hubbard modeldescribes hard-core bosons hopping on the triangular lat-tice formed by the graphene hexagon centers (Fig. 2). Us-ing the single-particle Wannier function approach, onecan also compute the nearest neighbor ( V ) and next-nearest neighbor ( V (cid:48) ) parameters directly from Eqs. (11)and (12) which lead to V W = 7540 K , V (cid:48) W = 638 K . (20)The resulting enormous energy scales associated withthese parameters are unphysical and suggest that thespatial extent of one-particle wave function is too large,and fails to capture the correct interaction physics. Thiscatastrophe originates from the fact that we study theadsorption of He atoms on a solid-state substrate andconsequently both the spatial extent of the one-particlewavefunction (determined by the graphene lattice struc-ture), and the most prominent (repulsive) part of theHe–He potential, vary on the same length scale, of orderof several ˚A.This is in contrast to the case of dilute cold atomicgases confined in optical lattices, where the confinementscale, as well the optical lattice wavelength (“lattice spac-ing”), are typically on the order of a µ m. These scales aremuch larger than the range of the interaction potential( ∼ (cid:6) A) such that the two-body interaction is contact-like and is assumed to be described by a δ -function.Thus in cold atom systems, which provide the mainexamples of Bose–Hubbard models and associated quan-tum phase transitions in nature, the one-body confine-ment scale and the two-body interaction scale are wellseparated. Due to this, the effective model is of the t − U type [28], with a finite Hubbard U and irrelevant (i.e.much smaller) additional interactions V, V (cid:48) .In the case where the 2D limit is reached via adsorp-tion on graphene, the very strong He–He repulsion onthe scale of the one-particle wave-function effectively pro-duces an infinite on-site Hubbard U and therefore it is thenearest-neighbor V and next-nearest neighbor V (cid:48) thatdetermine the relevant quantum phases of the system,leading to the hard-core t − V − V (cid:48) model consideredhere. Therefore, the determination of V, V (cid:48) presents con-siderable technical challenges and has to be done via so-phisticated techniques that take into account the cor-rect structure of the wave function which is modifiedby two-body interactions and at finite density deviatessignificantly from the one-particle results presented sofar. In this sense, our analysis is quite unconventionalcompared to the usual approaches to the Bose–Hubbardmodel. Because of the well-localized structure of themany-body wave functions (as will be clear from the re- sults of the next sections), the effective Bose–Hubbardmodel is still dominated by two-body (density–density)interactions, with the nearest-neighbor term being thelargest one ( V (cid:29) | V (cid:48) | ).The remainder of this section presents a number of dif-ferent approaches to gain access to the many-body wave-functions of He on graphene in order to compute V and V (cid:48) exemplifying the strong correlations in the problem. D. Hartree–Fock Approach to InteractionParameters Here we provide details on how the parameters V and V (cid:48) of the effective Bose–Hubbard model can be computedfrom an effective 2D model of the adsorbed layer. Since V is the energy of nearest-neighbor interaction, then forits computation, one needs to consider a helium layerwith a unit filling fraction. However, as noted in § III, He atoms at this density form two, not one, layers ontop of graphene. To resolve this issue, we will rely onan important result from our QMC simulations, which isdescribed in detail in § IV E. Namely, a quasi-2D, single-layer arrangement of helium over graphene is restoredwhen one imposes a confining potential in the z -direction.Importantly, the particle density in that direction ob-tained with the confinement is close to the density pro-file at filling fraction f = 1 / 3; see e.g. Figs. 12 and 13.This justifies the use of a 2D model for the approximatecomputation of nearest-neighbor He–He interactions.Let us stress again that the spatial extent of the max-imally localized Wannier functions found in the previ-ous subsection (well-suited for the description of an iso-lated helium atom), is on the order of the spacing be-tween the nearest graphene hexagon centers. Therefore,the standard approach of computing interaction param-eters in the Bose–Hubbard model via the overlap inte-gral Eq. (11) would give an unphysically large result.However, we note that the mutual repulsion of adjacenthelium atoms narrows their wavefunctions considerablycompared to the Wannier functions (see Fig. 7). We willnow show how these narrower wavefunctions are foundand then use them in the calculation of V and V (cid:48) viaEqs. (11) and (12).Such wavefunctions are obtained by numerically solv-ing a system of 2D Hartree–Fock equations [51]: − (cid:126) m ∇ r ψ i ( r ) + V He − s ( r ) ψ i ( r )+ (cid:88) i (cid:54) = j (cid:90) d r (cid:48) ψ ∗ j ( r (cid:48) ) V He − He ( r − r (cid:48) ) × [ ψ j ( r (cid:48) ) ψ i ( r ) + ψ i ( r (cid:48) ) ψ j ( r )] = ˜ E i ψ i ( r ) , (21)where the wavefunctions ψ i ( r ) ≡ ψ ( r − r i ) also satisfy theorthonormality constraint: (cid:90) d r ψ i ( r ) ψ j ( r ) = δ ij , (22)1with δ ij being the Kronecker delta.Equations (21), (22), were solved by the acceleratedimaginary-time evolution method (a variant of fixed-point iterations), whose general framework for systemsof equations subject to constraints were laid out in [80];its technical details will be described elsewhere. To es-timate the significance of the exchange interaction (i.e.,the last term, ψ i ( r (cid:48) ) ψ j ( r ), in Eq. (21)), we also simulatedthe Hartree approximation, obtained from Eq. (21) bydropping that term and not imposing the constraint inEq. (22).We performed simulations for the z -averaged poten-tials V He − s and V He − He , as described in Appendix A. Forthe potentials averaged with two different ρ ( z )’s: that de-fined in Appendix A and that found by QMC ( § IV E),Eq. (11) gives, respectively: V HF = 69 . . V He − He isreduced (smoothened) more by the more spread-out ρ ( z )obtained by the QMC. On the other hand, the contribu-tion of the difference between the two averaged V He − s ’sto the difference in the corresponding V ’s is negligible.In fact, we found that the effect of even larger — on theorder of 50% — changes in the magnitude of V He − s on V was well under 1%. For completeness, we also note thatwhen we used V He − s and V He − He averaged with ρ ( z ) de-fined in Appendix A but used the Hartree rather thanHartree–Fock approximation, we found V = 72 . V (cid:48) computed from Eq. (12) by anyof these approximations equals − . V HF = 69 . , V (cid:48) HF = − . 08 K . (23)For a system of localized bosons with strong short-rangeinteractions, the Hartree–Fock equations provide a veryaccurate description as many-particle correlations be-yond the scope of the method are expected to be weak. Inaddition, and quite reassuringly, we find that the aboveresults are similar to those obtained by the accuratemany-body quantum Monte Carlo technique. E. Quantum Monte Carlo At T = 0 K, the path integral ground state quantumMonte Carlo (QMC) algorithm [81–83] provides access toground state properties of a many-body system by sta-tistically sampling the imaginary time propagator e − βH .Starting from a trial wave function | Ψ T (cid:105) , in the longimaginary time limit β → ∞ , e − βH | Ψ T (cid:105) converges tothe exact ground state, | Ψ (cid:105) , provided (cid:104) Ψ | Ψ T (cid:105) (cid:54) = 0.Within this framework we can directly compute ground state properties by statistically sampling the expectationvalue of an observable O , via: (cid:104)O(cid:105) (cid:39) (cid:104) Ψ T | e − βH O e − βH | Ψ T (cid:105)(cid:104) Ψ T | e − βH | Ψ T (cid:105) . (24)We work in a first-quantized representation | R (cid:105) in 3 spa-tial dimensions where configurations are sampled fromthe 3 + 1 dimensional imaginary time worldlines of inter-acting particles. Appendix B provides additional detailson the convergence and scaling of our QMC approachand the source code can be found online [84].In the remainder of this subsection we discuss howQMC simulations of the 3D microscopic many-bodyHamiltonian in Eq. (3) can be analyzed in the contextof an emergent 2D Bose–Hubbard model. We begin byconfirming the single-particle description of the adsorbedmonolayer described in Section III which allows us tocompute an effective 2D potential that can be used to de-termine the hopping parameters t . We then proceed byreducing the size of the simulation cell in the z -directionwhere the extra dimensional confinement allows us to sta-bilize a monolayer at the large filling fractions needed todetermine the interaction parameters V and V (cid:48) . 1. Single Particle Properties: f = 1 /N s We begin with the simplest case of considering a single He atom proximate to the graphene surface at T = 0.The results of QMC simulations are shown in Figure 9for N = 1 with N s = 24 adsorption sites that arecommensurate in a cell with volume L x × L y × L z =9 . (cid:6) A × . (cid:6) A × . (cid:6) A = 1257 (cid:6) A . The cell has periodicboundary conditions in the x and y directions, while mo-tion in the z -direction is restricted through the graphenesheet at z = 0 and a hard wall at z = L z enforced by thepotential V wall ( z ) = V He − s ( r min )1 + e ( L z − r vdW − z )) / ∆ . (25)Here, r min = a ( √ / , / , 1) is located at z = a abovea carbon atom such that V He − s ( r min ) ∼ O(10 ) K setsthe scale of the repulsive potential, r vdW ≈ . (cid:6) A is thevan der Waals radius of helium, while ∆ = 0 . (cid:6) A de-fines the rapidness of its onset. The functional form ofEq. (25) and the choice of parameters are unimportantat filling fractions f (cid:46) / L z (cid:38) (cid:6) A. For thevalue L z = 10 (cid:6) A considered here, simulation results areindependent of L z and can be considered to be reflectiveof bulk adsorption phenomena.Figure 9(a) shows the particle density in the z -direction determined from the expectation value ρ ( z ) = (cid:68)(cid:80) Ni =1 δ ( z i − z ) (cid:69) ∝ (cid:82)(cid:82) d x d y | Ψ ( x, y, z ) | via Eq. (24)( N = 1 here). It has a well-defined peak near 2 . (cid:6) A anda corresponding sub- (cid:6) A width (shown as the full widthhalf maximum) demonstrating that adsorbed He atoms2 z / (cid:6) A0 . . . ρ ( z ) / (cid:6) A − (a) N = 1 , N s = 24 − y / ˚A ρ ( x, y ) (b) − hV He − s ( x, y ) i (c) − . . . x / ˚A0 . . ρ ( x , ) / (cid:6) A − (d) − . . . x / ˚A − − − h V H e − s ( x , ) i / K (e)lowhigh FIG. 9. Quantum Monte Carlo results for a single He atom( N = 1) above a graphene membrane with N s = 24 adsorp-tion sites. (a) The average density in the z -direction showinga well defined particle position a distance 2 . (cid:6) A above thesheet with a width of 0 . (cid:6) A. (b) The average density in the xy -plane showing the ability of a single particle to hop be-tween the sites of the triangular lattice. White dots showthe location of carbon atoms (not to scale). (c) The aver-age potential energy experienced by the He atom due to thegraphene sheet in the xy -plane. (d) A horizontal cut of theparticle density ρ ( x, y ) along the line y = 0. (e) A horizon-tal cut of the adsorption potential V He − s ( x, y ) along the line y = 0. indeed form a quasi-two dimensional layer. Panel (b) in-cludes the average particle density in the plane normal-ized such that N = (cid:82)(cid:82) d x d y ρ ( x, y ) and the existenceof density in each of the N s = 24 adsorption sites isevidence of particle hopping and an ergodic simulation.The lower panel (d) is a cut showing the scale of den-sity fluctuations. Panel (c) shows the average adsorptionpotential experienced by the He as it moves in 2D: (cid:104)V He − s ( x, y ) (cid:105) ≡ (cid:28) (cid:82) d z V He − s ( x, y, z ) ρ ( x, y, z ) (cid:82) d z ρ ( x, y, z ) (cid:29) (26)while (e) is a horizontal cut along the line y = 0highlighting that the minimum-to-saddle corrugation is V spHe − s − V minHe − s (cid:39) . V spHe − s − V minHe − s (cid:39) . z -direction and partial localization of the wavefunc-tion in the xy -plane.These QMC results for a single particle can be usedin conjunction with the band structure analysis intro-duced in § IV B to map the system to a non-interactingBose–Hubbard model. In particular, under the assump-tion that an adsorbed He atom is confined in a 2D layer,we employed (cid:104)V He − s ( x, y ) (cid:105) and extracted t from the re-sulting spectrum in Fig. 6. This is equivalent in princi-ple to using the overlap in Eq. (9) for a real wavefunc-tion | Ψ ⊥ ( x, y ) | ∝ ρ ( x, y ) where the QMC average hasbeen performed by exploiting translational invariance,i.e. moving from the Bloch to Wannier basis. The result-ing localized single particle wavefunction (labelled QMC)was previously shown in Fig. 7. We find: t QMC = 1 . . 2. Many-Body Adsorption: f > In order to investigate the effects of He–He interac-tions and thus determine the effective parameters V and V (cid:48) in the Bose–Hubbard model we need to increase thefilling fraction until He atoms occupy every site of thetriangular lattice defined by hexagon centers. However,as discussed in § III, as the density of helium atoms nearthe surface is increased, the strong repulsive interactionin Eq. (3) will cause layer completion and promote thegrowth of further layers such that the system can nolonger be considered within the 2D approximation.It is thus easier to first consider the case of f = 1 / N = 16 particles near N s = 48 adsorption sites yieldsthe 2D density profile ρ ( x, y ) shown in Fig. 10. Notethat in contrast to Fig. 9(b) for f = 1 /N s , here the localspread of the wavefunction around the hexagon centers inthe xy -plane is strongly reduced with vanishing densitybetween. The ground state is a stable solid and interac-tions are mediated through next-nearest neighbor sites ata distance of 3 a as indicated with dashed lines in anal-ogy with Fig. 2. In order to estimate the value of V (cid:48) fromthis data, we can compute the ground state energy in the2D Bose–Hubbard model in Eq. (1) for a Fock state char-acterizing the C1/3 phase, denoted by | (cid:105) , where the ki-netic energy and nearest neighbor interaction terms areidentically zero: (cid:104) | H BH | (cid:105) ≡ E BH (cid:12)(cid:12)(cid:12)(cid:12) f =1 / = 3 N V (cid:48) = N s V (cid:48) , (27)where we are neglecting the effects of even further V (cid:48)(cid:48) in-teraction terms. Thus, measuring the total contributionof the interaction potential to the ground state energy in3 − . − . . . . x / (cid:6) A − − − − y / (cid:6) A L z = 10 (cid:6) A 0 . . . . ρ ( x , y ) / (cid:6) A − FIG. 10. The two-dimensional density of particles ρ ( x, y ) ob-tained from ground state quantum Monte Carlo simulationsfor a simulation cell with L x × L y × L z = 14 . × . × . (cid:6) A corresponding to the C1/3 phase with f = 1 / N = 16 He atoms on N s = 48 adsorption sites. Finitesize effects in the spatial wavefunction are negligible beyond N s = 12. QMC, (cid:104)V He − He (cid:105) and equating this with E BH , we identify: V (cid:48) QMC = 1 N s (cid:104)V He − He (cid:105) f =1 / (28)and estimate: V (cid:48) QMC = − . V (cid:48) He − He = V He − He ( | r | = 3 a ) (cid:39) − . V ,we need to hinder the formation of multiple layers whichcan be accomplished by restricting our simulation cell inthe z -direction using (25). However, it is not clear whichvalue of L z will (1) maintain the existence of a singlewell-defined 2D monolayer as the filling is increased past f (cid:39) . f = 1 / f = 1 / , L z ∈ [4 . , . 5] and N s = 24 , , 96. Finite size effects in N s were negligiblefor the density profiles in the z -direction, and we showsimulation results for N s = 48 in Fig. 11. Here panels . . . . . . ρ ( z ) / N / (cid:6) A − N s = 48 , f = 1(a) f = 1 / L z = 10 ˚A1 . . . . . . z / (cid:6) A0 . . . . . . ρ ( z ) / N / (cid:6) A − N s = 48 , f = 1 / f = 1 / L z = 10 ˚A4 . . . . . L z / (cid:6) A FIG. 11. The density profile (per particle) of the adsorbedlayer(s) for different vertical box sizes L z enforced through thepotential in Eq. (25) for helium above a graphene sheet with N s = 48 adsorption sites such that L x × L y = 14 . 757 07 (cid:6) A × . (cid:6) A. Panels correspond to (a) filling f = 1 and (b) f = 1 / L z = 5 . (cid:6) A wasdetermined to be the optimal value (see text). The dashed lineindicates the density profile for a “bulk” cell with L z = 10 (cid:6) Aat f = 1 / Heatoms in the simulation can be determined from N = fN s . correspond to different filling fractions and colors to dif-ferent values of L z . In panel (a) at unit filling ( f = 1) weobserve that the density ρ ( z ) smoothly evolves as a func-tion of L z from one that contains a single well-definedlayer for small box sizes ( L z (cid:46) (cid:6) A), to a profile with twopeaks in the density for L z (cid:38) (cid:6) A. In order to quantifythese two regimes and determine at which value of L z we should analyze the system, we performed additionalsimulations at f = 1 / f = 1 / L z = 10 (cid:6) A beyond L z (cid:38) . (cid:6) A as indicated by thedashed line. This data can then be exploited by search-ing for the value of L z at unit filling that produces a den-sity profile most similar to that of the bulk monolayer at f = 1 / z -spread of thewavefunction. To proceed, we search for a minimum in4 . . . . . L z / (cid:6) A510152025 χ / (cid:6) A − N s = 24 N s = 48 2 4 z / (cid:6) A0 . . . ρ ( z ) / N / (cid:6) A − f = 1 / L z = 10 ˚A f = 1 L z = 5 . 05 ˚A FIG. 12. The squared deviation between solid and dashedcurves in Fig. 11(a) as a function the box size in the z -direction as quantified in Eq. (29). The minimum at L z =5 . (cid:6) A is independent of the size of the graphene sheet, wheredata for N s = 24 and 48 are shown. The inset shows a com-parison of the density profiles in the z -direction for this valueat N s = 48. The dashed line is a guide to the eye. the squared deviation of densities: χ ( L z ) ≡ (cid:88) i ρ ( z i ) N s (cid:12)(cid:12)(cid:12)(cid:12) L z f =1 − ρ ( z i ) N s / (cid:12)(cid:12)(cid:12)(cid:12) Lz =10 (cid:6) A f =1 / (29)where i runs over all spatial positions in z where densitydata has been obtained. The results of this procedure areshown in Fig. 12 and indicate a quadratic dependence on L z with the minimum occurring at L z = 5 . (cid:6) A. At thisvalue of L z , the inset shows a comparison of the twodensity profiles from Fig. 11. Finite size effects in N s were not found to alter the optimal value of L z .Recall that the goal of this procedure was to stabilizea single monolayer at filling fraction f = 1 in order todetermine the effects of nearest and next-nearest neigh-bor interactions between He atoms without substan-tially distorting the physics of the adsorbed phase. As anadditional check, we have computed the equation of stateat L z = 5 . (cid:6) A and compared it with that determined forthe unrestricted bulk cell with L z = 10 (cid:6) A for a systemwith N s = 24 adsorption sites. The results, shown inFig. 13, demonstrate that the additional confinement po-tential in Eq. (25) does not alter the ground state prop-erties of the adsorbed monolayer for f (cid:46) . 6. The in-sets show that the monolayer profile remains mostly un-changed for filling fraction f = 1 / 3. At unit filling with f = 1, while the confined box with L z = 5 . (cid:6) A still ex-hibits only a single layer, the unbounded cell can nowaccommodate an energetically favorable second layer.Thus, provided we are interested in constructing alow energy effective model at lower filling for a mono-layer, we conclude that L z = 5 . (cid:6) A is an appropri-ate choice for simulations at f = 1. In this case, theground state is an insulator as seen in the 2D density . . . . . . f = N/N s E / N − E / K N s = 24 L z = 10 (cid:6) A L z = 5 . (cid:6) A0 5 z / (cid:6) A ρ ( z ) / N / (cid:6) A − f = 1 / f = 1 / z / (cid:6) A f = 1 f = 1 FIG. 13. Equation of state (energy per particle as a functionof filling fraction) for N s = 24 adsorption sites for two boxsizes with L z = 5 . , . (cid:6) A. The curves have been shiftedby the energy for a single particle N = 1 corresponding to afilling fraction-independent value of ∼ V wall . The insets show the density of particles along the z -direction at filling fractions f = 1 / L z = 10 (cid:6) A can accommodate a second layer. in Fig. 14. Here, the wavefunction is strongly local-ized near the center of a graphene hexagon, with a cutalong y = 0 having been previously shown in Fig. 7.Following similar logic to that employed for the insulat-ing phase at f = 1 / 3, Eq. (28), we examine the Bose–Hubbard model on the triangular lattice at f = 1 where (cid:104) | H BH | (cid:105) ≡ E BH (cid:12)(cid:12) f =1 = 3 V N + V (cid:48) N and compute V QMC = 13 N (cid:104)V He − He (cid:105) f =1 − V (cid:48) QMC . (30)The results are shown in Fig. 15 as a function of L z andwe identify: V QMC = 54 . L z = 5 . (cid:6) A, where the uncertainty in the last digitarises from a combination of stochastic errors and finitesize effects. This value is larger than an estimate ob-tained from the bare interaction potential for two he-lium atoms separated by the nearest-neighbor distance: V He − He = V He − He (cid:0) | r | = √ a (cid:1) (cid:39) 31 K. While there arevery limited finite size effects in N s , the chosen value of L z does have an effect on the value of V QMC , reducingit from 61 . L z = 4 . (cid:6) A to 49 . L z = 5 . (cid:6) A.For larger values of L z , there is no longer a single well-defined monolayer and the rapid reduction in V observedin Fig. 15 can be attributed to the promotion of a secondlayer where He atoms can now move to larger valuesof z to minimize their repulsive interaction (as seen inFig. 11).5 − . − . . . . x / (cid:6) A − − − − y / (cid:6) A L z = 5 . (cid:6) A 012345 ρ ( x , y ) / (cid:6) A − FIG. 14. The two-dimensional density of particles ρ ( x, y ) ob-tained from ground state quantum Monte Carlo simulationsfor a simulation cell with L x × L y × L z = 14 . × . × . (cid:6) A corresponding to the fully filled phase with f = 1for N = 48 He atoms on N s = 48 adsorption sites. Finitesize effects in the spatial wavefunction are negligible beyond N s = 12. Nearest neighbor ( V ) and next-nearest neighbor( V (cid:48) ) couplings in the effective Bose–Hubbard description areindicated with solid and dashed lines, respectively. F. Density Functional Theory We performed DFT calculations with the PBE(Perdew–Burke–Ernzerhof) generalized gradient approx-imation [85] for the exchange–correlation functional andprojector augmented wave (PAW) [86] pseudopotentials(PPs), as implemented in the Quantum Espresso elec-tronic structure package [30, 31]. For He–graphene cal-culations, one or two He atoms are placed at a speci-fied distance from a periodic graphene sheet consistingof 6 × (cid:6) A. For He–He calculations, twoHe atoms are placed at a specified distance within a cu-bic simulation cell of 30 (cid:6) A. PAW PPs for C and He wereobtained from the standard solid-state PP library [87–89]. We applied the DFT-D4 semi-empirical dispersioncorrection [32–34] when computing single point energiesand structural optimizations to account for long-rangeelectronic–correlation effects. The energy cut-off forwavefunctions was 50 Ry (680 eV) and 360 Ry (4900 eV)for the charge density and potential. The Brillouin zoneis sampled using a Monkhorst–Pack grid with 6 × × . . . . . L z / (cid:6) A30405060 V Q M C / K V QMC = 54 . N s = 24 N s = 48 FIG. 15. The effective nearest neighbor interaction parameterof the Bose–Hubbard model computed from Eq. (30) via quan-tum Monte Carlo for simulation cells with N s = 24 , 48 as afunction of the cell size in the z -direction, L z . The indicatedvalue of V QMC = 54 . L z = 5 . (cid:6) Aas described in the text. The semi-transparent symbols for L z ≥ . (cid:6) A indicate cell sizes which allowed the nascent for-mation of a second layer, where the mapping of the micro-scopic Hamiltonian to the 2D Bose–Hubbard model breaksdown. of the potential (between centers of neighboring latticesites and passing through the saddle point), the positionof the He atom is fixed in the plane of the sheet and theoptimal distance from the sheet is then found at eachpoint to compute the energy along the minimum energysurface (see Fig. 18). We followed the same approach tofind the maximum value of he potential (centered at theposition of a C atom) with the results shown in Table II.For He–He on graphene calculations, two He atomsare placed at the centers of various lattice sites and atan optimal distance from the sheet, obtained beforehandfor a single He atom ( z opt (cid:39) . (cid:6) A, see Fig. 4). The re-sulting interaction (relative to non-interacting adsorbedatoms) provides an estimate for the nearest and next-nearest neighbor values: V DFT = 21 . V (cid:48) DFT = − . 36 K . (31) G. Møller–Plesset Perturbation Theory Because the He–He and He–graphene interactions aredominated by dispersion terms which require accuratetreatment of the correlation energy [36, 37], second-orderMøller–Plesset (MP2) [35] perturbation theory calcula-tions, which in most cases capture ca. 95% of the corre-lation energy [36], were performed using Gaussian 09 [90]utilizing Pople-type [91] bases sets up to 6-31++G(d,3p),which include diffusion of all orbitals, and polarizationfunctions d and p for all atoms.Such a high-order basis set was needed to obtain theHe–He interactions in vacuum to reasonable accuracy6(Fig. 4(a)). To model the interaction of He atom(s)with graphene (and possible modifications of the He–Hepotential on graphene), a sequence of increasing aromaticmolecules was considered (benzene, coronene, hexaben-zocoronene, circumcoronene — the latter with 54 C and18 H atoms). The energy of the system was computedfor different values of z between the He atom(s) and theC plane, and the asymptotic energy for z → ∞ was re-moved as a baseline. To reproduce graphene, the aro-matic molecules were constructed with C–C distancesconstrained to a = 1 . (cid:6) A, and only the coordinates ofterminating H atoms were optimized. Figure 4(b) showsthe potential energy vs. height for a single He atom abovethe center of a circumcoronene molecule. We observedthat the calculations converge after hexabenzocoroneneand there was a relatively small “radial dependence” of V He − s ( r , z ) for other hexagon centers, making this a rea-sonable model for He on graphene.We also performed scans of the potentials over differ-ent positions over the circumcoronene. Figure 8 depictsthe dependence of V He − s ( x, y, z ), i.e., the lateral de-pendence of the He–graphene minimum energy surface(values in Table II) which allows for the calculation of t MP2 reported in § IV B.Additionally, we performed calculations for the energyfor two He atoms adsorbed onto various hexagon centers.After removing the baseline 2 × V He − s terms, we find aremnant V He s − He s ( r ) which remains strongly repulsivefor nearest neighbors ( r = 2 . 46 ˚A) and attractive fornext-nearest neighbors and beyond ( r ≥ . 26 ˚A): V MP2 = 51 . , V (cid:48) MP2 = − . 97 K . (32) V. DISCUSSION Our main result is the construction of a reliable andconsistent description of the effective two-dimensionaladsorption problem of helium-4 on graphene in the lan-guage of the hard-core Bose–Hubbard Model, Eq. (1).The relevant hopping and interaction parameters com-puted via different techniques are summarized in TableI. The differences can be intuitively understood by ex-amining Fig. 4. For example, density functional theorypredicts a deeper (more attractive) He–graphene and He–He potential (compared to the empirical, Lennard–Jonesparametrized potential). The resulting corrugation of theadsorption potential is enhanced (see Table II) leading toa suppressed hopping t , while the increased two-body at-traction between He atoms results in more spatially local-ized wavefunctions (in the xy -plane) that yield a stronglyreduced V , as the effects of the hard-core overlap are sup-pressed.On the other hand, the Møller–Plesset perturbativemethod gives the strongest He–graphene interaction,leading to a smaller hopping t , while the He–He inter-action is close to the empirical one, and it gives similarvalues of V . Perhaps most importantly, the quantum Monte Carlo and Hartree–Fock methods, both based onthe empirical potentials, lead to similar results for allparameters.Overall, a remarkably consistent picture emerges. Allof the above methods take into account the strong many-body He–He correlations on the scale of several (cid:6) A, com-parable to the localization length of the one-particlewavefunctions in the lattice field of graphene. Thesimple one-particle Wannier description fails completelyin this case, and thus the many-body techniques de-scribed in this work are essential to capture the self-consistent, interaction-driven adjustment of the one-particle orbitals, in turn leading to significant changesin the effective He–He interactions on the lattice scale.To the best of our knowledge, this represents a uniquecase of a Bose–Hubbard model construction outside theusual examples that involve cold atom systems in opti-cal lattice potentials. Moreover, the Bose–Hubbard classof models that appear in cold atoms are much simplerto define and parametrize due to the diluteness of theatomic gases involved, which implies that the details ofthe atom–atom interactions on short scales are not im-portant. By contrast, for the case of helium on graphene,the details of the small distance He–He potential on thescale of the graphene lattice are extremely important,and consequently the effective Bose–Hubbard parametersare very sensitive to the microscopic form of the poten-tial employed. Due to the strong short-range repulsion,our model describes hard-core bosons ( U ≈ ∞ ), with fi-nite nearest neighbor repulsion V > V (cid:48) < ACKNOWLEDGMENTS We dedicate this paper to our late colleague, Dr. Dar-ren Hitt, former director of the VT Space Grant Con-sortium. Darren’s leadership, encouragement, and vi-sion to expand the scope of space grant activities inVermont was essential to forming our interdisciplinarycollaboration. This work was supported, in part, underNASA grant number 80NSSC19M0143. Computationalresources were provided by the NASA High-End Com-puting (HEC) Program through the NASA Advanced Su-percomputing (NAS) Division at Ames Research Center. Appendix A: Dimensional Reduction of theAdsorption Layer to 2D The He–graphene potential in Eq. (4) can be writtenas V He − s ( r ) = V ( z ) + V ⊥ ( r , z ) , (A1)where V ( z ) and V ⊥ ( r , z ) are the g = 0 and g (cid:54) = 0 terms,respectively. A justification for such a splitting was pre-sented in § III. We seek the single-particle minimizer ofthe Hamiltonian in Eq. (3), i.e. the solution of: − (cid:126) m ∇ Ψ( r ) + V He − s ( r )Ψ( r ) = E Ψ( r ) , (A2)as a (truncated) expansion over eigenfunctions { φ n ( z ) } of the 1D potential V ( z ):Ψ( r ) ≡ Ψ( r , z ) = (cid:88) n χ n ( r ) φ n ( z ) ; (A3a)where φ n satisfy − (cid:126) m φ (cid:48)(cid:48) n + V φ n = (cid:15) n φ n , (cid:104) φ n | φ n (cid:105) = 1 , (A3b)and (cid:104) . . . (cid:105) stands for integration over z . SubstitutingEq. (A3a) into Eq. (A2), multiplying by (cid:104) φ n | and inte-grating over z , one obtains a system of coupled equationsfor χ n ( r ). Such a system can, in principle, be solved bythe same numerical method as described in § IV B. To fo-cus on the conceptual consequences of the z -spread of thewavefunction rather than on finer details, we proceed bytruncating the expansion in Eq. (A3a) at lowest order:Ψ( r , z ) ≈ χ ( r ) φ ( z ) . (A4a) The quantity χ ( r ) plays the role of an effective “reduced”2D wavefunction and satisfies the Schr¨odinger equation: − (cid:126) m ∇ r χ + (cid:101) V He − s χ = ˜ Eχ , (A4b) (cid:101) V He − s ( r ) ≡ (cid:104) φ |V ( z ) + V ⊥ ( r , z ) | φ (cid:105) , (A4c)where ∇ r is the Laplacian in r . We have absorbed aconstant (cid:104) φ |V ( z ) | φ (cid:105) into both sides of (A4b) to ob-tain correspondence with the quantum Monte Carlo z -averaging results described in § IV E. Equation (A4b)represents the 2D reduction of the 3D one-particle model,which is solved in § IV B.We note that the particle density corresponding to theapproximation (A4a) is: ρ ( r , z ) = | ψ ( r ) | ρ ( z ) , ρ ( z ) ≡ | φ ( z ) | . (A5)Substituting this into (26), one finds that the expressionthere coincides with (cid:101) V He − s in (A4c). The approximateparticle density in the z -direction, ρ ( z ), can be foundby solving Eq. (A3b) with n = 0 by, e.g. , the shootingmethod. The result is shown in Fig. 5.The spread of the single-particle density ρ ( z ) in the z -direction also affects the computation of the effectivenearest neighbor V in the Bose-Hubbard model. Thisspread leads to nearest-neighbor He atoms having a dis-tribution of z -values relative to one another [61] in adja-cent graphene hexagons: γ ( δ ) = (cid:90) ρ ( z ) ρ ( z + δ ) dz . (A6)Here we have assumed, in agreement with the finding byQMC simulations, that there is no correlation between z -values of the centers of nearest-neighbor He atoms.This quantity can then be used to estimate the effect ofthe relative z -spread δ of nearest-neighbor He atoms ontheir interaction via: (cid:101) V He − He ( r ) = (cid:90) V He − He (cid:0)(cid:112) | r | + δ (cid:1) γ ( δ ) dδ . (A7)It leads to some “softening” of the He–He interactionpotential. In § IV D we show how much this softeningaffects V . Arguably, to more properly account for thespread in the z -direction, instead of using the ρ ( z ) de-fined in (A5) one would need to use ρ ( z ) obtained byMonte Carlo simulations ( § IV E). The results of usingboth approaches are compared in § IV D. Appendix B: Simulation Details and Scaling In this appendix, we provide details on the quantumMonte Carlo method used in § IV E. Access to the em-ployed software can be obtained via Ref. [84].By choosing a suitably small imaginary time step τ (polynomial scaling) and long enough imaginary time8projection length β (exponential dependence) it is possi-ble to ensure that any systematic errors inherent in thechoice of an approximate propagator ρ τ [92, 93] are madesmaller than any statistical errors in our ground statequantum Monte Carlo scheme. At the additional com-putational expense of requiring a potentially larger valueof β to obtain convergence, we have chosen to employthe constant trial wavefunction Ψ T ( R ) = 1 to preventthe breaking of translational symmetry of the adsorbedphase in order to explore the effects of particle tunneling.Figure 16 shows the convergence of the energy as afunction of τ and β . We choose τ (cid:39) . 003 13 K cor- h K / N i / K (a) N s = 4 , f = 1 QMC K + c e − c β − − h V H e − s / N i / K (b) 0 . 25 0 . 50 0 . 75 1 . 00 1 . 25 1 . β / K − . . . . h V H e − H e / N i / K (c) 0 . 00 0 . 01 0 . 02 0 . τ / K − . − . − . V He − s , + c τ + c τ . 00 0 . 01 0 . 02 0 . τ / K100150 FIG. 16. Convergence of the kinetic (a), adsorption (b), andinteraction potential (c) energy per particle with the imagi-nary time projection length β and imaginary time step τ (in-sets) for a system with N s = 4 adsorption sites at unit filling.All quantities measured via quantum Monte Carlo (symbols+ errorbars) are observed to converge to their ground statevalue with the expected exponential dependence on β or poly-nomial dependence on τ as quantified via the indicated fits(lines). β scaling was performed with τ = 1 / 319 K (verticaldashed line in insets) while τ scaling had β (cid:39) . − . Asubscript 0 indicates the ground state value. responding to 319 discrete imaginary time steps and β (cid:39) . − to 1 . − for all simulations presented inthis work.While the spatial extent of the simulation cell in the x and y directions had a minimal effect on most observables(see discussion and results in § IV E) there was some ob- served dependence of V (cid:48) on N s as shown in Fig. 16. This /N s − . − . − . − . V Q M C / K L z = 10 (cid:6) AQMCfit: − . . /N s FIG. 17. The finite size scaling of the next nearest neigh-bor interaction V (cid:48) extracted from quantum Monte Carlo as afunction of the number of graphene hexagon adsorption sites.The solid line shows a linear fit to the data for N s ≥ 48 whichallows extrapolation to the thermodynamic limit. is likely due to the finite size configuration at f = 1 / N s = ∞ to obtain the reported result for V (cid:48) . Appendix C: Quasiclassical approximation The one-dimensional (1D) WKB approach is a semi-classical approach which can portray tunnel splitting andis well defined in 1D double-well or periodic potentials.In order to obtain an intuitive and simple estimate ofthe hopping t we will apply the 1D approach along thepath passing through the saddle point in our 2D poten-tial, Fig. 8. This path connects two adjacent minima(Fig. 18) and leads to the largest hopping amplitude.In this quasi 1D limit the energy dispersion along anyof the three triangular lattice directions has the form: ε ( k ) − ε = − t cos ( ka ). The quasiclassical expressionfor the hopping is known [78, 94] to be: t = (cid:126) ω π exp (cid:32) − (cid:126) (cid:90) x c x c (cid:112) m ( V He − s ( x ) − (cid:126) ω / dx (cid:33) , (C1)where the classical turning points x c , satisfy V He − s (cid:0) x c , (cid:1) = (cid:126) ω / ω is the frequency of smallamplitude oscillations in the wells which are fitted toparabolic (harmonic oscillator) form. This expressionis valid as long as the potential barrier is high enoughand the exponential tunneling factor is small, whichis only approximately satisfied in our system. Overallthe hopping is a product of the tunneling factor andthe attempt frequency ω , leading to a finite number9 − − x / (cid:6) A0204060 V H e − s / K MP2DFTQMC FIG. 18. A one dimensional cut through V He − s for a He atom above graphene corresponding to the spatial paththrough the saddle point of the full 2D potential shown inFig. 8. The energy barriers that inhibit tunneling between dif-ferent graphene adsorption sites originate from three differentmethods: Møller–Plesset perturbation theory (MP2, green),density functional theory (DFT, blue), and quantum MonteCarlo (QMC, red) as described in the text, where the curveshave been shifted with respect to the bottom of the poten-tial well. The classical turning points x c , for each potentialbarrier appear at the intersection of the dashed ( E = (cid:126) ω / (cid:126) ω / (K) t / (K) t/V QMC 15.1 2.22 0.041DFT 19.6 1.25 0.058MP2 27.3 0.43 0.008TABLE III. Ground state energy in the potential wells in theharmonic approximation ( (cid:126) ω / 2) and quasi-classical hoppingparameters t , Eq. (C1), derived for He–graphene interactionsobtained by various methods, Fig. 18. expected to provide a good numerical estimate. Thehopping parameter derived from He–graphene interac-tions obtained by various methods (Fig. 18) are shownin Table III. It is clear that the results from this simpleapproximation provide quite reasonable estimates asthey are comparable to the numbers and tendenciesfrom the full 2D calculations whose results are displayedin Table I. [1] M. Bretz and J. Dash, Quasiclassical and Quantum De-generate Helium Monolayers, Phys. Rev. Lett. , 963(1971).[2] D. S. Greywall, Heat capacity and the commensurate-incommensurate transition of He4 adsorbed on graphite,Phys. Rev. B , 309 (1993).[3] J. G. Dash, M. Schick, and O. E. Vilches, Phases of he-lium monolayers: search and discovery, Surf. Sci. ,405 (1994).[4] F. M. Gasparini, M. O. Kimball, K. P. Mooney, andM. Diaz-Avila, Finite-size scaling ofHe4at the superfluidtransition, Rev. Mod. Phys. , 1009 (2008).[5] L. Reatto, D. E. Galli, M. Nava, and M. W. Cole,Novel behavior of monolayer quantum gases on graphene,graphane and fluorographene, J. Phys.: Condens. Mat. , 443001 (2013).[6] T. Makiuchi, M. Tagai, Y. Nago, D. Takahashi, andK. Shirahama, Elastic anomaly of helium films at a quan-tum phase transition, Phys. Rev. B , 235104 (2018).[7] J. Saunders, Realizing quantum materials with helium:Helium films at ultralow temperatures, from strongly cor-related atomically layered films to topological superflu-idity, in Topological Phase Transitions and New Develop-ments (World Scientific, 2018) p. 165.[8] J. Saunders, B. Cowan, and J. Ny´eki, Atomically layeredhelium films at ultralow temperatures: Model systems forrealizing quantum materials, J. Low Temp. Phys. ,615 (2020).[9] A. Thomy and X. Duval, Adsorption de mol´ecules sim-ples sur graphite, J. Chim. Phys. , 1966 (1969).[10] P. A. Crowell and J. D. Reppy, Superfluidity and filmstructure in He4 adsorbed on graphite, Phys. Rev. B , 2701 (1996).[11] J. Ny´eki, R. Ray, G. Sheshin, V. Maidanov, V. Mikheev,B. Cowan, and J. Saunders, Structure and superfluidityof 4He films on plated graphite, J. Low Temp. Phys. ,379 (1997).[12] F. Abraham and J. Broughton, Phases of helium ad-sorbed on graphite: A Feynman path-integral MonteCarlo study, Phys. Rev. Lett. , 64 (1987).[13] M. Pierce and E. Manousakis, Path-integral Monte Carlosimulation of the second layer of 4He adsorbed ongraphite, Phys. Rev. B , 3802 (1999).[14] M. C. Gordillo and J. Boronat, Superfluid and SupersolidPhases of He4 on the Second Layer of Graphite, Phys.Rev. Lett. , 205301 (2020).[15] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, The electronic properties ofgraphene, Rev. Mod. Phys. , 109 (2009).[16] N. S. Nichols, A. D. Maestro, C. Wexler, and V. N. Kotov,Adsorption by design: Tuning atom-graphene van derWaals interactions via mechanical strain, Phys. Rev. B , 205412 (2016).[17] M. Gordillo and J. Boronat, 4He on a Single GrapheneSheet, Phys. Rev. Lett. , 085303 (2009).[18] J. Happacher, P. Corboz, M. Boninsegni, and L. Pollet,Phase diagram of 4He on graphene, Phys. Rev. B ,094514 (2013).[19] K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H.Castro Neto, 2D materials and van der Waals het-erostructures, Science , 461 (2016).[20] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea,and A. H. Castro Neto, Electron-electron interactions ingraphene: Current status and perspectives, Rev. Mod. Phys. , 1067 (2012).[21] L. Bruch, M. W. Cole, and H.-Y. Kim, Transitions ofgases physisorbed on graphene, J. Phy.: Condens. Mat. , 304001 (2010).[22] M. C. Gordillo, C. Cazorla, and J. Boronat, Supersolidityin quantum films adsorbed on graphene and graphite,Phys. Rev. B , 121406 (2011).[23] L. V. Marki´c, P. Stipanovi´c, I. Beˇsli´c, and R. E. Zil-lich, 4He clusters adsorbed on graphene, Phys. Rev. B (2013).[24] Y. Kwon and D. M. Ceperley, 4He adsorption on a singlegraphene sheet: Path-integral Monte Carlo study, Phys.Rev. B , 224501 (2012).[25] M. C. Gordillo and J. Boronat, Zero-temperature phasediagram of the second layer of 4He adsorbed on graphene,Phys. Rev. B , 195457 (2012).[26] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, andP. Zoller, Cold Bosonic Atoms in Optical Lattices, Phys.Rev. Lett. , 3108 (1998).[27] D. Jaksch and P. Zoller, The cold atom Hubbard toolbox,Ann. Phys. , 52 (2005).[28] I. Bloch, J. Dalibard, and W. Zwerger, Many-bodyphysics with ultracold gases, Rev. Mod. Phys. , 885(2008).[29] R. Walters, G. Cotugno, T. H. Johnson, S. R. Clark, andD. Jaksch, Ab initio derivation of Hubbard models forcold atoms in optical lattices, Phys. Rev. A , 043613(2013).[30] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz-covitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of ma-terials, J Phys Condens Matter , 395502 (2009).[31] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau,M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavaz-zoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carn-imeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A.DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo,R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia,M. Kawamura, H.-Y. Ko, A. Kokalj, E. K¨u¸c¨ukbenli,M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L.Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto,S. Ponc´e, D. Rocca, R. Sabatini, B. Santra, M. Schlipf,A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thon-hauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Ad-vanced capabilities for materials modelling with Quan-tum ESPRESSO., J Phys.: Condens. Matt. , 465901(2017).[32] E. Caldeweyher, C. Bannwarth, and S. Grimme, Exten-sion of the D3 dispersion coefficient model, J. Chem.Phys. , 034112 (2017).[33] E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer,S. Spicher, C. Bannwarth, and S. Grimme, A generallyapplicable atomic-charge dependent London dispersioncorrection, J. Chem. Phys. , 154122 (2019).[34] E. Caldeweyher, J.-M. Mewes, S. Ehlert, and S. Grimme,Extension and evaluation of the D4 London-dispersionmodel for periodic systems, Phys. Chem. Chem. Phys. , 8499 (2020).[35] C. Møller and M. S. Plesset, Note on an ApproximationTreatment for Many-Electron Systems, Phys. Rev. ,618 (1934).[36] C. Cramer, Essentials of Computational Chemistry: The-ories and Models (Wiley, 2013).[37] R. J. Bartlett, Many-Body Perturbation Theory andCoupled Cluster Theory for Electron Correlation inMolecules, Ann. Rev. Phys. Chem. , 359 (1981).[38] C. Becker, P. Soltan-Panahi, J. Kronj¨ager, S. D¨orscher,K. Bongs, and K. Sengstock, Ultracold quantum gasesin triangular optical lattices, New J. Phys. , 065025(2010).[39] J. Iba˜nez-Azpiroz, A. Eiguren, A. Bergara, G. Pettini,and M. Modugno, Tight-binding models for ultracoldatoms in honeycomb optical lattices, Phys. Rev. A (2013).[40] R. A. Aziz, V. P. S. Nain, J. S. Carley, W. L. Taylor, andG. T. McConville, An accurate intermolecular potentialfor helium, J. Chem. Phys. , 4330 (1979).[41] R. A. Aziz, F. R. McCourt, and C. C. Wong, A newdetermination of the ground state interatomic potentialfor He2, Mol. Phys. , 1487 (1987).[42] R. A. Aziz, A. R. Janzen, and M. R. Moldover, Ab InitioCalculations for Helium: A Standard for Transport Prop-erty Measurements, Phys. Rev. Lett. , 1586 (1995).[43] M. Przybytek, W. Cencek, J. Komasa, G. (cid:32)Lach,B. Jeziorski, and K. Szalewicz, Relativistic and Quan-tum Electrodynamics Effects in the Helium Pair Poten-tial, Phys. Rev. Lett. , 183003 (2010).[44] W. Cencek, M. Przybytek, J. Komasa, J. B. Mehl,B. Jeziorski, and K. Szalewicz, Effects of adiabatic, rel-ativistic, and quantum electrodynamics interactions onthe pair potential and thermophysical properties of he-lium, J. Chem. Phys. , 224303 (2012).[45] K. T. Tang, J. P. Toennies, and C. L. Yiu, AccurateAnalytical He-He van der Waals Potential Based on Per-turbation Theory, Phys. Rev. Lett. , 1546 (1995).[46] R. E. Grisenti, W. Sch¨ollkopf, J. P. Toennies, G. C.Hegerfeldt, T. K¨ohler, and M. Stoll, Determination of theBond Length and Binding Energy of the Helium Dimerby Diffraction from a Transmission Grating, Phys. Rev.Lett. , 2284 (2000).[47] W. L. McMillan, Ground State of Liquid He4, Phys. Rev. , 442 (1965).[48] P. A. Whitlock, D. M. Ceperley, G. V. Chester, and M. H.Kalos, Properties of liquid and solid He4, Phys. Rev. B , 5598 (1979).[49] Y. Lutsyshyn, Coordinated wave function for the groundstate of liquid He4, Phys. Rev. B , 214507 (2015).[50] E. Kaxiras and J. D. Joannopoulos, Quantum Theory ofMaterials (Cambridge University Press, 2019).[51] C. J. Pethick and H. Smith, Bose-Einstein condensationin dilute gases , 2nd ed. (Cambridge University Press,Cambridge, 2008).[52] N. Gheeraert, S. Chester, M. May, S. Eggert, and A. Pel-ster, Mean-Field Theory for Extended Bose-HubbardModel with Hard-Core Bosons, in Selforganization inComplex Systems: The Past, Present, and Future of Syn-ergetics , edited by G. Wunner and A. Pelster (SpringerInternational Publishing, 2016).[53] G. Murthy, D. Arovas, and A. Auerbach, Superfluids andsupersolids on frustrated two-dimensional lattices, Phys.Rev. B , 3104 (1997). [54] S. Wessel and M. Troyer, Supersolid Hard-Core Bosonson the Triangular Lattice, Phys. Rev. Lett. , 127205(2005).[55] J.-Y. Gan, Y.-C. Wen, and Y. Yu, Supersolidity andphase diagram of soft-core bosons on a triangular lattice,Phys. Rev. B , 094501 (2007).[56] X.-F. Zhang, R. Dillenschneider, Y. Yu, and S. Eggert,Supersolid phase transitions for hard-core bosons on atriangular lattice, Phys. Rev. B , 174515 (2011).[57] (2021), All code, scripts and data used inthis work are included in a GitHub reposi-tory: https://github.com/DelMaestroGroup/papers-code-BoseHubbardModelHeAdsorptionGraphene ,permanent link: https://doi.org/10.5281/zenodo.4553595 .[58] W. A. Steele, The physical interaction of gases with crys-talline solids, Surf. Sci. , 317 (1973).[59] W. E. Carlos and M. W. Cole, Anisotropic He-C pairInteraction for a He Atom Near a Graphite Surface, Phys.Rev. Lett. , 697 (1979).[60] W. E. Carlos and M. W. Cole, Interaction between a heatom and a graphite surface, Surf. Sci. , 339 (1980).[61] G. Vidali and M. W. Cole, Effective interaction betweenHe atoms on a graphite surface, Phys. Rev. B , 4661(1980).[62] F. Pirani, D. Cappelletti, and G. Liuti, Range, strengthand anisotropy of intermolecular forces in atom-moleculesystems: an atom-bond pairwise additivity approach,Chem. Phys. Lett. , 286 (2001).[63] F. Pirani, M. Alberti, A. Castro, M. M. Teixidor, andD. Cappelletti, Atom-bond pairwise additive representa-tion for intermolecular potential energy surfaces, Chem.Phys. Lett. , 37 (2004).[64] L. Bruch, M. Cole, and E. Zaremba, Phys. Adsorption:Forces and Phenomena , Dover Books on Physics (DoverPublications, 2007).[65] T. L. Badman and J. M. McMahon, On the Phase Dia-grams of 4He Adsorbed on Graphene and Graphite fromQuantum Simulation Methods, Crystals , 202 (2018).[66] M. C. Gordillo, Diffusion Monte Carlo calculation of thephase diagram of 4He on corrugated graphene, Phys.Rev. B , 155401 (2014).[67] R. Burganova, Y. Lysogorskiy, O. Nedopekin, andD. Tayurskii, Adsorption of Helium Atoms on Two-Dimensional Substrates, J. Low Temp. Phys. , 1(2016).[68] C. E. Campbell, F. J. Milford, A. D. Novaco, andM. Schick, Helium-Monolayer Completion on Graphite,Phys. Rev. A , 1648 (1972).[69] P. A. Whitlock, G. V. Chester, and B. Krishnamachari,Monte Carlo simulation of a helium film on graphite,Phys. Rev. B , 8704 (1998).[70] G. Zimmerli, G. Mistura, and M. Chan, Third-soundstudy of a layered superfluid film, Phys. Rev. Lett. ,60 (1992).[71] G. Zimmerli and M. H. W. Chan, Complete wetting ofhelium on graphite, Phys. Rev. B , 8760 (1988).[72] Y. Shibayama, H. Fukuyama, and K. Shirahama, Tor-sional oscillator studies for possible supersolidity in two-dimensional 4He solid, J. Phys.: Conf. Ser. , 032096(2009).[73] S. Nakamura, K. Matsui, T. Matsui, and H. Fukuyama,Possible quantum liquid crystal phases of helium mono-layers, Phys. Rev. B , 180501 (2016). [74] J. Ny´eki, A. Phillis, A. Ho, D. Lee, P. Coleman, J. Parpia,B. Cowan, and J. Saunders, Intertwined superfluid anddensity wave order in two-dimensional 4He, Nature Phys. , 455 (2017).[75] M. E. Pierce and E. Manousakis, Role of substrate cor-rugation in helium monolayer solidification, Phys. Rev.B , 5228 (2000).[76] J. Ahn, H. Lee, and Y. Kwon, Prediction of stable C7/12and metastable C4/7 commensurate solid phases for 4Heon graphite, Phys. Rev. B , 064511 (2016).[77] L. V. Marki´c, P. Stipanovi´c, I. Beˇsli´c, and R. E. Zil-lich, Solidification of 4He clusters adsorbed on graphene,Phys. Rev. B , 045428 (2016).[78] E. Lifshitz and L. P. Pitaevskii, Statistical Physics , Part2 (Pergamon Press, 1980).[79] G. D. Mahan, Many-Particle Physics (Plenum Press,1990).[80] T. I. Lakoba, Convergence conditions for iterative meth-ods seeking multi-component solitary waves with pre-scribed quadratic conserved quantities, Math. Comput.Simul. , 1572 (2011).[81] A. Sarsa, K. E. Schmidt, and W. R. Magro, A path in-tegral ground state method, J. Chem. Phys. , 1366(2000).[82] J. E. Cuervo, P.-N. Roy, and M. Boninsegni, Path inte-gral ground state with a fourth-order propagator: Appli-cation to condensed helium, J. Chem. Phys. , 114504(2005).[83] Y. Yan and D. Blume, Path integral Monte Carlo groundstate approach: formalism, implementation, and applica-tions, J. Phys. B: At., Mol. Opt. Phys. , 223001 (2017).[84] A. Del Maestro, Available online, Del Maestro GroupCode Repository (2021), https://code.delmaestro.org .[85] J. Perdew, K. Burke, and M. Ernzerhof, Generalized Gra-dient Approximation Made Simple, Phys. Rev. Lett. ,3865 (1996).[86] P. Bl¨ochl, Projector augmented-wave method, Phys. Rev.B , 17953 (1994).[87] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet,and N. Marzari, Precision and efficiency in solid-statepseudopotential calculations, npj Comput. Mater. , 72(2018).[88] A. Dal Corso, Pseudopotentials periodic table: From Hto Pu, Comput. Mater. Sci. , 337 (2014).[89] M. Schlipf and F. Gygi, Optimization algorithm for thegeneration of ONCV pseudopotentials, Comput. Phys.Comm. , 36 (2015).[90] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria,M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone,B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Cari-cato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino,G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toy-ota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima,Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Mont-gomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark,J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov,R. Kobayashi, J. Normand, K. Raghavachari, A. Ren-dell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi,N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B.Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts,R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi,C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma,V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannen- berg, S. Dapprich, A. D. Daniels, ¨O. Farkas, J. B. Fores-man, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaus-sian 09 Revision E.01 (2009), Gaussian Inc. WallingfordCT 2009.[91] W. J. Hehre, R. F. Stewart, and J. A. Pople, Self-Consistent Molecular-Orbital Methods. I. Use of Gaus-sian Expansions of Slater-Type Atomic Orbitals, J.Chem. Phys. , 2657 (1969). [92] S. A. Chin, Symplectic integrators from composite oper-ator factorizations, Phys. Lett. A , 344 (1997).[93] S. Jang, S. Jang, and G. A. Voth, Applications of higherorder composite factorization schemes in imaginary timepath integral simulations, J. Chem. Phys. , 7832(2001).[94] B. R. Holstein, Semiclassical treatment of the periodicpotential, Am. J. Phys.56