Two-dimensional Bose-Hubbard model for helium on graphene

An exciting development in the field of correlated systems is the possibility of realizing two-dimensional (2D) phases of quantum matter. For a system of bosons, an example of strong correlations 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\approx \infty )$, repulsive nearest-neighbor $(V>0)$ and small attractive $(V^{\prime }<0)$ next-nearest-neighbor interactions. The mapping onto the Bose-Hubbard model is accomplished by a variety of many-body techniques which take into account the strong He-He correlations on the scale of the graphene lattice spacing. Unlike the case of dilute ultracold atoms where interactions are effectively pointlike, the detailed microscopic form of the short-range electrostatic and long-range dispersion interactions in the helium-graphene system is crucial for the emergent Bose-Hubbard description. The result places the ground state of the first layer of ${}^4\mathrm{He}$ adsorbed on graphene deep in the commensurate solid phase with $1/3$ of the sites on the dual triangular lattice occupied. Because the parameters of the effective Bose-Hubbard model are very sensitive to the exact lattice structure, this opens up an avenue to tune quantum phase transitions in this solid-state system.

Figure 1

A schematic of the $\sqrt{3}\times \sqrt{3}$ commensurate solid phase for ${}^4\mathrm{He}$ atoms (blue) adsorbed on graphene (gold) showing $1/3$ of the sites filled where the He atoms are localized on the sites of a triangular lattice (see shadow). The size of the indicated region of graphene is $L_x\times L_y\simeq 22\AA \times 17\AA$.

Figure 2

The triangular lattice defined by hexagon centers of the graphene lattice (shown in grey) with the lattice and basis vectors in Eq. (2) shown in the upper right. The nearest-neighbor hopping $t$ and interaction $V$ are experienced between sites separated by $\sqrt{3}{a}_0$ (blue) where $a_0\simeq 1.42\text{Å}$ is the length of a hexagon side. The dashed red lines indicate the next-nearest-neighbor interaction $V^{\prime }$ between sites separated by $3a_0$ which form a triangular superlattice at $1/3$ filling (open circles) which corresponds to the adsorbed C1/3 solid depicted in Fig. 1. Lattice vectors ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$ are shared between the graphene and triangular lattices.

Figure 3

The mean-field phase diagram for hard-core bosons on the triangular lattice with nearest-neighbor interactions $V$ and hopping $t$ with density controlled by the chemical potential $\mu$. Identified phases include commensurate solids (at fillings $f=1/3\text{and}2/3$), superfluid, and supersolid (a superfluid that breaks triangular lattice symmetries). Solid lines indicate discontinuous (first-order) transitions, while continuous (second-order) transitions occur across dashed lines. The data points in the lower left-hand corner represent the major results of this paper, indicating that the ground state of a single layer of ${}^4\mathrm{He}$ on graphene resides deep in the commensurate solid phase at $1/3$ filling. For these data points, the chemical potential has been chosen such that $\mu /V$ has the same value as the tip of the first lobe.

Figure 4

The interaction (a) and adsorption (b) potentials in Eq. (3). The three curves in both panels indicate three different approaches to the potentials. Empirical indicates the formulas discussed in the text where ${\mathcal{V}}_{\mathrm{He}-\mathrm{He}}$ is taken from Ref. [46] while is determined from Eq. (4) at $\mathcal{r}=(0,0)$. The potentials labeled DFT and MP2 are extracted from the minimal energy surfaces determined by those ab initio methods as described in Secs. 4f and 4g.

Figure 5

The particle density, $\rho \left(z\right)\propto |{\varphi }_0{\left(z\right)|}^2$ in Eq. (A5) obtained via the shooting method for computed for the three different approaches to the potentials described in Secs. 2a and 2b. The table indicates the value of $z$ at which the maximum density occurs, the root mean squared value of $z$ $\mathrm{\Delta }{z}_{\mathrm{rms}}=\sqrt{\langle z^2\rangle -{\langle z\rangle }^2}$, and the difference between the zero point energy and the minimum of for each potential. All three methods yield density profiles for the adsorbed layer that are effectively two-dimensional with sub-angstrom widths.

Figure 6

The band structure obtained using as introduced Eq. (A4c) along a high-symmetry path in the first Brillouin zone as shown in the inset. The first band (with $n=1$) is well separated from the higher excited bands and thus the low energy properties of the system are determined by the lowest band. The dashed line shows the tight-binding dispersion from Eq. (15), in excellent agreement with the continuum model supporting the use of an effective 2D lattice model. The shaded region corresponds to the bandwidth equal to $9t$ by symmetry.

Figure 7

Spatial dependence of the density $\rho (x,y)={\left|\psi (x,y)\right|}^2/N$ of an adsorbed ${}^4\mathrm{He}$ atom on graphene for a $y=0$ cut in the $xy$-plane corresponding to the lattice path shown in the top inset. (a) A comparison of the localized Wannier function defined in Eq. (13) to that computed via quantum Monte Carlo (QMC) for a single particle $N=1$, by slightly biasing a single site at the level of the trial wave function as discussed in Sec. 4e. The wave function strongly penetrates into neighboring lattice sites, leading to the breakdown of Wannier theory for the computation of interaction parameters. (b) The density at unit filling as computed via QMC, and within the Hartree-Fock and Hartree approximations showing the tendency towards exponential localization on a single site. In both panels, the 2D normalization is computed over the full graphene sheet.

Figure 8

The effective 2D potential used to calculate the hopping $t$ can be reconstructed from MP2 and DFT calculations by determining three values corresponding to the minimum, maximum, and saddle-point values as indicated. The resulting scale along the black line can be seen in Fig. 18.

