# In search of a spatially inhomogeneous strongly interacting superfluid:
Study of submonolayer 4He on Fluorographene, Hexagonal Boron-Nitride, and
Graphene

Saverio Moroni  CNR-IOM Democritos and SISSA, Via Bonomea 265, 34136 Trieste,
Italy  Francesco Ancilotto  Dipartimento di Fisica e Astronomia ’Galileo
Galilei’ and CNISM, Università di Padova, via Marzolo 8, 35122 Padova, Italy
CNR-IOM Democritos, Via Bonomea 265, 34136 Trieste, Italy  Pier Luigi
Silvestrelli  Dipartimento di Fisica e Astronomia ’Galileo Galilei’ and CNISM,
Università di Padova, via Marzolo 8, 35122 Padova, Italy  CNR-IOM Democritos,
Via Bonomea 265, 34136 Trieste, Italy  Luciano Reatto  Dipartimento di Fisica,
Università degli Studi di Milano,via Celoria 16, 20133 Milano, Italy

###### Abstract

We study a sub monolayer 4He adsorbed on fluorographene (GF) and on hexagonal
boron nitride (hBN) at low coverage. The adsorption potentials have been
computed ab-initio with a suitable density functional theory including
dispersion forces. The properties of the adsorbed 4He atoms have been computed
at finite temperature with path integral Monte Carlo and at $T=0$ K with
variational path integral. From both methods we find that the lowest energy
state of 4He on GF is a superfluid. Due to the very large corrugation of the
adsorption potential this superfluid has a very strong spatial anisotropy, the
ratio between the largest and smallest areal density being about 6, the
superfluid fraction at the lowest $T$ is about 55% and the temperature of the
transition to the normal state is in the range 0.5-1 K. Thus, GF offers a
platform for studying the properties of a strongly interacting highly
anisotropic bosonic superfluid. At a larger coverage 4He has a transition to
an ordered commensurate state with occupation of 1/6 of the adsorption sites.
This phase is stable up to a transition temperature located between 0.5 and 1
K. The system has a triangular order similar to that of 4He on graphite but
each 4He atom is not confined to a single adsorption site and the atom visits
also the nearest neighboring sites giving rise to a novel three–lobed density
distribution. The lowest energy state of 4He on hBN is an ordered commensurate
state with occupation of 1/3 of the adsorption sites and triangular symmetry.
A disordered state is present at lower coverage as a metastable state. In
presence of an electric field the corrugation of the adsorption potential is
slightly increased but up to a magnitude of 1 V/Å the effect is small and does
not change the stability of the phases of 4He on GF and hBN. We have verified
that also in the case of graphene such electric field does not modify the
stability of the commensurate $\sqrt{3}\times\sqrt{3}R30^{\circ}$ phase.

###### pacs:

67.25.bh,02.70.Ss,71.15.Mb

##  I Introduction

Bosons moving in a periodic external potential $V_{\rm ext}({\bf r})$ can be
found in two different states: depending on the amplitude of the modulation of
the external potential either the bosons are localized at the minima of
$V_{\rm ext}({\bf r})$ or the bosons are delocalized and superfluid at low
temperature. For very strong inter-boson interaction an incommensurate solid
can be present. These two regimes, localized or delocalized, have been
achieved with cold bosonic atoms moving in the periodic potential generated by
optical standing waves greiner_2002 . Also a submonolayer film of light bosons
like 4He atoms adsorbed on a crystalline substrate is expected to show one of
these two regimes depending on the character of the adsorption potential. In
practice so far the only substrate that can approach such ideal situation is
graphite because this material can be obtained with an extended almost perfect
surface on which one can study the properties of the adsorbed species
brunch_2007 . It turns out that the He-graphite adsorption potential is
characterized by a corrugation that is large enough so that the He atoms are
localized around the preferential adsorption sites and in fact experiment
brunch_2007  and theory pierce_1999 ; manousakis_recente  agree that the
ground state of a monolayer 4He is an ordered commensurate state in which the
4He atoms occupy one third of the adsorption sites so the ad atoms have a
crystalline triangular symmetry, the so called
$\sqrt{3}\times\sqrt{3}R30^{\circ}$ phase. Therefore 4He on graphite is non-
superfluid, superfluidity appearing only when at least two layers are present
reppy ; nyeki_2017 .

On theoretical basis a similar behavior is expected for the adsorption of 4He
on graphene (G) because its adsorption potential turns out to be very similar
to that of graphite. Other substrates commonly used in adsorption studies are
intrinsically disordered like that of a glass or are in any case too
disordered to be relevant on this issue. It would be of great interest to find
other materials with extended crystalline surfaces with a corrugation of the
adsorption potential smaller than that of graphite because this would give the
possibility of a new superfluid state of strongly interacting particles that
will be spatially anisotropic due to the influence of the substrate potential.
It has been of interest that recent theoretical studies nava_2012 ;
reatto_2013  found that for two derivatives of graphene, fluorographene (GF)
and graphane (GH), the adsorption potential is very different from that of
graphene and of graphite and it turns out that the 4He atoms are delocalized
and the ground state of monolayer 4He on GF and on GH was claimed to be
superfluid. That study was based on a semi empirical adsorption potential.
These results have been corroborated by a more recent study silvestrelli_2019
in which some of us have developed an adsorption potential based on ab intio
methods: even with such adsorption potential it was found that the 4He atoms
on GF are delocalized and the ground state of submonolayer 4He is a
superfluid. Interestingly, that work silvestrelli_2019  has shown that also
the monolayer of 4He on hexagonal-Boron Nitride (hBN) turned out to be a
superfluid.

Studies of Refs. nava_2012 ; reatto_2013 ; silvestrelli_2019  are based on
state of the art many body computations that should be able to provide exact
results for the system. However quantum simulations of particles in a highly
structured potential can be tricky due to the multiple energy scales that are
present and we decided to perform a new investigation of the adsorption of 4He
atoms on GF, hBN, and G. In the present study we tested the adsorption
processes by also adopting an alternative formulation of the ab initio theory.
Moreover, we have derived the adsorption potential in presence of an external
electric potential, that is a possible way to alter the corrugation of the
adsorption potential.

