1. Introduction
The physics of thin helium films varies greatly, depending on the substrate upon which they are adsorbed [1]. The equilibrium phase of 4He in two dimensions (2D) is a superfluid liquid, and indeed on the weakest substrate for which a stable 2D film forms, namely lithium, a stable superfluid monolayer is predicted [2] and observed [3], with continuous growth of film as a function of the chemical potential [4]. On the other hand, on stronger substrates such as graphite, the phase diagram features a dazzling variety of phases, including crystalline ones, either commensurate or incommensurate with the underlying substrate.
An interesting theoretical question is whether a phase simultaneously displaying diagonal (i.e. crystalline) and off-diagonal (i.e. superfluid) long-range order may occur at sufficiently low temperature, in some well-defined range of 4He coverage. Some experimental claims have been made of concurrent superfluidity and crystalline order in the second layer of 4He adsorbed on graphite, specifically at, or in the vicinity, of a solid phase registered with the underlying first (solid) layer [5–8]. Though the denomination supersolid is, strictly speaking, not applicable to a system of this type [9, 10], the observation of such a phase would nonetheless be of great significance, as it would help elucidate the interplay between these two seemingly antithetic types of order. It seems fair to state, however, that at this time all the evidence is indirect and/or inconclusive; indeed, the interpretation of the findings has been questioned and is not supported by reliable microscopic theoretical calculations [11–15].
Graphene, which consists of a sheet of graphite, is also regarded as a potentially interesting adsorber. A single layer of carbon atoms is about 10% less attractive than graphite; the ensuing enhancement of the quantum excursions of helium atoms away from the basal plane acts to soften the repulsive interaction among them at short distances. All of this may ultimately increase the mobility of 4He, possibly stabilizing novel phases displaying a more quantal character (e.g. superfluid phases with an unusual degree of atomic localization). Early microscopic calculations for helium on graphite showed that a mere 10% reduction in the strength of the attractive part of the interaction potential results in no significant change in the overall physical behavior of an adsorbed helium film [11], a conclusion confirmed by subsequent calculations for 4He on graphene. In particular, there is general agreement on the existence of the same commensurate crystalline phase of a 4He monolayer (known as C1/3) which forms on graphite. Indeed, the basic features of the phase diagram of an adsorbed 4He monolayer on graphene can be understood on the basis of a highly simplified description of the system, making use of a lattice Hamiltonian of hard core bosons [16].
There is, however, some controversy concerning a possible finite superfluid response of such a crystalline phase, at sufficiently low temperature, as well as on the existence of superfluid phases at lower coverage, either liquid or crystalline (doped with vacancies).
Finite-temperature calculations have yielded no evidence of a nonzero superfluid response of the C1/3 phase [17, 18], while the claim of a small superfluid signal in the ground state has been made [19]. Furthermore, some studies [18, 20] have given evidence of a possible low-lying (metastable) superfluid phase at sufficiently low temperature, one that could conceivably be long-lived, and therefore observed experimentally.
In this work, we revisit the issue of superfluidity at low temperature of a monolayer of 4He absorbed on a single graphene sheet, at and below commensurate coverage, by means of quantum Monte Carlo simulations at zero and finite temperature. We make use of a microscopic model of the system that is consistent with most previous studies, explicitly accounting for the corrugation of the graphene substrate. We monitor the superfluid response and the emergence of diagonal long-range order at low temperature.
Altogether, our results confirm the absence of superfluidity in this system at commensurate coverage, all the way to temperature T = 0. No inconsistency is observed between finite-temperature and ground-state simulations, i.e. even the latter yield no evidence of any finite superfluid response. In particular, our ground state estimates result from projecting the lowest-energy component from an initial trial wave function that does not break translational invariance, i.e. with no built-in information about the presence of a corrugated substrate. Yet, crystalline long-range order quickly builds in, even for relatively short projection imaginary-time intervals, disproving the contention put forth in [20] about the existence of a superfluid state energetically competitive with the crystalline ground state, based on a similar projection procedure.
We find no evidence of a thermodynamically stable, uniform superfluid phase for coverage below commensuration; indeed, our results are in disagreement with the contention made in [17] that structural and superfluid properties of the system below commensuration are greatly affected by the choice of potential describing the interaction of a helium and a carbon atom, specifically an isotropic form of the potential leading to a superfluid uniform phase below commensurate coverage. On the contrary, we find the physical picture to be the same as that predicted with an anisotropic potential, with commensurate crystalline order persisting in the presence of vacancies, and with the formation of solid clusters at low coverage, analogously to what has been predicted to occur on graphite. Finally, we see no evidence of any metastable superfluid phase at commensurate coverage.
This manuscript is organized as follows: we introduce the microscopic model of the physical system for which we have carried out our calculations and describe the computational methodology adopted in this work in section 2 , present our results in detail in section 3 , and outline our conclusion in section 4 .
2. Model and methodology
We consider an assembly of N 4He atoms, regarded as point-like identical particles of mass m and spin S = 0, moving in three physical dimensions. The system is enclosed in a cell of volume Ω = Lx × Ly × Lz, with Lx = 34.08 Å, Ly = 29.514 Å, and Lz = 40 Å. Periodic boundary conditions are used in the three directions (but the boundary condition along the z direction is immaterial, as helium atoms remain confined in the vicinity of the substrate). The cell sizes are chosen to accommodate an integer number of graphene unit cells in both directions, without introducing frustration. The graphene sheet is modeled as an ideal, 2D (honeycomb) lattice; our simulation cell accommodates Nc = 384 carbon atoms at fixed positions on the z = 0 plane, with a carbon-carbon bond length a = 1.42 Å. Regarding the carbon atoms as pinned at classical lattice sites is an assumption that has been made in all other comparable studies, and one that seems justified by their relatively large mass, compared to that of the helium atoms. The 4He coverage (i.e. 2D density) is θ = N/A, where A = Lx × Ly.
The quantum-mechanical many-body Hamiltonian reads as follows:
$\begin{eqnarray}\hat{H}=-\displaystyle \sum _{i}\left[\lambda {{\rm{\nabla }}}_{i}^{2}+U({{\boldsymbol{r}}}_{i})\right]+\displaystyle \sum _{i\lt j}v({r}_{{ij}}),\end{eqnarray}$
where the first (second) sum runs over all particles (pairs of particles), λ ≡ ℏ2/2m = 6.05964 KÅ2, rij ≡ ∣ri − rj∣, v(r) denotes the pairwise interaction between two 4He atoms and $\begin{eqnarray}U({\boldsymbol{r}})=\displaystyle \sum _{\alpha =1}^{{N}_{c}}u(| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{\alpha }| ),\end{eqnarray}$
the sum running over all carbon atoms and u being a potential describing the interaction between a helium and a carbon atom. Both v and u are assumed to depend only on the relative distance between two particles; this is a standard assumption for the helium pair potential v, justified by the essentially spherical shape of a helium atom in its ground state. For the helium–carbon part, there is still uncertainty on whether a spherically symmetric potential u is adequate.Because the individual carbon atoms are explicitly included in equation (1 ), the model is expected to capture the most important contribution of the corrugation of the graphene substrate, which is necessary to stabilize a commensurate crystalline layer [21]. However, Kwon and Ceperley [17] have contended that the physics of the system below the commensurate C1/3 coverage θ0 = 0.0636 Å−2 is significantly affected by the choice of u; in particular, a superfluid phase (either liquid or crystalline, doped with vacancies) exists if an isotropic form of the potential u is used, while, on using an anisotropic form, no thermodynamically stable phase is found for coverage θ < θ0. This assertion is at variance with the subsequent findings of Happacher et al [18], who carried out grand canonical simulations using an isotropic potential u and found no evidence of any thermodynamically stable phase of coverage less than θ0.
The choice made here for u and v is motivated both to ensure consistency with (most of) the existing calculations, as well as to clarify some of the above, still unresolved aspects of the physics of the system for the simplest choice of helium–carbon interaction, i.e. that based on an isotropic form for u. We adopt therefore the same version of the Aziz potential utilized in [19, 20], though it should be noted that the differences between the various versions of the Aziz potential are quantitatively small and that the helium-graphene interaction dominates the energetics in this system. The carbon–helium pair-wise interaction u is a Lennard–Jones potential with parameters ε = 16.25 K and σ = 2.74 Å taken from [22], again for consistency with previous calculations. In our simulation, the potentials are cut off and smoothly shifted to zero at a distance rc = 14.75 Å; we estimate the energetic contribution arising from interactions between particles located at greater distances to amount to 0.36 K per 4He atom.
The low-temperature phase diagram of the thermodynamic system described by equation (1 ) as a function of coverage and temperature has been studied in this work by computer simulations at zero and finite temperature. In the majority of our calculations the coverage θ = θ0, i.e. it is N = 64. For this coverage, the C1/3 crystalline monolayer 4He film commensurate with the underlying graphene substrate is expected to form. We also carried out some simulations at lower coverage, to investigate the possible occurrence of a thermodynamically stable superfluid phase, either liquid or a vacancy-rich commensurate crystal.
Finite-temperature simulations have been carried out in the temperature range 0.25 ≤ T ≤ 4 K. They are based on the continuous-space Worm algorithm [23, 24]. Since this technique is well-established, and extensively described in the literature, we shall not review it here; we used the canonical variant of the algorithm, in which the number of particles N is fixed [25, 26]. The most important aspects of this methodology, in the context of the study carried out here are that it allows one to obtain unbiased estimates of energetic and structural properties of the system, essentially with no approximation, as well as of the superfluid fraction, based on the winding number estimator [27].
Ground state simulations, on the other hand, were performed using the related Path Integral Ground State (PIGS) method [28], which projects the true ground state out of an initial trial wave function. Details of both simulations are standard; in both cases, we made use of a fourth-order approximation for the short imaginary time (τ) propagator (see, for instance, [29, 30]), and all of the results presented here are extrapolated to the τ → 0 limit. We generally found numerical estimates for the physical properties of interest here obtained with a value of the time step τ ∼ 1.5 × 10−3 K−1 to be indistinguishable from the extrapolated ones, within the statistical uncertainties of the calculation (the only exception is the energy in ground state calculations, for which a smaller value of the time step was needed. We discuss this aspect in greater depth below).
3. Results
As a preliminary comment, which applies to both finite-temperature and ground-state results presented here, it needs to be mentioned that in all cases they are independent of the initial atomic configuration utilized in the simulations.
3.1. Finite-temperature results
We begin with a discussion of the low-temperature energetics of the system. Figure 1 shows the computed energy per 4He atom e(T) as a function of temperature, at the commensurate coverage θ0. As shown in figure 1, the estimates are unchanged below T = 1 K, i.e. they can be regarded as ground-state results. The data can be fitted with the expression e(T) = e0 + αT3, yielding e0 = −129.550(25) K. The first observation is that energy estimates at finite temperature (for T ≤ 1 K), for the particular coverage considered here are significantly lower (by about 0.25 K) than the ground state estimate quoted in [20], based on a calculation making use of the same potentials utilized in this work.1(1For comparison, we have carried out the same calculation with an earlier version of the Aziz helium pair potential [31] and obtained a ∼0.05 K higher ground state energy estimate. No significant change is observed in any of the other physical observables considered in this work if the different helium potential is used.)
Figure 1. Energy per 4He atom (in K) computed by simulation as a function of the temperature at a coverage θ = 0.0636 Å−2. Statistical errors are smaller than symbol sizes. The solid line represents a fit to the data obtained with the expression e(T) = e0 + αT3. The estimate shown with a diamond at T = 0 is that obtained with PIGS (see section |
It is worth noting that the difference between our estimate and that of [20] is three times greater than that between the two results provided in [20], projecting the lowest-energy state out of two initial trial wave functions, one possessing crystalline long-range order and the other not breaking translational invariance. It was argued therein that the relatively small energy difference arising from the calculation carried out with the two different trial wave functions pointed to the existence of a low-lying superfluid phase, energetically competitive with the insulating, crystalline ground state. The results obtained in this work suggest instead that the (Diffusion Monte Carlo, DMC) projection utilized in [20] actually failed to reach convergence to the ground state with either trial wave function.2(2It seems quite likely that the statistical errors have been underestimated in [20], i.e. no statistically significant difference exists between the energy estimates obtained with the two trial wave functions.) The same observation was made in the past [32], pointing out how path integral approaches, which are not affected by population control bias [33–35] are typically a more reliable option than DMC, when it comes to investigating the ground state of Bose systems.
The occurrence of crystalline long-range order at low temperature at coverage θ0 can be quantitatively established by the calculation of the static structure factor S(q), shown in figure 2 (left). The occurrence of a sharp (Bragg) peak at q = 1.7 Å−1 is what one expects for the commensurate ($\sqrt{3}\times \sqrt{3}$) commensurate crystal. The right part of figure 2 shows instead the average 4He density profile as a function of the distance from the graphene plane. The results shown for both quantities are for a temperature T = 0.5 K, but they are found to be essentially temperature-independent below T = 2 K.
Figure 2. Left: static structure factor S(q) for an adsorbed film of 4He on graphene, computed by simulation at temperature T = 0.5 K, for wave vectors all lying in the graphene plane. Right: average 4He density as a function of the distance from the graphene plane, computed at the same temperature. |
No evidence of superfluidity is seen at the lowest temperature considered here, namely 0.5 K, with exceedingly infrequent atomic exchanges, never extending beyond a four-atom ring permutation. Obviously, however, the onset of superfluidity should be expected at a considerably lower temperature than 0.5 K. Indeed, the estimate of the superfluid fraction proposed in [20] would suggest a superfluid transition temperature of the order of 5 mK, based on the well-known ‘universal jump’ condition [36]. We come back to a critical re-examination of this, and other claims made in [20] below, when discussing our ground state results. For the moment, we simply note the paucity of atomic exchanges, which instead occur quite frequently at this temperature in 2D, as well as that a previous study of this system, based on the same computational methodology utilized here yielded no evidence of superfluidity of the C1/3 phase down to temperatures as low as 10 mK [18].
3.1.1. Low coverage results
In [17] the claim was made that, upon using an isotropic helium–carbon potential in model (1 ), one could observe a (metastable) superfluid phase at coverage below θ0, typically by introducing a few vacancies, or by reducing the 4He coverage to approximately 0.058 Å−2. Our study yielded no evidence of any low-coverage superfluid liquid phase. Specifically, our simulations at coverage θ ≤ θ0 unambiguously show crystalline order, consistently with the coexistence of a crystal (of coverage θ0) and vapor, which is in agreement with the findings of [18]. The claim of a liquidlike phase of [17] is based on the observation of a sudden drop of the peak of the static structure factor, on lowering the coverage from θ0 to approximately 0.9 θ0, and from the visual inspection of the configuration generated in the course of a Monte Carlo simulation, suggesting loss of local crystalline order at the lower coverage. On the other hand, if an anisotropic pair potential is utilized the static structure factor does not drop as significantly, and local crystalline order is retained. The results of our simulations show that, while, on the one hand, the anisotropic potential can plausibly be more effective at ‘pinning’ helium atoms in place than the isotropic one, the main physical behavior is actually unchanged in the case of an isotropic potential. Figure 3 shows a comparison of the static structure factor computed at T = 1 K (the main results and observations do not change at lower temperatures, down to the lowest considered here, namely T = 0.25 K) for the case of coverage θ0 (circles) and θ = 0.0577 Å−2 (diamonds). The peak of S(q) is lower at the lower coverage (though the difference is nowhere near as large as that reported in [17]), but this is insufficient to conclude that no crystalline order is present at the lower coverage. For example, the calculation of [17] at coverage θ0 yields a ∼50% higher peak of S(q) if an anisotropic He–C pair potential is used, but in both cases the crystalline character of the monolayer is not in question. Furthermore, although indeed the main peak is depressed by approximately a third (left part of figure 3), nevertheless it remains relatively sharp, suggesting that crystalline order is largely retained, as confirmed by snapshots of instantaneous many-particle world line configurations (right part of figure 3).
Figure 3. Left: static structure factor S(q) for an adsorbed film of 4He on graphene, computed by simulation at temperature T = 1 K, for wave vectors all lying in the graphene plane, at coverage θ0 = 0.0636 Å−2 (circles) and θ = 0.0577 Å−2 (diamonds). Statistical errors are smaller than symbol sizes. Right: snapshot of many-atom configuration (particle world lines) corresponding to the lower coverage. Solid circles represent the carbon atoms. |
It is contended in [17] that, at coverage less than θ0, permutations of identical helium atoms begin to occur, leading at sufficiently low temperature to a finite superfluid response, which for sufficiently low coverage is interpreted as the stabilization of a superfluid liquid-like state. In our simulations we also observe such permutations, in fact even at temperatures as high as 1 K (as shown in figure 3). They typically involve atoms that lack some nearest neighbors, and in a small simulation cell cycles involving several particles occasionally appear, connecting two opposite sides of the cell and resulting in a significant superfluid signal, which, for relatively small system sizes (e.g. ∼30 atoms) may appear as a bulk superfluid response. This is a finite-size effect, however; in the thermodynamic limit coexistence of a crystal with a low-density vapor is expected. In all of our simulations, the estimate of the superfluid fraction obtained as an average over sufficiently long runs is zero within statistical errors (i.e. the average value is typically much smaller than the statistical error). This is true both at commensurate and at lower coverage.
Altogether, therefore, although the anisotropic He–C potential has a stronger pinning effect on the helium atoms, nonetheless the basic physics of the system below commensurate coverage is the same, regardless of whether an isotropic or an anisotropic potential is used. The physical picture that emerges from our simulations below commensurate coverage is the formation of solid 4He clusters, analogously to what is observed on graphite [37].
3.1.2. Search for metastable superfluid state
In order to explore the possible existence of a metastable, long-lived liquid-like superfluid phase of coverage θ0 we have carried out the same ‘computer experiment’ of [38, 39]. Specifically, we initially prepared the system in a superfluid phase by simulating a 4He monolayer adsorbed on a Li substrate (using the same simulation cell and number of atoms). The Li substrate is modeled as flat (i.e. smooth, featureless); on it, at a temperature T = 0.5 K, a 4He monolayer at coverage θ0 forms a nearly 2D superfluid, with a superfluid fraction close to 100% [2]. Upon establishing that thermal equilibrium was reached, the lithium was replaced with the graphene substrate, i.e. we reverted to the Hamiltonian equation (1 ), to assess the resilience of the superfluid phase. One would expect that if a low-lying metastable superfluid phase existed, the simulation algorithm may remain ‘stuck’ in it, finding the true equilibrium (crystalline) phase requiring the disentanglement [40] of the many-particle paths present in the low-temperature superfluid phase stabilized on the Li substrate.
But what is observed, instead, is that as long as the all-important ‘swap’ moves [24] are attempted sufficiently often, the replacement of the Li substrate with the graphene causes the existing long atomic exchange cycles to disappear rather quickly, with a corresponding decrease of the initial superfluid signal. As most long permutation cycles have disappeared, the system spontaneously begins to develop solid order. All of this is illustrated in figure 4, which displays instantaneous many-particle configurations (world lines) for the simulations conducted on a Li substrate (left part of figure 4) and on the graphene substrate (right part of the figure). On a Li substrate, particle world lines entangle, as expected of a superfluid, and even telling atoms apart is a difficult proposition; on the graphite substrate, on the other hand, after a relatively short time individual atoms begin to be clearly identifiable, as exchanges become suppressed and a regular (lattice) arrangement of atoms emerges. Concurrently, the superfluid fraction decreases (it is essentially zero, within statistical fluctuations, by the time the arrangement shown on the right side of figure 4 appears).
Figure 4. Snapshots of many-particle configurations (world lines) corresponding to simulations of 4He on a flat Li substrate (left) and on the graphite substrate (right, small circles represent the C atoms). Both simulations are at T = 0.5 K. The one on Li is such that the system is entirely superfluid, within the statistical errors of the calculation. The simulation on graphene is restarted from a configuration for the system on the left, after thermodynamic equilibrium is reached, simply with the replacement of the substrate. |
It is possible, and we have observed this in our simulations too, to have at times few isolated and resilient permutation cycles, which in a relatively small system can wind around the periodic boundaries, giving rise to a finite superfluid signal. We attribute to this effect the observation of possible metastable superfluid phases in the grand canonical simulations of [18]. However, such a superfluid signal is spurious and disappears in the thermodynamic limit. A useful diagnostic tool is the calculated histogram of the frequency P(n) of occurrence of exchange cycles including up to n 4He atoms. In the superfluid phase, this is a smoothly (exponentially) decaying function. On the other hand, if P(n) = 0 for n greater than a few (with the possible exception of isolated peaks at specific numbers, merely signaling the inability of the underlying sampling procedure to remove all permutation cycles), one is in the presence of a (crystalline) insulator. Our results clearly point to the latter scenario, i.e. they do not support the existence of a low-lying superfluid liquid-like phase, at variance with the contention of [20]. We re-examine critically such a contention when discussing our ground state results in section 3.2 .
3.2. Ground state results
In the case of ground state simulations, we restricted ourselves to the case of commensurate coverage θ0. The study of the ground state of the system described by equation (1 ) was carried out using the following trial wave function:4 ) are α = 19 and β = 0.12 Å−5.
$\begin{eqnarray}{{\rm{\Psi }}}_{T}({{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2},\ldots {{\boldsymbol{r}}}_{N})=\displaystyle \prod _{i\lt j}\ \exp \left[-w({r}_{{ij}})\right],\end{eqnarray}$
where $\begin{eqnarray}w(r)=\displaystyle \frac{\alpha }{1+\beta {r}^{5}}.\end{eqnarray}$
In the PIGS method [29], the true ground state Φ0 is approached by projecting it out of the initial trial wave function $\Psi$T as Φ0 = lim${}_{{\rm{\Lambda }}\to \infty }{\rm{\Phi }}({\rm{\Lambda }})\equiv {\mathrm{lim}}_{{\rm{\Lambda }}\to \infty }{{\rm{e}}}^{-{\rm{\Lambda }}\hat{H}}{{\rm{\Psi }}}_{T}$. The values of the variational parameters in equation (The most important aspect of this trial wave function is that it only includes pairwise correlations among helium atoms, accounting for the presence of a repulsive core at short interparticle distances in the interatomic potential v. Although it would be in principle possible to construct a trial wave function with additional terms reflecting the presence of a corrugated substrate, in a way that retains the intrinsic indistinguishability of helium atoms [41], the wave function utilized here has no built-in information concerning the presence of an attractive substrate. In other words, if one were to attempt a variational calculation the helium atoms would not remain close to the substrate, i.e. the whole system would simply evaporate. At the variational level, this is remedied by including, for example, a one-body term, confining helium atoms to the proximity of the substrate. In this work, we chose not to include any such term, relying instead on the projection algorithm to stabilize a 4He monolayer; we find that if Λ ≥ 0.125 K−1 a stable monolayer forms with no detectable atomic evaporation, and the 4He density profile in the direction perpendicular to the substrate is essentially indistinguishable from that computed at finite temperature, shown in figure 2.
The reason for using a trial wave function of such a simple form is that it does not break translational invariance, and thus, by definition, it represents a fully superfluid state. Since we aim to explore the possible presence of a finite superfluid signal in the ground state, we choose as the starting point of the projection a state that embodies the most favorable conditions for superfluidity, searching for a residual superfluid signal in the ground state. As the true ground state breaks translational invariance due to the presence of the corrugated substrate, it cannot be fully superfluid [42].
As stated in section 2 , we have observed numerical convergence of the estimates for all physical quantities, within statistical errors, for a value of the time step equal to that used in finite temperature simulations, namely τ = 1.5 × 10−3 K−1. For the ground state energy, on the other hand, we found that the value of the time step needed to observe convergence is as much as 10 times smaller. Our result for a projection time Λ = 1 K−1 is shown in figure 1, and is consistent with the finite temperature estimate, within the statistical errors of the calculation (it is worth noting that the energy estimate yielded by PIGS for a finite projection time is a strict upper bound on the true ground state energy).
At this point, a general remark is in order. The ground state wave function of a Bose system, for the case of a Hamiltonian not breaking time-reversal symmetry, is real and positive. Consequently, a projection algorithm such as PIGS or DMC necessarily must converge to the true ground state, given a sufficiently long projection time, upon starting from a positive-definite trial wave function. One could imagine, in the case in which the ground state were, e.g. crystalline, that the projection time required might become exceedingly long if one started from a liquid-like trial wave function, and if there existed an excited state not possessing crystalline order, with energy very close to that of the ground state. One way to explore such an occurrence is by computing relevant quantities, indicative of the presence of order, as a function of projection time.
Figure 5 shows the static structure factor S(q) computed by PIGS for three different values of the projection time Λ. For a relatively short Λ, the projected state Φ(Λ) retains the physical properties of the trial wave function $\Psi$T, including the lack of crystalline order. However, as Λ increases a sharp peak in correspondence of the $\sqrt{3}\times \sqrt{3}$ crystalline order appears, with the height of the peak monotonically growing with Λ.
Figure 5. Static structure factor computed by PIGS for three different projection times, for wave vectors in the plane of the substrate, showing the emergence of a sharp (Bragg) peak as the projection time is increased, signaling the onset of $\sqrt{3}\times \sqrt{3}$ crystalline order. |
In other words, the projected state Φ(Λ) develops crystalline long-range order irrespective of the fact that no such order is present in the trial wave function. This is an all-important point, as it invalidates the contention made in [20] that a low-lying liquid state very close in energy to the ground state exists, based on the presumption that a state projected out of a trial wave function not possessing crystalline order should retain the same physics, i.e. remain disordered. The result of figure 5 clearly shows that this is not the case; on the contrary, crystalline order is robust at low temperature, in this system, and no energetically competitive fluid state exists, a conclusion consistent with our findings at finite temperature.
We now discuss the possible presence of a small but finite superfluid response in the ground state of the system at coverage C1/3, to address the claim made in [19] of a superfluid fraction ρS = 0.0067 at zero temperature. No evidence of superfluidity has been seen in any finite temperature study, including this one [17, 18], and even though the superfluid transition temperature can be expected to be relatively low, finite temperature simulations of small-size systems approaching that temperature failed to yield even a hint of developing superfluid order [18].
Within a ground state simulation method like PIGS, the superfluid fraction ρs can be obtained from the long imaginary-time diffusion of particle world lines [43], specifically as5 ) trivially yields a value of 1, but if translational invariance is broken, as is the case in equation (1 ), then ρS < 1 at T = 0.
$\begin{eqnarray}\begin{array}{rcl}{\rho }_{S} & = & {\mathrm{lim}}_{t\to \infty }D(t),\ \mathrm{with}\\ D(t) & = & \displaystyle \frac{N}{2d\lambda }\displaystyle \frac{\langle {\left[{{\boldsymbol{R}}}_{\mathrm{CM}}(t)-{{\boldsymbol{R}}}_{\mathrm{CM}}(0)\right]}^{2}\rangle }{t}\end{array}\end{eqnarray}$
where ⟨...⟩ stands for average value, d is the dimensionality (in this case d = 2 as the system is essentially two-dimensional) and with ${{\boldsymbol{R}}}_{\mathrm{CM}}(t)=(1/N){\sum }_{i\,=\,1}^{N}{{\bf{r}}}_{i}(t)$, ri(t) being the position of the ith particle at imaginary time t along the world line, with 0 ≤ t ≤ Λ. For a system that enjoys translational invariance, the estimator equation (Since the projection time Λ is finite, the limit in equation (5 ) must be carried out by extrapolation, based on a fit of the diffusion curve D(t) for a finite Λ interval, in the t → Λ limit. We used here the formula proposed in [43], i.e.
$\begin{eqnarray}D(t)={\rho }_{S}+\displaystyle \frac{a}{t}\left(1-\exp (-{bt})\right),\end{eqnarray}$
where a, b and ρS are fitting parameters, ρS representing the extrapolated superfluid fraction at infinite projection time. For a sufficiently large projection time Λ, one expects the values of the fitting parameters not to change, within their statistical uncertanties.We carried out such a fitting procedure by means of a Metropolis random walk through parameter space using the value of χ2 as the ‘energy’ function; this approach offers access to the probability distributions of the values of the fitting parameters and allows for a reliable estimation of their uncertainties. Figure 6 illustrates the result of this procedure for ρS in correspondence of a projection time Λ = 4 K−1.
Figure 6. Left: fit to the diffusion curve at projection time Λ = 4 K−1 using equation ( |
Although the expression equation (6 ) aims at providing a good fit to the data only in the t → Λ limit, it turns out to be possible to obtain an acceptable fit, essentially over the whole imaginary time range, as shown by the left side of figure 6. The probability distribution for the value of the parameter ρS of the fit, shown on the left side of figure 6 places our estimate of the superfluid fraction in the ground state at less than ∼0.005 with 99.99% confidence, which is already well over an order of magnitude lower than the estimate of [19]. We repeated this analysis for different values of Λ but we never obtained a lower bound for ρS greater than zero, with an acceptable statistical confidence. Thus, we conclude that our results are consistent with a null value of ρS in the ground state, in accord with all finite temperature studies. Because this conclusion is established on a finite system, a fortiori it must hold the thermodynamic limit.
Once again, we wish to point out that by choosing a trial wave function that is translationally invariant we started out from the most favorable conditions for superfluidity. The fact that, given a sufficiently long projection time, the superfluid signal becomes essentially not measurable constitutes in our view strong evidence against the presence of superfluidity in the ground state. We attribute the disagreement between this conclusion and that of [19] to likely remnant variational bias in the DMC projection, quite likely arising from the finite population size, a problem already extensively documented in the literature [35].
4. Conclusion
We investigated the physics of a 4He monolayer adsorbed onto graphene substrate via quantum Monte Carlo simulations both at zero and finite temperature. The goals of this study were to clarify the possible presence of (metastable) superfluid phases at commensurate and below commensurate (C1/3) coverage and to clarify outstanding discrepancies between results and predictions yielded by different numerical studies.
The main conclusions of our study are that there are no (metastable) superfluid phases of the system at low temperature, either at C1/3 coverage or below. In particular, our study lends no support to the contention, put forward in [17], that the superfluid properties of the system below commensurate coverage are significantly affected by the choice of potential describing the interaction between a helium and a carbon atom (specifically, an isotropic potential stabilizing a superfluid phase below commensurate coverage). We find that the physics of the system is essentially the same for both types of potentials, and mimics that on a graphite substrate, i.e. no vacancy-induced supersolidity [44] at low doping, with the formation of solid clusters at coverages below C1/3.
Our study also confirms the absence of any superfluid signal at commensuration, in agreement with all similar studies and with general arguments to the effect that no superfluid signal is seen at commensuration [45]. Indeed, it is now accepted that only in the presence of multiple occupation of the unit cell (i.e. in cluster crystal) can supersolidity occurs at commensuration [46]. We attribute the finite superfluid response reported at commensuration in [19] to the (now understood) limitations of the technology utilized therein.
This work was supported by the Natural Sciences and Engineering Research Council of Canada. The computer codes utilized to obtain the results can be obtained by contacting the authors. The authors declare no conflict of interest.