Figure 9

Quantum Monte Carlo results for a single ${}^4\mathrm{He}$ atom ($N=1$) above a graphene membrane with adsorption sites. (a) The average density in the $z$ direction showing a well defined particle position a distance $2.71\text{Å}$ above the sheet with a width of $0.64\text{Å}$. (b) The average density in the $xy$ plane showing the ability of a single particle to hop between the sites of the triangular lattice. White dots show the location of carbon atoms (not to scale). (c) The average potential energy experienced by the ${}^4\mathrm{He}$ atom due to the graphene sheet in the $xy$ plane. (d) A horizontal cut of the particle density $\rho (x,y)$ along the line $y=0$. (e) A horizontal cut of the adsorption potential along the line $y=0$.

Figure 10

The two-dimensional density of particles $\rho (x,y)$ obtained from ground state quantum Monte Carlo simulations for a simulation cell with $L_x\times L_y\times L_z=14.75707\times 17.04\times 10.0\text{Å}$ corresponding to the C1/3 phase with $f=1/3$ for $N=16$ ${}^4\mathrm{He}$ atoms on adsorption sites. Finite size effects in the spatial wave function are negligible beyond

Figure 11

The two-dimensional density of particles $\rho (x,y)$ obtained from ground state quantum Monte Carlo simulations for a simulation cell with $L_x\times L_y\times L_z=14.75707\times 17.04\times 5.05\text{Å}$ corresponding to the fully filled phase with $f=1$ for $N=48$ ${}^4\mathrm{He}$ atoms on adsorption sites. Finite size effects in the spatial wave function are negligible beyond nearest-neighbor ($V$) and next-nearest-neighbor ($V^{\prime }$) couplings in the effective Bose-Hubbard description are indicated with solid and dashed lines, respectively.

Figure 12

The effective nearest-neighbor interaction parameter of the Bose-Hubbard model computed from Eq. (29) via quantum Monte Carlo for simulation cells with and 48 as a function of the cell size in the $z$ direction, $L_z$. The indicated value of $V_{\mathrm{QMC}}=54.3\left(1\right)\mathrm{K}$ was computed at $L_z=5.05\text{Å}$ as described in the text. The semi-transparent symbols for $L_z\ge 5.25\text{Å}$ indicate cell sizes which allowed the nascent formation of a second layer, where the mapping of the microscopic Hamiltonian to the 2D Bose-Hubbard model breaks down.

Figure 13

Convergence of the kinetic (a), adsorption (b), and interaction potential (c) energy per particle with the imaginary time projection length $\beta$ and imaginary time step $\tau$ (insets) for a system with adsorption sites at unit filling. All quantities measured via quantum Monte Carlo (symbols $+$ errorbars) are observed to converge to their ground state value with the expected exponential dependence on $\beta$ or polynomial dependence on $\tau$ as quantified via the indicated fits (lines). $\beta$ scaling was performed with $\tau =1/319\mathrm{K}$ (vertical dashed line in insets) while $\tau$ scaling had $\beta \simeq 0.5{\mathrm{K}}^{-1}$. A subscript 0 indicates the ground state value.

Figure 14

The finite size scaling of the next nearest-neighbor interaction $V^{\prime }$ extracted from quantum Monte Carlo as a function of the number of graphene hexagon adsorption sites. The solid line shows a linear fit to the data for which allows extrapolation to the thermodynamic limit.

Figure 15

The density profile (per particle) of the adsorbed layer(s) for different vertical box sizes $L_z$ enforced through the potential in Eq. (25) for helium above a graphene sheet with adsorption sites such that $L_x\times L_y=14.75707\text{Å}\times 17.04\text{Å}$. Panels correspond to (a) filling $f=1$ and (b) $f=1/3$ where statistical uncertainties are indicated by the shaded envelope. The thicker curve in (a) for $L_z=5.05\text{Å}$ was determined to be the optimal value (see text). The dashed line indicates the density profile for a “bulk” cell with $L_z=10\text{Å}$ at $f=1/3$ that is used for comparison. The number of ${}^4\mathrm{He}$ atoms in the simulation can be determined from

Figure 16

The squared deviation between solid and dashed curves in Fig. 15 as a function the box size in the $z$ direction as quantified in Eq. (B1). The minimum at $L_z=5.05\text{Å}$ is independent of the size of the graphene sheet, where data for and 48 are shown. The inset shows a comparison of the density profiles in the $z$ direction for this value at . The dashed line is a guide to the eye.

Figure 17

Equation of state (energy per particle as a function of filling fraction) for adsorption sites for two box sizes with $L_z=5.05\text{and}10.0\AA$. The curves have been shifted by the energy for a single particle $N=1$ corresponding to a filling fraction-independent value of $\sim 6K$ due to the presence of ${\mathcal{V}}_{\mathrm{wall}}$. The insets show the density of particles along the $z$ direction at filling fractions $f=1/3$ and 1. For unit filling, the cell with $L_z=10\text{Å}$ can accommodate a second layer.

Figure 18

A one-dimensional cut through for a ${}^4\mathrm{He}$ atom above graphene corresponding to the spatial path through the saddle point of the full 2D potential shown in Fig. 8. The energy barriers that inhibit tunneling between different graphene adsorption sites originate from three different methods: Møller-Plesset perturbation theory (MP2, green), density functional theory (DFT, blue), and quantum Monte Carlo (QMC, red) as described in the text, where the curves have been shifted with respect to the bottom of the potential well. The classical turning points $x_{c_{1,2}}$ for each potential barrier appear at the intersection of the dashed ($E_0=\hslash {\omega }_0/2$) and solid lines.