We have studied 4He on GF and hBN at $T=0$ K with the variational path
integral (VPI) method ceperley_1995 , also known as path integral ground state
(PIGS) sarsa_2000 , and at finite temperature with path integral Monte Carlo
(PIMC) method ceperley_1995 . In a VPI computation the quantum state is
obtained by projecting an assumed initial state with the imaginary time
evolution operator. For large enough propagation time $\beta$ one gets an
unbiased sampling rossi_2009  of the properties of the exact ground state of
the system if the initial state is not orthogonal to the ground state. How
large $\beta$ has to be must be found empirically in terms of convergence. Our
evidence is that the results in nava_2012 ; reatto_2013 ; silvestrelli_2019
were not fully converged. We are confident that the present VPI results are at
convergence also because they are in agreement, as expected, with those of
PIMC at low temperature.

In the case of 4He on GF we confirm the main result of the previous studies
nava_2012 ; silvestrelli_2019 : The ground state of 4He on GF is a superfluid
that has a strong spatial modulation. At variance with the result of nava_2012
; silvestrelli_2019 , we find that at coverage larger than that of the
equilibrium state there is a first order phase transition to a commensurate
state with triangular symmetry at a coverage corresponding to the occupation
of 1/6 of the adsorption sites.

On the other hand, in the case of 4He on hBN, we do not confirm the previous
result of the existence of a superfluid state and our results show that, in
the ground state and at low temperature, 4He atoms on hBN form a non-
superfluid commensurate state very similar to that of 4He on graphite.

For all the considered substrates we also find that the change of the
adsorption potential due to the application of an external electric field E is
not large enough to alter in a significant way the properties for E=0, at
least for the strength of E allowed by our approximations (see below).

The contents of the paper are as follows. In Section II we present the results
for the adsorption potential of He on GF, hBN, and G with and without an
external electric field. The quantum simulations of our study are described in
Section III. A summary and our conclusion are contained in Section IV.
Technical details of the quantum simulations are given in the Appendix.

##  II Potentials

Monolayer 4He films on GF and GH have been proposed recently as novel
superfluid systems characterized by strong in-plane anisotropies. This
remarkable prediction nava_2012 ; Nav12 ; Nav13 ; reatto_2013  was based on
quantum simulations where an essential ingredient is an accurate description
of the interaction between He atom and the substrates. Specifically, the He-
substrate interaction potential was modeled using a traditional semiempirical
approachSto80 , where the potential energy of a single He atom near the
surface is written as a sum of pair potential interactions made of a repulsive
part proportional to the local electron density (Hartree-Fock repulsion), and
an attractive part, in the form of a sum of damped He atom van der Waals (VdW)
interactions and polarization interaction due to the surface electric field
nava_2012 ; Nav12 ; Nav13 .

These effective potentials are known to be affected by quite large
uncertainties in the empirical coefficients used to model the interaction. The
importance of a precise knowledge of the adatom-surface interaction potentials
to make quantitative prediction on the adsorption properties of surfaces
cannot be overlooked: the wetting properties of rare-gas atoms on solid
surfaces, for instance, are known to strongly depend on the strength and
corrugation of the adatom-substrate potential.

For this reason, we decided to investigate from first principles the
interaction of He atoms with GF and hBN. using state-of-the-art Density
Functional Theory (DFT) functionals specifically designed to describe the weak
VdW interactions, with the goal of providing a more accurate description of
the interaction of He atoms with these surfaces. Recent applications of vdW-
corrected DFT schemes to the problem of atoms/molecules-surface interactions
have proven the accuracy of such methods in the calculation of both adsorption
distances and adsorption energies, as well as the high degree of its
reliability across a wide range of adsorbates. In particular, we have computed
the He atom adsorption energies on different surface sites and the potential
energy corrugations along the plane, which are the most crucial ingredients
for accurate quantum simulations of the adsorption of 4He. Our calculations
have been performed with the Quantum-ESPRESSO ab initio packageESPRESSO . A
single He atom, per supercell, is considered and we model the substrates
adopting periodically repeated orthorhombic supercells, with a $4\times 2$
structure, in the case of G with 64 C atoms, in the case of GF with 32 C atoms
plus as many F atoms, while in the case of hBN the substrate is formed by 16 B
and 16 N atoms. The lattice constants correspond to the equilibrium state of
the substrates. Repeated slabs were separated along the direction orthogonal
to the surface by a vacuum region of about 24 Å to avoid significant spurious
interactions due to periodic replicas. The Brillouin Zone has been sampled
using a $2\times 2\times 1$ $k$-point mesh. Electron-ion interactions were
described using ultrasoft pseudopotentials and the wavefunctions were expanded
in a plane-wave basis set with an energy cutoff of 51 Ry.

The calculations have been performed by adopting the rVV10Sabatini  DFT
functional (this is the revised, more efficient version of the original VV10
schemeVydrov ), where vdW effects are included by introducing an explicitly
nonlocal correlation functional. rVV10 has been found to perform well in many
systems and processes where vdW effects are relevant, including several
adsorption processesSabatini ; psil15 ; psil16 . This DFT functional is able
to well reproduce the reference structural data of GF, hBN, and G including
the “buckling displacement” in GFnava_2012 ; Nav12 ; Nav13 .

We have also investigated the effect of the application of an external uniform
electric field. Since in supercell calculations periodic boundary conditions
are imposed on the electrostatic potential, an external electric field is
simulated by adding a saw-like potential to the bare ionic potentialNeugebauer
, also including a dipole correction, according to the recipe of
BengtssonBengtsson . This represents the proper way to simulate an external
electric field in surface calculations with a slab geometry, provided the
electrostatic potential discontinuity falls in the middle of the vacuum
regionNeugebauer ; Bengtsson ; Meyer . A positive value means that the
external electric field points away from the surface in the positive $z$
direction (the adsorbed He atom is located in the positive $z$ region). The
additional external potential also leads to changes in the total energy of the
system and the Hellmann-Feynman forcesMeyer .

An upper limit for the amplitude of the external electric field, which can be
considered in the simulations, exists, that is determined by the thickness of
the vacuum region. In fact a kind of “quantum well” is formed in the vacuum
regionMeyer ; if the electrostatic potential of this quantum well drops below
the Fermi level, it can become populated by the transfer of electrons from the
slab region: the threshold for the occurrence of this unwanted, artificial
behavior depends on both the width of the vacuum region and the strength of
the electric fieldMeyer . In our applications this upper limit turns out to be
around 1.0 V/Å, which is a value of the same order (or even larger) than the
maximum electric field achievable in actual experiments (for instance, in
electrowetting applicationsOstrowski ).

We computed the He-substrate interaction (with and without electric field) for
a selected set of non-equivalent sites in the primitive surface unit cell. In
particular, for He on GF, we chose the following points: $H$ (hollow site),
$TC$ (site on top of C atom), $TF$ (site on top of F atom), $B$ (bridge site
between C and F); for He on hBN we chose: $H$ (hollow site), $TB$ (site on top
of B atom), $TN$ (site on top of N atom), $B$ (bridge site between B and N);
finally, for He on G, we considered the points: $H$ (hollow site), $TC$ (site
on top of C atom), and $B$ (bridge site between 2 adjacent C atoms). The sites
for GF and hBN are shown in Fig. 1.

Figure 1: Sites chosen for the calculation of the He-GF (left panel) and He-
hBN (right panel) potentials. GF has a three-layer structure: The central
layer is occupied by the carbon atoms as in a graphene sheet except for a
small vertical offset between the two triangular sublattices, the upper plane
and the lower one contain the fluorine atoms. In the left panel C atoms are
shown as grey circles and F atoms as green circles. hBN is a single layer
structure isomorphous to graphene. In the right panel B atoms are shown as
yellow circles and N atoms as blue circles.

Besides the lowest-energy configurations for a given investigated adsorption
site, we have also computed the dependence upon the normal coordinate $z$ of
the He-substrate interaction potentials above those sites, shown in Fig. 2.
Our goal was to provide a reliable three-dimensional He-substrate potential
function to be used for numerical simulations based on Quantum Monte Carlo
methods, that will be discussed in the following.

Figure 2:  Calculated He-hBN and He-GF interaction potentials along the $z$ direction above selected sites. GF (from top to bottom): $TF$, $B$, $H$, $TC$; hBN (from top to bottom): $TN$, $B$, $TB$, $H$. The points show the ab-initio values, the lines are the results of 8-parameter curve fitting with the form $\sum_{i=1}^{2}a_{i}\,exp(-b_{i}z)-\sum_{i=1}^{4}c_{i}/z^{2i+2},$ whose parameters are given in Table 1.  Table 1:  Best-fit parameters for the $z$-dependence of the He-substrate potentials shown in Fig. 2.  |  | $a_{1}$ | $b_{1}$ | $a_{2}$ | $b_{2}$ | $c_{1}$ | $c_{2}$ | $c_{3}$ | $c_{4}$
---|---|---|---|---|---|---|---|---|---
GF: | $TF$ | $2.53645\times 10^{7}$ | $1.55209$ | $1.07222\times 10^{6}$ | $0.842514$ | $-4.49616\times 10^{6}$ | $7.37512\times 10^{8}$ | $-5.63025\times 10^{9}$ | $1.19823\times 10^{10}$
| $B$ | $-4.25772\times 10^{8}$ | $2.63907$ | $-1.0177\times 10^{7}$ | $1.61871$ | $-339981$ | $7.38605\times 10^{7}$ | $-3.58557\times 10^{9}$ | $1.1976\times 10^{10}$
| $H$ | $-2.73175\times 10^{6}$ | $1.91148$ | $-2.77692\times 10^{6}$ | $1.91101$ | $42978.4$ | $8.18279\times 10^{6}$ | $-3.9442\times 10^{8}$ | $1.54338\times 10^{9}$
| $TC$ | $526386$ | $0.975559$ | $20543.4$ | $0.389085$ | $4.68048\times 10^{6}$ | $8.32948\times 10^{6}$ | $-4.34654\times 10^{8}$ | $1.61137\times 10^{9}$
hBN: | $TN$ | $-244960$ | $1.94545$ | $-231170$ | $1.94594$ | $93994$ | $-1.97338\times 10^{6}$ | $1.6644\times 10^{6}$ | $2.86811\times 10^{6}$
| $B$ | $2.15693\times 10^{7}$ | $3.62009$ | $2.03143\times 10^{7}$ | $3.61956$ | $87334$ | $-2.0872\times 10^{6}$ | $2.18302\times 10^{7}$ | $-3.33307\times 10^{7}$
| $TB$ | $-208875$ | $1.90869$ | $-218649$ | $1.90819$ | $99505.8$ | $-2.25015\times 10^{6}$ | $4.40446\times 10^{6}$ | $-1.32199\times 10^{6}$
| $H$ | $1.75506\times 10^{7}$ | $3.56209$ | $1.7902\times 10^{7}$ | $3.56201$ | $89581.7$ | $-2.20804\times 10^{6}$ | $2.3581\times 10^{7}$ | $-3.80861\times 10^{7}$

We approximate such potential by using a truncated Fourier expansion over the
first three stars of the two-dimensional reciprocal lattice associated with a
triangular lattice with a two-atom basis (two C atoms in the case of G, one C
and one F atom in the case of GF, one B and one N atom in the case of h-BN).
The Fourier components can be obtained from the calculated $z$-dependence of
the various symmetry sites described above.nota

Table 2: Binding energy in the lowest-energy configuration (with and without external electric field), $E_{b}$, distance $d$ of He from the reference plane (defined by averaging over the $z$ coordinates of the C atoms for G and GF, and of the B and N atoms for hBN), maximum corrugation, $\Delta_{\rm{max}}$, minimum intersite energy-barrier, $\Delta_{\rm{min}}$ (see text for the definitions), for He-G, He-GF, and for He-hBN, using the rVV10 DFT functional. system | electric field (V/Å) |  $E_{b}$(K) | $d(\rm{\AA})$ |  $\Delta_{\rm max}$(K) |  $\Delta_{\rm min}$(K)
---|---|---|---|---|---
He-GF | 0.0 | -261 | 4.10 | 54 | 11
He-GF | 1.0 | -278 | 4.10 | 55 | 10
He-hBN | 0.0 | -214 | 2.96 | 36 | 16
He-hBN | 1.0 | -235 | 2.96 | 43 | 19
He-G | 0.0 | -298 | 2.96 | 50 | 47
He-G | 0.6 | -290 | 2.96 | 47 | 44

Our numerical results for the adsorption of He on the different substrates we
have considered are summarized in Table 2. In particular, we report the
distance $d$ of He from the substrate and the binding energy $E_{b}$ of the
the lowest-energy configuration. This configuration is found above the hollow
site $H$ for hBN and G, and above the $TC$ site for GF (see Fig. 2). For GF a
nearly isoenergetic configuration is found also above the hollow site
nava_2012 ; silvestrelli_2019 . Therefore on GF two essentially equivalent
kinds of adsorption sites are present and the overall number of adsorption
sites is double that on hBN and G. In Table 2 we also report two other
energetic parameters: the ”maximum corrugation”, $\Delta_{\rm{max}}$, defined
as the difference between the binding energy of He on top of C, F, and N atom
(which represents the less-favored configuration for G, GF, and hBN,
respectively) and the binding energy of the lowest-energy configuration, and
the ”minimum intersite energy barrier”, $\Delta_{\rm{min}}$, which is given by
the minimum energy barrier that the He atom must overcome to be displaced from
an optimal adsorption site to another, namely from $H$ to $H$ for hBN and G
and from $H$ to $TC$ for GF. This latter quantity has been evaluated by
monitoring the binding energy corresponding to a reaction path generated by
constraining the planar $x$, $y$ coordinates of the He atom and optimizing the
vertical $z$ coordinate only. In the case of hBN, $\Delta_{\rm{min}}$
($\Delta_{\rm{max}}$) corresponds to the difference between the binding energy
of He on top of the B(N) atom, $TB$($TN)$ site, and the binding energy of the
lowest-energy configuration, $H$.

For the case of no external electric field, the results for the corrugation of
the adsorption potential agree with those reported in previous
studies:nava_2012 ; silvestrelli_2019 ; nuova_nota  the most striking
difference between the case of He-G and of He-GF is that in He-G
$\Delta_{\rm{max}}$ and $\Delta_{\rm{min}}$ are comparable (the difference is
about 7%) while, on the contrary, in He-GF $\Delta_{\rm{min}}$ is smaller than
$\Delta_{\rm{max}}$ by a factor 5. This confirms that the adsorption potential
of He on GF is characterized by narrow ”canyons” between adsorption sites,
with a much larger anisotropy in the corrugation and a relatively low energy
barrier compared to G. A large difference between $\Delta_{\rm{min}}$ and
$\Delta_{\rm{max}}$ (about a factor 2) is found also for hBN so that the
corrugation is larger than in the case of G. Another difference between hBN
and G is that each adsorption site in hBN is surrounded by 3 saddle points and
not 6 as in the case of G and graphite. It should be noticed that the
adsorption energy of He on GF found with the semiempirical approach in Ref.
nava_2012,  appears to be strongly overestimated, we find that this energy is
about 10% smaller than that on G and not much larger as reported in Ref.
nava_2012, .

When an external, uniform electric field is applied (with the maximum strength
of $E$ allowed by our approximations, see discussion above), there is a change
of the adsorption potential, leading to moderate quantitative changes in the
quantities reported in Table 2, which, however, are no so significant to alter
the basic properties found at vanishing electric field.

##  III Quantum simulations

We use the rVV10 potentials of Sec. II to calculate unbiased equilibrium
thermal averages and ground-state properties of 4He adsorbed on GF, hBN and G
with quantum Monte Carlo simulations. For $N$ He atoms the Hamiltonian is

$H=-\frac{\hbar^{2}}{2m}\nabla^{2}+\sum_{i=1}^{N}+v_{1}({\bf
r}_{i})+\sum_{i<j}v_{2}(r_{ij}),$ (1)

where $v_{1}$ is the appropriate He-substrate one-body potential with ${\bf
r}_{i}$ the position of the $i$th 4He atom, and $v_{2}$ is the Aziz HFDHE2 He-
He pair potentialaziz  with $r_{ij}=|{\bf r}_{i}-{\bf r}_{j}|$. The substrate
is placed at $z=0$ with the 4He atoms in the positive $z$ semispace, and
periodic boundary conditions are applied in the $x$ and $y$ directions.

The finite temperature simulations are performed using the PIMC method
ceperley_1995  with the worm algorithm boninsegni_2006 . In this approach the
matrix element $\langle R|e^{-\beta H}|PR\rangle$, where $\beta$ is the
inverse temperature, $R=\\{{\bf r}_{1},\ldots,{\bf r}_{N}\\}$, and $PR$ is a
permutation thereof, is represented as a real-space convolution of $M$ high-
temperature density matrices $\langle R^{\prime}|e^{-\tau
H}|R^{\prime\prime}\rangle$, where $\tau=\beta/M$, which are in turn
approximated by suitable closed-form expressions. The results are unbiased in
the limit $\tau\to 0$. We use the primitive approximationceperley_1995  with
$\tau=0.002~{}{\rm K}^{-1}$. In the Appendix we show that the finite-$\tau$
systematic error is negligible for our purposes. We work in the grand
canonical ensemble, i.e. at fixed volume $V$, temperature $T$ and chemical
potential $\mu$.

The ground state simulations are performed with a projection Monte Carlo
method known as VPI ceperley_1995  or PIGS sarsa_2000 . This approach uses the
same real-space convolution representation of the operator $e^{-\beta H}$
mentioned above for PIMC to calculate expectation values on the quantum state
$\Psi_{\beta}=e^{-\beta H}\Psi$, where the trial function $\Psi$ is a closed-
form approximation to the ground state. We use a trial function of the form

$\Psi(R)=\exp\left[-\sum_{i}u_{1}(z_{i})-\sum_{i<j}u_{2}(r_{ij})\right],$ (2)

where the one-body He-substrate pseudopotential $u_{1}$ and the two-body He-He
pseudopotential $u_{2}$ are parametrized with 6 variational parameters each,
and optimized by energy minimization using the stochastic reconfiguration
methodrocca . For $\beta\to\infty$ the state $\Psi_{\beta}$ approaches the
exact ground state. Thus, in addition to the finite–$\tau$ bias, one has to
control the error incurred using a necessarily finite value of $\beta$. We use
$\beta=2~{}{\rm K}^{-1}$; examples of convergence in $\beta$ are shown in the
Appendix.

For all substrates we use either “small” or “large” simulation cells,
containing 180 and 720 adsorption sites, respectively. For GF we will also
present a few results for a larger cell with 2880 sites.

###  III.1 GF

The initial configuration for our PIMC simulations in the grand canonical
ensemble is the empty cell. After equilibration, the number of 4He atoms
fluctuates around a stationary average value $N$ which depends on the input
chemical potential $\mu$. Figure 3 shows a plot of $N(\mu)$ at various
temperatures.

Figure 3:  4He adsorbed on GF: number of atoms as a function of the chemical
potential calculated for various temperatures, as indicated, for the small
cell (main figure) and the large cell (inset). The dashed lines indicate the
numbers at filling 1/6.

For the small cell (main figure) the data for $T=2$ K follow a smooth curve,
corresponding to a fluid phase extending over the whole range of $\mu$ shown
in the plot. Upon lowering the temperature below $T=1$ K, a flat region at
$N(\mu)=30$ develops between $\mu=-185$ and -184 K. Here the 4He atoms form a
commensurate solid phase with areal density $\rho=0.0574$ Å-2, at 1/6 of the
adsorption sites. A fluid phase is found in a narrow region around $\mu=-186$
K, and for lower values of the chemical potential the density drops abruptly
to zero. At low temperature, these results depict the system as a modulated
self–bound superfluid which undergoes a first–order phase transition into the
1/6 commensurate crystal upon increasing the density. A similar behavior is
observed for the large cell (inset of Figure 3), apart from details in the
somewhat irregular increase of $N$ after the flat region, with a tendency to
form stripes in the small cell and domain walls in the large cell. We believe
that these differences are finite size effects caused through steric
constraints by the dimension and shape of the cells. We do not investigate
this aspect further because our main interest here is in the modulated liquid
and the commensurate crystal phases. For $T=0.25$ K we also simulate an even
larger cell with 11520 adsorption sites. Starting from the empty cell, for
$\mu=-186$ K we find the same liquid phase at areal density $\rho=0.050$ Å-2
as for the other cells, but for higher chemical potential it takes too long to
equilibrate a monocrystal at 1/6 covering. However, a starting crystal
configuration with 1920 atoms remains stable even after very long runs for
$\mu$ between -185 and -184 K.

Figure 4:  4He adsorbed on GF: two–dimensional structure factor calculated in
the small cell for various temperatures, as indicated in the body of the
figure, at chemical potential $\mu=-186$ K (left panels) and $\mu=-185$ K
(right panels). The highest peaks are 6.2 tall.

These phase assignments are supported by structural data. Figure 4 shows the
change in the structure factor $S({\bf k})$ across the liquid–solid
transition, both in temperature and chemical potential. The six highest peaks
in each panel correspond to the first shell of reciprocal lattice vectors,
$|{\bf k}|=1.618$Å-1, of a triangular crystal which occupies 1/6 of the
adsorption sites.

For $\mu=-186$ K, $S({\bf k})$ features the characteristic ridge nava_2012  of
a strongly modulated liquid. At all temperatures six peaks are present at
larger wavevectors with $|{\bf k}|=2.802$ Å-1; these peaks arise because the
fluid has a density modulation due to the substrate adsorption potential. We
note that the structure in the liquid phase is most pronounced at $T=1$ K; for
higher temperatures it is reduced by thermal fluctuations, for lower ones by
Bose exchanges ceperley_1995 . The liquid has a finite superfluid fraction
$\rho_{s}$ below a critical temperature located between $T=0.5$ and 1K. This
is shown in Figure 5 (left panel, blue symbols), where $\rho_{s}$ is seen to
exceed 50% in a wide range of system sizes for $T\leq 0.5$ K, while for $T=1$
K it has a much lower value which further drops toward zero as the number of
particles increases. Even at the lowest $T$, $\rho_{s}$ is less than 100% as
it is expected on general grounds for a non–uniform superfluid Leggett .

Figure 5:  4He adsorbed on GF: dependence on the system size of the superfluid
fraction and of the main peak of $S(k)$ for two areal densities, $\rho\sim
0.050$ Å-2 (blue) and $\rho\sim 0.057$ Å-2 (green). Left panel: $\rho_{s}$ as
a function of $T$ for various numbers $N$ of particles; for the higher density
the data are multiplied by a factor of ten. Right panel: main peak of $S(k)$
as a function of $N$ for $T=0.25$ and 1 K.

For $\mu=-185$ K, at low temperature the system is in the commensurate crystal
phase: the main peaks of $S({\bf k}$ soar to a value of 6.2 (Figure 4), which
further increases linearly with the number of particles, as shown for $T=0.25$
K by the green circles in the right panel of Figure 5. On the other hand for
$T=1$ the peak is still rather high (Figure 4), but it does not increase
further for larger systems (green asterisks in the right panel of Figure 5).
This indicates a melting temperature between 0.5 and 1 K. In addition to the
main peaks, one notices additional weak peaks at larger $k$, related to the
periodicity of the adsorption potential, as found in the superfluid state.

Figure 6:  One-body areal density profile $\rho(x,y)$ of 4He adsorbed on
various substrates at $T=0.25$ K. Left panels: from top to bottom, superfluid
at equilibrioum density on GF; commensurate solid on GF; commensurate solid on
hBN; commensurate solid on G in the presence of an external electric field of
1 V/Å. The contour level increment is 0.01, and the lowest level is 0.02 in
the top panel and 0.01 in the other ones. For GF, the adsorption sites
correspond to the maxima (black contour lines) in the top left panel; for hBN,
the (unoccupied) adsorption sites are marked with plus signs in the third left
panel. Right panels: $\rho(x,y)$ along the dotted line shown in the
corresponding left panel.

The commensurate crystal is not supersolid: the small but finite superfluid
fraction calculated in this phase with the small cell drops to zero as the
system size increases, even at the lowest temperatures considered (green
symbols in the left panel of Figure 5).

Further details on the areal density profile $\rho(x,y)$ in the superfluid and
the commensurate solid phases are given in Fig. 6. In the liquid phase (top
panels) the 4He atoms distribute rather uniformly along the potential valleys
connecting the adsorption sites. In the solid phase (second top panels)
$\rho(x,y)$ features strongly anisotropic peaks centered on an adsorption site
with shoulders on its its three nearest neighbours and faint tails on its six
second nearest neighbours. Such peaks are broader and lower than for
commensurate crystal phases of 4He on other substrates such as hBN or G. The
local density is spatially very anisotropic, the ratio between the largest and
the lowest local density is about 40.

Figure 7:  4He adsorbed on GF: energy per particle as a function of the areal
density $\rho$ at $T=0$ calculated in the small cell with projection time
$\beta=2$ K-1 (blue). The line is a polynomial fit, restricted to the data in
the liquid phase. For areal density $\rho$=0.0574 Å-2, the system is in the
commensurate solid phase, in agreement with the finite temperature
calculations. The inset shows the structure factor along the $y$ axis of
reciprocal space for $\rho=0.0536$, 0.0555 and 0.0574 Å-2, in order of
increasing peak heigth. For comparison, we include the energy and the
structure factor obtained with $\beta=0.25$ K-1 (green).

For $T=0$, in agreement with the low temperature results as expected, we find
a transition between a modulated liquid and a commensurate solid. The energy
per particle calculated for the small cell using the VPI method with a
projection time $\beta=2$ K-1 is shown in Figure 7 with blue points. For
$\rho<0.0547$ Å-2 we can fit the Monte Carlo data with a cubic polynomial,
yielding an equilibrium density $\rho_{0}=0.044\pm 0.001$ Å-2. For larger
densities, the energy suddenly drops slightly below the fitted curve. In
particular, for coverage 1/6, the structure factor is very similar (with minor
quantitative differences discussed in the Appendix) to that shown in the top
right panel of Figure 4, representative of the commensurate solid. The
estimated coexistence region extends between $0.0535\pm 0.005$ and 0.0574 Å-2.

We note that a large projection time is required to converge to the solid
solution starting from the trial function of Eq. (2), which represents a
liquid. As shown by the green symbols and curves in Figure 7, a value of
$\beta$ similar to that employed in Ref. nava_2012,  is perfectly adequate for
the liquid phase, but not sufficient for coverages 1/6 and higher.

###  III.2 hBN

hBN has one adsorption site per substrate unit cell, located at the hollow
point $H$ (see Figs. 1 and 2). In the $xy$ plane, this corresponds to a
minimum of the He-hBN potential, surrounded by saddle points at the $TB$ and
maxima at the $TN$ sites. Therefore the corrugation of the absorption
potential on hBN around an adsorption site differs from that on G and on
graphite: the angular periodicity is 120∘ for hBN and not 60∘ as for G.

Figure 8:  4He adsorbed on hBN: number of atoms as a function of the chemical
potential calculated for various temperatures, as indicated, for the small
cell. We found similar results for the large cell at $T=0.25$ K. The dashed
line indicates the number of particles at $\rho=0.0606$ Å-2, corresponding to
1/3 of the adsorption sites. Inset: contour plot of the density profile in the
normal fluid phase at $T=3$ K, $\rho=0.0606$ Å-2; the contour level increment
and the lowest level are both 0.01.

For sufficiently high temperature a monolayer of adsorbed 4He forms a normal
fluid phase. The density profile is directly shaped by the mentioned features
of the corrugation potential: there is a peak centered on every $H$ site and
elongated towards the three nearest $TB$ points, as shown in the inset of Fig.
8 for $T=3$ K.

Figure 9:  4He adsorbed on hBN: two–dimensional structure factor calculated in
the small cell for various temperatures, as indicated in the body of the
figure, at $\rho=0.0606$ Å-2 ($\mu=-146$ K, see Fig. 8). The highest peaks are
10.4 tall.

The dependence of the number of 4He atoms on the chemical potential, obtained
with PIMC simulations in the small cell for several temperatures, is shown in
Fig. 8. For $T<2$ K, we find that the average areal density stays constant
over a wide range of $\mu$, with a value $\rho=0.0606$ Å-2. This corresponds
to a commensurate solid, with 4He atoms occupying 1/3 of the adsorption sites.
The analysis of the structure factor across the melting temperature (which
turns out to be about 2 K, Fig. 9) and its dependence on the system size
(similar to that shown in Fig. 5 for GF) support this assignment.

Figure 10:  4He adsorbed on hBN: energy per particle as a function of the
areal density $\rho$ calculated in the small cell with VPI simulations. In
order to illustrate the convergence toward the unbiased results for projection
time $\beta\to\infty$, we report energies calculated with $\beta$=0.2, 0.3 and
2.0 K-1.

From the heigth of the peaks of $S({\bf k})$ and from the density profile,
plotted in the third–row panels of Fig. 6, we see that the commensurate solid
on hBN is significanlty more localized than on GF. Analogously to the normal
fluid phase at higher temperature, the density peaks in the solid are
elongated in the direction of the three closest saddle points of the
adsorption potential.

A more striking difference with the GF substrate is that, in PIMC simulations
carried out at low temperature ($T<0.5$ K), the system jumps directly from an
empty cell to a commensurate crystal upon increasing the chemical potential.
In particular, no (super)fluid phase should appear in the ground state of the
system.

This is confirmed by VPI calculations at $T=0$, which show that the
commensurate solid is by far the lowest–energy state of the system (Fig. 10)
and a fluid state is present at a lower density as a metastable state. Just
like in the case of GF, we note that this VPI result, matching the findings of
PIMC calculations at low temperature, requires larger projection times than
used in Ref. silvestrelli_2019, .

###  III.3 Effect of an external electric field

We have performed PIMC simulations of 4He at $T=0.25$ K in the presence of an
external electric field $E=1$ V/Å on hBN and $E=0.6$ V/Å on G, the two
substrates with the larger dependence on $E$ of the adsorption potential
corrugation (see Table 1).

Figure 11:  4He adsorbed on hBN: variation of the density profile
$\Delta\rho(x,y)$ induced by the external electric field. Left panel,
greyscale map ranging from -0.08 (black) to 0.08 Å-2 (white); each of the
pictured objects with three negative and three positive lobes is centered on a
density peak of the commensurate solid (see Fig. 6, third left panel). Right
panel, variation $\Delta\rho(x,y)$ along the dotted line shown in the left
panel; the two structures displayed are not identical because of statistical
noise (note the small scale of the plot).

The external field tends to delocalize the adsorbed atoms: for hBN, the range
of chemical potential where 4He is stuck in the commensurate solid phase
shrinks by $\sim 3$ K; furthermore the peaks in the density profile get
further elongated towards the nearest saddle points of the He–substrate
potential. This can be seen by comparing the variation $\Delta\rho(x,y)$ of
the density profile, shown in Fig. 11, with $\rho(x,y)$ itself (Fig. 6, third
left panel): the directions of the positive (white) lobes of $\Delta\rho$
match those of the triangular stretching of the peaks.

However the effect is small (in fact, not visible on the scale of Fig. 6), and
for both hBN and G we find that the commensurate solid remains stable, with no
superfluid phases at lower densities.

The density profile of the commensurate solid on G with external field is
shown in Fig. 6; for a study of 4He adsorbed on G without external field see
Ref. kwon, ).

##  IV Summary and conclusion

In this paper we have revisited the adsorption at low coverage of 4He on
fluorographene and on hexagonal boron nitride in the sub monolayer regime and
we have also studied such adsorption in presence of an electric field for GF
and hBN and also for graphene. The motivation was the search of a superfluid
state of monolayer 4He adsorbed on very regular substrates, as can be
experimentally obtained with the above-mentioned substances, because such
superfluid state should be characterized by a strong spatial anisotropy, a
regime not yet explored with 4He. For instance, in this regime rotons should
be anisotropic with an energy depending on the direction of the wave vector
Nav13  and a vortex excitation should not be translationally invariant, but it
should have preferential sites for the location of the vortex core Gal2020 .

We have developed new adsorption potentials for such substrates with ab-initio
method and we have studied the properties of 4He at finite temperature with
PIMC and at $T=0$ K with VPI. In the case of GF we confirm previous results
nava_2012 ; silvestrelli_2019  that sub monolayer 4He in its ground state and
at low $T$ is a low-coverage self-bound superfluid. A commensurate state is
present at higher coverage in which one sixth of the adsorption sites are
occupied. This commensurate state was not detected in the previous studies. In
the present work we have not studied the system at higher coverage. The
superfluid fraction at the lowest $T$ is about 55% and not 100%, a depletion
expected even at $T=0$ K for a non-uniform superfluid Leggett . We estimate a
transition temperature to the normal state between 0.5 and 1.0 K. We have
characterized the structural properties in direct and in reciprocal space
confirming the extremely large spatial anisotropy of the superfluid. The
superfluid fills the bonds of a honeycomb lattice of the adsorption sites with
a ratio of about 6 between the largest and the smallest areal density: it is
like the superfluid were moving in a multi-connected space. The commensurate
state at coverage 1/6 of the adsorption sites is a triangular lattice similar
to the commensurate state of 4He on graphite at coverage 1/3 but special
features characterize the commensurate state on GF. In fact, even at the
lowest $T$ a 4He atom is not constrained to remain at a single adsorption site
but there is a sizeable probability of occupation of the neighboring sites
giving rise to a three-lobed density distribution. As far as we know this is
the first instance of such behavior among the ordered adsorbates on regular
substrates and it will be of interest to determine if such feature is
modifying the order-disorder phase transition compared to that of 4He on
graphite reatto_2013 . Other features of this ordered state of 4He on GF are a
mean square deviation from the equilibrium site more than double that of 4He
on graphite and a particularly low transition temperature to a disordered
state between 0.5 and 1 K, less than one third of the order-disorder
transition of 4He on graphite. This commensurate state is not supersolid.

We find that the lowest energy state of 4He on hBN at $T=0$ K and at low $T$
is a commensurate triangular solid in which one third of the adsorption sites
are occupied, a state isomorphous to that of 4He on graphite and graphene. The
temperature of the order-disorder transition is about 2 K. The present result
modifies the conclusion of a previous study silvestrelli_2019  in which it was
found that the ground state was a superfluid. We have verified that this
discrepancy is due to the use of a too short projection time in the earlier
$T=0$ K study. At coverage below this commensurate state a fluid state is
present as a metastable one. We have shown that an electric field up to 1.0
V/Å tends to delocalize the adsorbed atoms, but the effect is small and the
electric field does not modify the phase behavior of the system. We find that
also in the case of G the ground state of sub monolayer 4He remains an ordered
state commensurate with the adsorption sites at coverage 1/3.

Figure 12:  The $x,y$ classically allowed regions for a 4He atom on GF, hBN
and G at a distance $d$ from the substrate corresponding to the minimum value
$V_{\rm min}$ of the adsorption potential. The color maps represent the
adsorption potential $V(x,y)-V_{\rm min}$ from zero to the kinetic energy of a
single 4He atom on the substrate calculated with PIMC at $T=0.25$ K (i.e.
27.0, 22.6 and 28.4 K for GF, hBN and G respectively). The color scale ranges
from 0 (black) to 28.4 K (yellow). The plots are centered on the less favored
site (F atom for GF, N atom for hBN and C atom for G). In the superfluid phase
of 4He on GF (top panels of Fig. 6), the adsorbate atoms cover rather
uniformly the interconnected regions of low potential energy.

In conclusion, on present theoretical evidence the only regular substrate that
supports a superfluid ground state is GF. Possibly this is also true for
graphane nava_2012 , a compound isomorphous to GF with the fluorine atoms
replaced by the hydrogen ones, but experimentally it is difficult to produce
regular substrate with full stoichiometry and we have not investigated this
system with ab-initio methods. The different nature of the lowest energy state
of 4He on GF compared to that on hBN, G and graphite is mainly due to the
number of adsorption sites per surface unit cell, twice as many for GF as for
the other substrates. As seen in Fig. 12 this entails a larger, interconnected
region of favorable, low potential energy available to the Helium atoms. The
size and shape of such region stabilize the superfluid phase and induce the
three-lobed distortion of the peaks in the commensurate solid.

Other commensurate states might be present at larger coverage between the
density range of the present study and the promotion coverage to the second
layer, this is an interesting topic for future studies. Another interesting
development is the study of the adsorption of the fermionic 3He on GF and on
hBN. Due to the smaller mass of 3He it is very unlikely that the ground state
of 3He on GF is a commensurate ordered state so that the low coverage state
should be a liquid or an unbounded gas. Some evidence for a liquid state was
found in the earlier study nava_2012  with the semiempirical adsorption
potential so further study of the fermionic system with the new ab-initio
adsorption potential is of great interest because, in any case, 3He on GF
opens the possibility of studying a strongly interacting Fermi system with
large spatial inhomogeneity. We have shown that the most stable state of 4He
on hBN is ordered and a fluid state is only a metastable one. The smaller mass
of 3He might alter the balance between these two states. In fact, we have
performed PIMC computations for bosonic m=3 and indeed we find that with this
mass the fluid state if the stable one. Of course, the Fermi statistics might
alter this result, and this is left for future study. On this issue it will be
important to perform new experiments on 3He on hBN to verify the earlier
measurements crane_2000  that gave some evidence for a localized commensurate
state.

## Acknowledgments

SM acknowledges support from the European Centre of Excellence in Exascale
Computing TREX, funded by the European Union’s Horizon 2020 - Research and
Innovation program - under grant agreement no. 952165.

## Appendix

The PIMC simulations ceperley_1995 ; boninsegni_2006  give unbiased thermal
averages for $\tau\to 0$. We use the primitive action with $\tau=0.002$ K-1
which entails a negligible error in the quantities of main interest here,
namely the two–dimensional structure factor $S({\bf k})$ and the superfluid
fraction $\rho_{\rm s}$. This is shown in Figure 13 for $N=30$ 4He atoms
adsorbed on GF in the 1/6 commensurate solid phase (here the non–zero value of
the superfluid fraction is a finite–size effect).

Figure 13:  The structure factor along the $y$ axis of reciprocal space for
$N=30$ 4He atoms adsorbed on GF in the 1/6 commensurate solid phase,
calculated with PIMC at $T=0.25$ K. Data for different values of $\tau$ are
slightly offset in $k$ for clarity. The inset shows the superfluid fraction as
a function of $\tau$.

The VPI simulations ceperley_1995 ; sarsa_2000  give unbiased ground state
results for time step $\tau\to 0$ and projection time $\beta\to\infty$. We use
the primitive action with the same value $\tau=0.002$ K-1 as in PIMC, and
$\beta=2$ K-1. The convergence in $\beta$ is shown in Figure 14 for the
structure factor and the energy.

Figure 14:  The structure factor along the $y$ axis of reciprocal space for
$N=30$ 4He atoms adsorbed on GF at the density of the 1/6 commensurate phase,
calculated with VPI. Data for different values of $\beta$ are slightly offset
in $k$ for clarity. The PIMC result for $T=0.25$ K is also shown for
comparison. The inset shows the energy per particle as a function of $\beta$,
relative to its value at $\beta=2$ K-1.

In the range of $\beta$ between 0.25 and 2 the energy changes merely by $\sim
50$ mK, but the structure factor develops large peaks at the first reciprocal
lattice vector of the commmensurate solid. The convergence in $\beta$ of the
main peaks is not complete, because their heigth should be closer to that
obtained with PIMC at the lowest temperatures (also shown in Figure 14 for
comparison). Nevertheless, a projection time of $\beta=2$ K is clearly
sufficient to turn the liquid state represented by the trial function of Eq.
(2) into a solid.

## References

  * (1) M. Greiner, O. Mandel, T. Esslinger, T.W. Hansch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature 415, 39 (2002).
  * (2) Sec. 6.1 of L. W. Brunch, M. W. Cole, and E. Zaremba, Physical Adsorption: Forces and Phenomena (Dover Publishing, Mineola, NY, 2007).
  * (3) M.E. Pierce, and E. Manousakis, Phys. Rev. Lett. 83, 5314 (1999).
  * (4) P. Corboz, M. Boninsegni, L. Pollet, and M. Troyer, Phys. Rev. B 78, 245414 (2008).
  * (5) P. A. Crowell, and J. D. Reppy, Phys. Rev. B 53, 2701 (1996).
  * (6) J. Nyeki et al., Nat. Phys. 13, 455 (2017)
  * (7) M. Nava, D.E. Galli, M.W. Cole, L. Reatto, Phys. Rev. B 86, 174509 (2012).
  * (8) L. Reatto, D.E. Galli, M. Nava, M.W. Cole, J. Phys. Condens. Matter 25, 443001 (2013)
  * (9) P.L. Silvestrelli, M. Nava, F. Ancilotto, and L. Reatto, J. Low Temp. Phys. 196, 42 (2019).
  * (10) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
  * (11) A. Sarsa, K.E. Schmidt, W.R. Magro, J. Chem. Phys. 113, 1366 (2000).
  * (12) M. Rossi, M. Nava, L. Reatto, D.E. Galli, J. Chem. Phys. 131, 154108 (2009).
  * (13) L. Reatto, M. Nava, D.E. Galli, C. Billman, J.O. Sofo, and M.W. Cole, J. Phys.: Conf. Ser. 400, 012010 (2012).
  * (14) M. Nava, D.E. Galli, M.W. Cole, and L. Reatto, J. Low Temp. Phys. 171, 699 (2013).
  * (15) M.J. Stott, and E. Zaremba, Phys. Rev. B 22, 1564 (1980); G. Vidali, M.W. Cole, and C. Schwartz, Surf. Sci. 87, L273 (1979); K. T. Tang, and J.P. Toennies, J. Chem. Phys. 80, 3726 (1984).
  * (16) P. Giannozzi et al., J.Phys.: Condens. Matter 21, 395502 (2009).
  * (17) R. Sabatini, T. Gorni, S. de Gironcoli, Phys. Rev. B 87, 041108(R) (2013).
  * (18) O. A. Vydrov, T. Van Voorhis, J. Chem. Phys. 133, 244103 (2010).
  * (19) P. L. Silvestrelli, A. Ambrosetti, Phys. Rev. B 91, 195405 (2015).
  * (20) P. L. Silvestrelli, A. Ambrosetti, J. Low. Temp. Phys. 185, 183 (2016).
  * (21) J. Neugebauer, and M. Scheffler, Surf. Sci. 287/288, 572 (1993).
  * (22) L. Bengtsson, Phys. Rev. B 59, 12301 (1999).
  * (23) B. Meyer, and D. Vanderbilt, Phys. Rev. B 63, 205426 (2001).
  * (24) J. H. J. Ostrowski, and J. D. Eaves, J. Phys. Chem. B 118, 530 (2014).
  * (25) Three dimensional tabulations of the He-substrate potentials are available upon request from one of the authors (F.A.).
  * (26) In particular, minor numerical differences in the adsorption potential with respect to the data reported in Ref. silvestrelli_2019,  using the same functional are due to slightly better convergence in energy and in the optimization of the equilibrium structure, and in the use of a wider vacuum region separating the periodically repeated slabs.
  * (27) R. A. Aziz, V. P. S. Nain, J. S. Carley, W. L. Taylor, and G. T. Conville, J. Chem. Phys. 70, 4330 (1979).
  * (28) M. Boninsegni, N. V. Prokof’ev, B. V. Svistunov, Phys. Rev. E 74, 036701 (2006).
  * (29) S. Sorella, M. Casula, and D. Rocca, J. Chem. Phys. 127, 014105 (2007).
  * (30) A. J. Leggett, J. Stat. Phys. 93, 927 (1998).
  * (31) Y. K. Kwon, and D. M. Ceperley, Phys. Rev. B 85, 224501 (2012).
  * (32) A. Gallemi, S. M. Roccuzzo, S. Stringari, and A. Recati, Phys. Rev. A 102, 023322 (2020); F. Ancilotto, M. Barranco, M. Pi, and L. Reatto, Phys. Rev. A 103, 033314 (2021).
  * (33) T. P. Crane, and B. P. Cowan, Phys.Rev.B62,11359(2000); Y. Tang, and N. S. Sullivan, J. Phys. Conf. Ser. 568, 012018 (2014).