Table Of Content

Mott transition in bosonic systems: Insights from the variational approach Manuela Capello,1 Federico Becca,2,3 Michele Fabrizio,2,3,4 and Sandro Sorella,2,3 1 Laboratoire de Physique Théorique, Université Paul Sabatier, CNRS, 31400 Toulouse, France 2 International School for Advanced Studies (SISSA), I-34014 Trieste, Italy 3 CNR-INFM-Democritos National Simulation Centre, Trieste, Italy. 4 International Centre for Theoretical Physics (ICTP), P.O. Box 586, I-34014 Trieste, Italy (Dated: February 2, 2008) WestudytheMotttransitionoccurringforbosonicHubbardmodelsinone,two,andthreespatial dimensions,bymeansofavariationalwavefunctionbenchmarkedbyGreen’sfunctionMonteCarlo calculations. We show that a very accurate variational wave function, constructed by applying a 8 long-range Jastrow factor to the non-interacting boson ground state, can describe the superfluid- 0 insulator transition in any dimensionality. Moreover, by mapping the quantum averages over such 0 a wave function into the the partition function of a classical model, important insights into the 2 insulating phase are uncovered. Finally, the evidence in favor of anomalous scenarios for the Mott n transition in twodimensionsarereported wheneveradditional long-range repulsiveinteractions are a added to theHamiltonian. J PACSnumbers: 74.20.Mn,71.10.Fd,71.10.Pm,71.27.+a 7 ] l I. INTRODUCTION scenario has been confirmed mostly in one and two di- e mensions by using different numerical techniques, such - r The recent advances achieved on cold atoms trapped as quantum Monte Carlo and density-matrix renormal- t s inopticallatticeshavegeneratedanincreasinginterestin ization group.5,6,7,8,9 In particular, it has been verified . t thecondensedmattercommunity,sincetheyallowexper- that in one dimension at ρ = 1 there is a Kosterlitz- a m imental realizations of simple lattice models.1,2 A great Thouless transition and the estimation of the critical advantage of these systems is the possibility to have a valueoftheon-siteinteractionrangesbetweenUc/t 1.7 d- directcontroloftheparameters,suchasthewidthofthe and Uc/t 2.3.5,7 Instead, a second-order phase tr∼ansi- ∼ n bands and the strength of the interactions, that can be tion is claimed to occur for Uc/t 8.5 in two dimen- o manipulated by varying the depth of the optical poten- sions,9 and for Uc/t 14.7 in thre∼e dimensions.10. Ac- ∼ c tial.3 Therefore, cold atoms on optical lattices give the curate estimations for the critical values can be also ob- [ uniqueopportunitytomakeaclosecontactwiththeoret- tained from a strong-coupling expansion.11,12 1 ical models and to examine the origin of the fundamen- Besides the numerically exact techniques, important v tal physical phenomena that occur in crystalline materi- insights into the various phases can be obtained by con- 2 als. Inparticular,oneofthemostspectacularexampleis sidering simplified variational wave functions. The sim- 5 9 given by the superfluid-insulator transition in a system plestexampleisgivenbythecelebratedGutzwillerstate, 0 of interacting bosons: The so-called Mott transition.2 wheretheon-sitecorrelationtermallowsforthesuppres- . In this paper, we consider the superfluid-insulator sion of the energetically expensive charge fluctuations. 1 0 transition in a bosonic system with different kinds of in- Contrarytothefermioniccase,whenconsideringbosons, 8 teractions. The simplest model that contains the basic itispossibleto dealwiththis wavefunctionwithoutany 0 ingredients of strong correlations is the Hubbard model further approximation.13,14 Then, the Mott transition is : t U obtainedwithareasonableestimateofthecorresponding iv H=−2 b†ibj +h.c.+ 2 ni(ni−1), (1) criticalvalueUc/t,namelyUc/t=d(√nc+√nc+1)2 for X hXi,ji Xi integer fillings ρ = nc. The main drawback of this ap- ar where h...i indicates nearest-neighbor sites, b†i (bi) cre- pthreoatcrhanissitthioant,issirmeailcahrelydtwoitwhhaatvhaanpispheinnsgwkiinthetfiecremnieorngsy, ates (destroys)a bosononthe site i, andn =b†b is the i i i and the insulating state completely lacks chargefluctua- local density operator. Here, we consider N particles on tions,namelyallparticlesarelocalized,n ineachlattice a one-dimensional (1D) chain, a two-dimensional (2D) c site. Of course, this leads to a wrong description of the square lattice, and a three-dimensional (3D) cubic lat- insulator, whenever the local interaction is finite. tice with L sites and periodic-boundary conditions. At zero temperature and integer densities ρ = N/L there Followingtheideasofpreviousworksonfermionicsys- is a superfluid-insulator transition when the ratio be- tems,15,16,17 in a recent paper18 we have shown that, in tween the kinetic energy and the on-site interaction is order to correct this outcome and obtain a more accu- varied. Otherwise, for non-integer fillings, the ground rate description of the ground state, it is necessary to state is always superfluid. In a seminal paper,4 by using include a long-range Jastrow factor, whose singular be- a field-theoretical approach, Fisher and coworkers pro- havioratsmall momenta wasshownto be able to turn a posed that the transition of the d-dimensional clean sys- non-interacting bosonic state into an insulator that still tem belongs to the XY universality class in d+1. This contains density fluctuations. In this paper, we present 2 a more detailed study of the properties of that Jastrow condensed into the q =0 state wave function, as well as of its accuracy in comparison withGreen’sfunctionMonteCarlo(GFMC)calculations, which, because of the absence of the sign problem, pro- 1 Ψ =exp v (n 1)(n 1) Φ , (5) vide numerically exact results.19,20 | Ji −2 i,j i− j − | 0i The same Jastrow wave function turns out to be also  Xij  veryeffectivetodescribeHamiltoniansthatcontainlong-   range interactions. This case has not been considered where Φ = ( b†)N 0 is the non-interacting boson | 0i i i | i much in the literature, and we will show that different ground state with N particles, (n 1) is the variation i P − scenariosforthesuperfluid-insulatortransitioncouldoc- of the on-site density with respect to the average value cur. As a matter of fact, long-range interactions have ρ = 1, and v are translationally invariant parameters i,j beenstudiedmostlyinthecontinuum,whereatransition that can be optimized to minimize the variational en- betweenachargedbosonicfluidandaWignercrystalhas ergy.24 In order to get some physical insight from the beenfoundbyvaryingthedensity.21,22,23Inparticular,in variational state, it is more instructive to consider the the presence of a logarithmic repulsion, the 2D Bose liq- Fourier transformof the Jastrow parameters v . Indeed, q uidwasfoundto havenocondensatefraction,due tothe there is a tight connection between the small-q behavior predominance of long-wavelengthplasmon excitations.21 of v and the nature of the ground state. In particular, q Inthe lastpartofthis paper,we generalizethe Hubbard as we are going to discuss, v 1/q implies the exis- q ∼ | | model of Eq. (1) to tenceofsoundmodesinanydimensions,asexpectedina superfluid. On the contrary, to recoveran insulating be- t V HLR =−2 b†ibj +h.c.+ 2 Ω(Ri,Rj)ninj, (2) havior a much more singular vq for q → 0 is required.18 hXi,ji Xi,j The physical reason is that, in order to describe a re- alistic Mott insulating wave function that does include where Ω(R ,R ) is a long-range potential that only de- i j charge fluctuations, it is necessary to spatially correlate pends upon the relative distance R R between two | i− j| the latters. This is accomplished by a sufficiently sin- particles and V represents its strength. In particular, gular v that favors configurations where opposite-sign we will consider two possibilities for the long-range po- q fluctuations,(n 1)(n 1)<0,areclosetoeachother, tential. The first one is obtained by taking the Coulomb i− j− while equal-signones,(n 1)(n 1)>0,arefarapart. interactionbetween(charged)bosonsmovingina2Dlat- i− j− ticeembeddedinathree-dimensionalenvironment,which It should be mentioned that previous studies of leads to the following potential in reciprocal space fermionicsystems25 andmorerecentonesonthe bosonic Hubbard model,26 stressed the importance of a many- π Ω(qx,qy)= , (3) body term containing holon-doublon interactions for (cosq +cosq 3)2 1 x y − − nearest-neighbor sites: with a small-q behavpior Ω(q) 1/q . The second case ∼ | | consistsindirectlyconsideringthesolutionofthePoisson equation in 2D: Ψ =exp g n2+g ξ Φ , (6) | g,MBi − i MB i!| 0i 1 i i X X Ω(q ,q )= , (4) x y 2 (cosq +cosq ) x y − where g and g are variational parameters and the MB which for small momenta behaves like Ω(q) 1/q2, many-body operator is defined by ∼ leading to a logarithmic interaction in real space, i.e., Ω(r) log(r). In both cases, a uniform background is ∼− ξ =h (1 d )+d (1 h ), (7) consideredin orderto cancelthe divergentq =0 compo- i i i+δ i i+δ − − nent of the potential. δ δ Y Y The paper is organized as follow: In section II, we describe the variational wave function, in sections III, IV whereh =1(d =1)ifthe siteiis empty(doubly occu- i i andV, weshowthe results forthe 1D,2Dand3Dshort- pied) and 0 otherwise, δ = x, y; therefore, ξ counts i ± ± range models. Then, in section VI, we consider the thenumberofisolatedholons(emptysites)anddoublons 2D case with long-range interactions and finally, in sec- (doubly occupied sites). Even though this operator has tion VII, we draw our conclusions. been originally introduced for fermionic systems, where the maximum occupancy at each site is given by two electrons,it is usefulalsoforbosons,since inthe limitof II. THE VARIATIONAL WAVE FUNCTION large U/t the number of sites with an occupation larger than two is negligible. However, within this framework, A. The Jastrow wave function the evidence for a true Mott transition is rather contro- versial, and it is not clear if an insulating phase can be The variational wave function is defined by applying stabilized in the thermodynamic limit.26 a density Jastrow factor to a state with all the bosons Combining together the variational Eqs. (5) and (6), 3 we obtain the variational ansatz through a potential determined by v . In the presence q of the short-range many-body term, V (x) is no longer cl 1 Gaussian. However, when the density fluctuations are Ψ = exp v (n 1)(n 1) | J,MBi −2 i,j i− j − suppressed at large U/t, the quadratic term gives the  Xi,j most relevant contribution, hence the mapping onto a classical model of interacting oppositely charged parti-  + gMB i ξi)|Φ0i, (8) bcleeinsgtilrlephloalcdesdwbiythaVsclli(gxh)tl∼y diffqe6=re0nvtqefefffneqc(txiv)en−pqo(txe)n,tivaql X veff. In spite of the differencPes, it is plausible that the containing both the long-range Jastrow factor and a q small-q behavior of veff must follow the same singular short-range many-body term. As it will be shown in the q behavior of the Jastrow potential v . next sections, the presence of the latter term is impor- q Let us now discuss in more details the connection be- tantin2Dand3D,mainlytoincreasetheaccuracyinthe tweenthe formofv andthe low-energyexcitationspec- strong-coupling regime. Instead, in 1D, the many-body q trumofthe system. Bymeansofthe f-sumruleone can term does not improve the accuracy of the long-range show that18,27 Jastrow state and there is no appreciable difference between the wave function (5) and (8). In any case, we emphasize that, accordingto our calculations,the short- k D q dωωN(q,ω) range term alone cannot lead to an insulating behavior, Eq = 2 h i sin2 i = Z , (12) − DN 2 and the main ingredient to drive a superfluid-insulator q Xi=1 (cid:16) (cid:17) dωN(q,ω) transitionis the long-rangeJastrowfactor,parametrized Z by vq. where N(q,ω) is the dynamical structure factor, Nq = dωN(q,ω) the static one, and k the hopping energy h i per site. E can be interpreted as the averageexcitation q B. Mapping onto a classical model eRnergy of density fluctuations at momentum q. E 0 q → for q 0 is a sufficient but not necessary condition for → The variationalcalculationwith the wave function (8) the existence of gapless excitations. In particular, deep can be shown to correspond to a classical problem at fi- inside the superfluid phase, E must coincide with the q nite temperature. This correspondence provides many energy dispersion of the Bogoliubov sound mode, which insights into the properties of Ψ . To prove the carriesmostofthespectralweight. Althoughwecannot J,MB | i mapping,letusdenoteabosonicconfigurationbythepo- access directly dynamical properties by our variational sitions x of the particles and a generic quantum state wavefunction, stillwecanprovideavariationalestimate { } by Ψ . For all the operators θ diagonal in space coordi- of E , in the same spirit of the Feynmann’s construction | i q nates, the quantum average in liquid Helium.28 This amounts to use the variational valuesof k andofthestaticstructurefactorN ,defined q Ψθ Ψ h i θ = h | | i (9) as h i ΨΨ h | i Ψn n Ψ q −q N = h | | i, (13) can be written in terms of the classical distribution q ΨΨ Ψ(x)2 = xΨ 2/ x′ Ψ 2, as h | i | | |h | i| x′|h | i| with n the Fourier transform of the local density. Note q θ P= xθ x Ψ(x)2. (10) that the uncorrelated h i h | | i| | Xx N0 = hΦ0|n−qnq|Φ0i (14) Since Ψ(x)2 isapositivequantity,thereisaprecisecor- q Φ0 Φ0 | | h | i respondence between the wave function and an effective is constant for any q = 0 at ρ = 1. Whenever a weak- classical potential Vcl(x): couplingapproachint6heJastrowpotentialv ispossible, q the following relation holds:29 Ψ(x)2 =e−Vcl(x). (11) | | N0 1 tThheecehxopilciecitoffotrhme oJfastthreowpoftaecnttoiar,l Vwchl(ixle) |dΦe0pienddosesupnoont Nq =Xx nq(x)n−q(x)e−Vcl(x) ∼ 1+2Nq q0vq ∼ vq, contribute to it, since xΦ does not depend upon the (15) 0 h | i configuration x . In particular, whenever there is only the last equality following from the singular behavior of | i the two-body potential (i.e., g = 0) V (x) is Gaus- v that is expected both in the superfluid and in the MB cl q sian, i.e., V (x)= v n (x)n (x), being n (x) the Mott insulating phases. Eq. (15) shows that v 1/q cl q6=0 q q −q q q ∼ | | Fouriertransformofthelocaldensityoftheconfiguration allows to recover the correct behavior of N q and q P ∼ | | x . Inthiscase,thevariationalproblembecomesequiva- E q in the superfluid regime, which is what we in- q | i ∼ | | lenttosolveaclassicalmodelofoppositelychargedparti- deedfind, seefollowingsections. Inthe insulatingphase, cles (the holons and the doublons) mutually interacting should Eq. (15) be valid, we would expect v β/q2 q ∼ 4 to get the expected behavior N q2 and E finite for q ∼ q 0.25 small q. This would correspond through Eq. (11) to a Jastrow classical Coulomb gas model (CGM) with effective tem- Gutzwiller+MB peratureT =π/β.30In1D,foranyfinitetemperature, 0.2 Jastrow+MB eff the CGM is alwaysin a confined phase where oppositely charged particles are tightly bound in pairs, and with 0 0.15 exponential decaying correlation functions.31 This sug- E / E gests that the 1D CGM may indeed provide throughthe ∆ 0.1 mapping(11)agooddescriptionofa1DMottinsulating wave function. 0.05 In 2D the CGM undergoes a Berezinskii-Kosterlitz- Thouless phase transition at TCGM 1/4, between a confinedphase(stableatlowtemcperat∼ure)andaplasma 0 0 1 2 3 4 phase (stable at high temperature). Similarly to the 1D case, one would argue that the confined phase should U/t correspond to the 2D Mott insulator. However, the 2D confined phase of the CGM displays power-law decaying FIG. 1: Accuracy of different variational wave functions as a correlations30thatwouldcorrespondtopower-lawdecay- functionofU/tfor60sitesand60bosons. ∆E =E0−EVMC, ingequal-time correlationsofthe quantumgroundstate. where EVMC is thevariational energy and E0 is the ground- stateone,obtainedbyGFMC.ThestateofEq.(5)isdenoted This is not compatible with a genuine Mott insulator, by “Jastrow”, the one of Eq. (6) by “Gutzwiller+MB”, and whichalreadysuggeststhatthe insulating wavefunction theone of Eq. (8) by“Jastrow+MB”. must be characterized by a Jastrow potential v more q singularthan 1/q2 as q 0,as indeed we find. Inturns, → this impliesthat(15)mustnotbe valid,since,inspiteof vq q2 for q 0, we still expect Nq q2. 6 U/t=2 T×he b→rea∞kdown →of Eq. (15) becomes ev∼en more pro- U/t=2.4 U/t=2.5 nounced in 3D, where a potential vq 1/q2 cannot de- 5 U/t=3 ∼ scribe at all an insulator since it is not sufficiently sin- U/t=4 2 4 gular to empty the condensate fraction.32 In fact, we q| U/t=6 findthattheoptimizedvariationalwavefunctionshowsa x | 3 q more divergingv 1/q 3 in the 3D insulator,18 though v q ∼ | | 2 N q2. The properties connected to the insulating q ∼ phase can be again uncovered within a classical 3D gas 1 1D with a 1/q 3 potential, recently considered in Ref. 33. | | 0 Indeed, analogouslyto whathappens in 2D,alsoin3Da 0 π/4 π/2 3π/4 π system of charges interacting with a logarithmic poten- |q| tial in realspace admits a high-temperature fluid regime separatedbya low-temperaturedielectric phase. Within FIG.2: Jastrowparametersv multipliedbyq2 asafunction q this mapping, the 3D insulating state is found to cor- of |q| for 60 sites and 60 bosons. respond to the low-temperature phase of this classical model. Inthefollowing,wewillpresentourresultsobtainedby difference betweenthe wavefunctions (5) and(8), forall considering the quantum variational wave function and the values of the on-site interaction. By contrast, the wewillusetheclassicalmappingtogaininsightsintothe Gutzwiller state, also when supplied by the many-body ground-stateproperties. Moreover,inorderto verifythe term, is much less accurate by increasing U/t. accuracy of the variational calculations, we will perform Therefore, in the following, we will consider the state the numerically exact GFMC that allows one to obtain withlong-rangeJastrowfactor(5)alone,sincethemany- ground-state properties. bodytermmakesthealgorithmmuchsloweranddoesnot improvethe quality ofthe variationalstate. In Fig 2, we report the minimized Jastrow parameters v multiplied q III. THE 1D HUBBARD MODEL by q2 for different U/t: There is a clear difference in the small-q behavior for U/t.2.4, where v α/q and for q ∼ | | Here, we consider a chain of L sites with periodic U/t&2.5,wherev β/q2. Atthevariationallevel,the q ∼ boundary conditions and N = L bosons. First of all, in change of the singular behavior of the Jastrow parame- Fig. 1 we compare the variational accuracy of the wave ters for U/t 2.45 marks the superfluid-insulator tran- ∼ functions(5),(6)and(8)fordifferentvaluesofU/t. Once sition. Indeed, as discussed in the previous paragraph, along-rangeJastrowfactorisconsidered,themany-body v α/q implies a gapless system, whereas v β/q2 q q ∼ | | ∼ term of Eq. (7) is irrelevant and there is no appreciable indicates a finite gap in the excitation spectrum. Let 5 VMC GFMC 50 0.4 U/t=2.4 U/t=2.5 0.35 40 0.4 0.3 L 30 / 0.2 0.25 0 q| nk n /|q 0.2 20 0 N 0 0.01 0.02 0.15 1/L 10 0.1 0.05 0 0 0 π/8 π/4 0 π/4 π/2 3π/4 π 0 π/4 π/2 3π/4 π |q| |q| |k| FIG.3: Density structurefactor N dividedby|q|calculated FIG.4: Variationalresultsforthemomentumdistributionn q k with variational Monte Carlo (left panel) and GFMC (right in 1D for L=60 (squares), 100 (circles), and 150 (triangles) panel) for different U/t and L = 60. From top to bottom across the transition (U/t=2.4 and 2.5). Inset: Size scaling U/t=1.6, 1.8, 2, 2.2, 2.4, 2.5, 3, and 4. of the condensate fraction n0/L for U/t = 2.4 (upper curve) and U/t=2.5 (lower curve). us now concentrate on the insulating phase. Here, the 1 Jastrow wave function (5) can be mapped onto the par- 40 tition function of an effective classicalCGM and β plays L the role of the inverse classicaltemperature β =π/T . 0.8 x eff s 20 In 1D the CGM is in the confined phase for any finite D temperature, with exponential correlations.31 This out- 0.6 1 0 comeisconsistentwiththefactofhaving,inthequantum s 1.5 2 2.5 3 D model, a finite gap in the excitation spectrum. Remark- U/t 0.4 0.5 ably closeto the Motttransitioninthe insulating phase, the value of β obtained from the optimized Jastrow po- tentialisverysmallandapproachestozerowhenU Uc 0.2 0 fromabove. Since the correlationlengthofthe 1DC→GM 0 0.02 0.04 1/L divergesforβ 0,ournumericaloutcomegivesastrong 0 → indicationinfavorofacontinuoustransitionbetweenthe 1 1.5 2 2.5 3 superfluid and the Mott insulating phase. U/t Let us now analyze the transition by considering the density structure factor N . In the small-q regime, we q FIG. 5: Superfluid stiffness D calculated by GFMC as a s can generally write that function of U/t for different sizes in the 1D boson Hubbard N =γ q +γ q2+O(q3), (16) model. Lowerinset: SizescalingofDsfordifferentU/t(same q 1 2 | | values of the main panel, with U/t increasing from top to where γ1 and γ2 depend upon the Jastrow parameters. bottom). Upper inset: Ds ×L as a function of U/t. The In analogy with spin systems, we have that γ = v χ, point where the different curves cross marks the transition 1 c withv andχthesoundvelocityandthecompressibility, point. c respectively. The fact of having γ =0 in the insulating 1 regime indicates that this state is incompressible (i.e., The superfluid-insulator transition can be also easily χ=0). From Fig. 3, we obtain that γ has a very sharp 1 detected by considering the momentum distribution: crossover from a finite value to zero across the transition, suggestive of a true jump in the thermodynamic Ψ b†b Ψ limit. This outcome is consistent with the fact that the n = h J| k k| Ji, (17) k compressibilityalsohasafinite jumpinthe1Dquantum ΨJ ΨJ h | i phase transition.34 Moreover, just above U in the insu- c latingregime,γ isverylargeforboththevariationaland whereb† isthecreationoperatorofabosonofmomentum 2 k theGFMCcalculations,indicatingthepeculiarcharacter k. This quantity has a radicallydifferent behaviorbelow ofthe1Dtransition. Byusingthesmall-qbehaviorofthe and above the transition: In the superfluid phase, it has density structure factor, the Mott transition can be lo- acuspatk =0,althoughthereisnocondensatefraction, catedatU /t 2.45 0.05forthevariationalcalculation, i.e., n /L 0 for L , while in the insulating phase c 0 ∼ ± → →∞ whereas the GFMC approach gives U /t 2.1 0.1. it is a smooth function of the momentum k, see Fig. 4. c ∼ ± 6 Finally, we want to conclude the 1D part by consider- 0.3 ing the superfluid stiffness Ds. In analogy to what has Jastrow been done by Pollock and Ceperley at finite tempera- Gutzwiller+MB 0.25 ture,35 this quantity can be also calculated directly at Jastrow+MB zero temperature by using the GFMC and the so-called 0.2 winding numbers (see Appendix) 0 E / 0.15 E Ψ0 W~ (τ)2 Ψ0 ∆ D = lim h || | | i, (18) s τ→∞ D Lτ 0.1 where Ψ istheground-statewavefunctionobtainedby 0.05 0 | i GFMC, D is the dimension of the system and W~ (τ) = [~r (τ) ~r (0)],~r (τ)beingthepositionofthei-thpar- 0 i i − i i 0 2 4 6 8 10 12 14 16 ticleafterevolvingitforadiffusiontimeτ fromtheinitial pPosition~r (0). The diffusion process must be done with- U/t i outconsideringperiodicboundaryconditions,namelyby increasing or decreasing the values of the coordinates of FIG. 6: Accuracy of different variational wave functions as function of U/t for the 10×10 cluster and 100 bosons. The a particle that crosses the boundaries of the cluster. In symbols are thesame as in Fig. 1. this way, non-zero winding numbers across the lattice can be detected. It should be stressed that, exactly at zero temperature, Ds can only give information about 25 the presence of a gap in the excitation spectrum, and, U/t=11,12,16 therefore, it can discriminate between conducting and 0 20 insulating phases. Our results showthat D is finite and s B -0.5 M large in the weak-coupling regime, whereas it vanishes g for U/t & 2.1, see Fig. 5. Again, in analogy with spin 2 15 -1 q| 8 12 16 20 systems andfromgeneralscaling argumentsvalid for1D x | U/t boson models, we expect a jump at the transition, that q 10 v howeveris veryhard to detect by consideringfinite clusters.34 Nevertheless, an accurate value of the transition point is obtained from the size scaling of D L, see 5 s × Fig. 5. 2D U/t=6,8,9,10 0 0 π/4 π/2 3π/4 π |q| IV. THE 2D HUBBARD MODEL FIG.7: Jastrowparametersv multipliedbyq2 asafunction q Letusnowturntothe2DHubbardmodelandconsider of |q| [along the (1,0) direction] for 20×20 (circles), 26×26 square clusters with L = l l sites and N = L bosons. (squares), and 30×30 (triangles) clusters. Inset: The many- × Theaccuracyofthethreewavefunctions(5),(6)and(8) body variational parameter g as a function of U/t. MB arereportedinFig.6. Thesituationisdifferentfromthe previous 1D case: The Gutzwiller state with the many- body term, which in 1D is not accurate for large U/t, in in the confined phase. In reality, on the sizes available 2D becomes competitive with the Jastrowwavefunction to our numerical approach, we can not firmly establish for U/t & 14. Moreover, the presence of the many-body whether a classical potential log(q)/q2 has indeed ex- − termstronglyimprovestheaccuracyoftheJastrowstate ponentialdecayingcorrelationfunctions. Nevertheless,it as soon as U/t & 8. Then, in the following, we will is remarkable and very encouraging that the variational consider the state (8) for both the variational and the optimization leads to such a singular v . q GFMC calculations. The evidence that the change in behavior of v does q In Fig. 7, we show the behavior of the optimized v as correspond to the superfluid-insulator transition comes q a function of the interaction strength. Similarly to what also by the small-q limit of the structure factor N . In q happens in the 1D system, we obtain that v α/q for Fig. 9, we show the results for the variational and the q U/t.10.5,whilev isbestfittedbyv lo∼g(q)/|q2| for GFMC calculations as a function of U/t. In both cases, q q ∼− U/t & 10.5 (see Fig. 8), corresponding to the superfluid we find a different small-q behavior for large and small and the Mott insulating phase, respectively. As we an- couplings. In the variationalcalculations, for U/t.10.3 ticipated in section IIB, the Jastrow potential in the in- the structure factor behaves as N γ q , while for q 1 sulatingphaseismoresingularthan1/q2,suggestiveofa U/t & 10.3 we get N γ q2. There∼fore,|a|t the varia- q 2 ∼ classicalmodelwithboundchargesbutpresumablywith- tional level, the superfluid-insulator transition is located out the power-law correlations displayed by the CGM around U /t 10.3 0.1; this value is slightly smaller c ∼ ± 7 25 2 π) 20 1.5 / q g( s v 2/lo 15 U/t=11,12,16 0 1 q| c | x 0.5 q v 10 - 0 0 1 2 3 4 5 6 7 8 5 0 π/4 |q| π/2 3π/4 U/t FIG. 8: Jastrow parameters v multiplied by q2/log(q/π) as FIG. 10: GFMC results for the sound velocity v obtained q c a function of |q| [along the (1,0) direction] for different U/t through a finite size scaling of the ground-state energy, see and the same sizes as Fig. 7. Eq.(19). The line is a guideto the eyes. VMC GFMC 80 0.02 0.07 U/t=10 U/t=12 0.06 60 0.2 0.015 0.05 L / 0.1 0 q| 0.04 k 40 n /|q 0.01 n 0 N 0.03 0 0.01 0.02 0.005 0.02 20 1/L 0.01 0 0 0 0 π/4 π/2 3π/4 π 0 π/4 π/2 3π/4 π 0 π/8 π/4 |q| |q| |k| FIG. 9: Left panel: Density structure factor N divided by q |q|calculatedbythevariationalMonteCarlofordifferentU/t FIG. 11: Variational results for the momentum distribution and L= 20×20. From top to bottom U/t = 10, 10.2, 10.4, n in 2D for 10 × 10 (squares), 16 × 16 (circles), 20× 20 k 11, and 12. Right panel: The same for the GFMC on the (upwardtriangles),and26×26(downwardtriangles)clusters L = 16×16 cluster. From top to bottom U/t = 8, 8.2, 8.4, with U/t = 10 and 12. Inset: Size scaling of the condensate 8.6, and 8.8. fractionn0/LforU/t=10(uppercurve)andU/t=12(lower curve). than the one obtained by the simple Gutzwiller wave function, for which U /t 11.65.14 The critical value of finite-size scaling of the exact ground-state energy by c ∼ the on-site interaction is rather different within GFMC, c v 0 s yielding Uc/t ∼ 8.5±0.1, in close agreement with the ǫ0(L)=ǫ0(∞)− l3 , (19) value found in the literature.9 Let us focus more deeply on the behavior of the structure factor N . Approach- whereǫ (L)istheground-stateenergypersiteforaclus- q 0 ingthetransitionfromtheweak-couplingregion,γ goes ter with L = l2 sites, ǫ ( ) is the extrapolated value 1 0 ∞ smoothly to zero, in contrast with the results of the 1D in the thermodynamic limit, and c is a given model- 0 model, where we observed an abrupt jump. Also con- dependent constant. Our results, shown in Fig. 10, trastingthe1Dcase,γ isfoundtobefiniteandcontinu- clearlyindicate that v staysfinite acrossthe superfluid- 2 s ous across the transition. These results, based upon the insulatortransition,thusimplyingavanishingcompress- variational wave function, are confirmed by the GFMC ibility when approaching the insulator. calculations, see Fig. 9. The vanishing linear coefficient The fingerprint of this transition is also given by the of N can be ascribed either to v or to χ. In order momentumdistribution, see Fig.11. For this quantity, a q s to clarify which one of these quantities goes to zero at strikingdifferenceisobservedbelowandaboveU . Inthe c the transition, we extractthe soundvelocity v from the formercase,acusp-likebehaviorwithafinitecondensate s 8 0.6 6x6 100 12x12 16x16 80 U/t=20,22 0.4 0.6 3q| 60 Ds x | 0.4 v q 40 U/t=8,12,16 0.2 Ds 0.2 20 3D 0 0 1/L 0.02 0 0 0 π/4 π/2 3π/4 π 6 6.5 7 7.5 8 8.5 9 |q| U/t FIG. 14: Jastrow parameters v multiplied by |q|3 as a func- q tion of |q| for 8×8×8 (circles), 10×10×10 (squares), and FIG.12: SuperfluidstiffnessD asafunctionofU/tfordiffer- s 12×12×12 (triangles) clusters [along the(1,0,0) direction]. entsizes. Inset: SizescalingofD fordifferentU/t(thesame s values of the main panel are reported, with U/t increasing from top to bottom). 0.0025 0.002 2 q| /| 0.0015 q 0.35 N Jastrow 0.001 0.3 Gutzwiller+MB 6x6x6 3D Jastrow+MB 8x8x8 0.25 0.0005 10x10x10 VMC 12x12x12 /E 0.2 00 π/4 π/2 3π/4 E ∆ 0.15 0.04 1e-03 0.1 0.03 2 0.05 q| /|q 0.02 5e-04 0 N 0 5 10 15 20 25 30 0.01 β =20 β =55 U/t 3D 3D 0 0e+00 0 π/4 π/2 3π/4 0 π/4 π/2 3π/4 FIG. 13: Accuracy of different variational wave functions as |q| |q| afunctionofU/tforthe10×10×10clusterand1000bosons. The symbols are thesame as in Fig. 1. FIG. 15: Upper panel: Variational results for the density structure factor N for 3D and U/t = 20. Lower panels: N q q for non-optimized wave functions with vq ∼β3D/|q|3 for two valuesof β3D. fraction is found, whereas in the latter case a smooth behavioris detected, with a vanishing n /L. Notice that 0 a vanishing condensate fraction in the thermodynamic limit immediately follows from N q2 by using the f- V. THE 3D HUBBARD MODEL q ∼ sum rule derived in Ref. 36. Let us move now to the 3D case. Here, we mostly re- Finally,wereportinFig.12thestiffnessD ,calculated strict our analysis to the variational method that allows s by using GFMC. In 2D, we obtain a different behavior ustoassessratherlargesizes. Westartbyconsideringthe withrespecttothe1Dcase,whereafinitejumpisrather accuracy of the different variational states as a function clear at the superfluid-insulator transition. Indeed, the of the interaction U/t, see Fig. 13. It turns out that in evaluation of the stiffness for different sizes confirms the the large-U regionthe mostaccuratestate containsboth absenceofthejumpin2D,asexpectedforasecondorder the long-range Jastrow term and the many-body term, phase transition.9 as it occurs in 2D. Moreover, in analogy to 1D and 2D, 9 the presence of a phase transition upon increasing U/t 25 is clearly signaled by the sudden change in the small-q 0 behavior of the Jastrow factor, see Fig. 14. Here, its V/t=9,10 B-0.4 M behavior changes drastically from v α/q to a more 20 g-0.8 q ∼ | | diverging vq 1/q 3 in the insulating regime. In partic- 2 4 6 8 10 ∼ | | 2 ular, the sudden change of behavior allows us to locate 3/ 15 V/t=8 V/t the transition around U /t 18, which is close to the q| critical value of recent Mcont≃e Carlo simulations.10 x | q 10 Focusing on the large-U region of the phase diagram, v V/t=6 let us discuss the implications ofthe singularv 1/q 3 Jastrow potential. First of all, following the saqm∼e ar|gu|- 5 V/t=4 ments of Ref. 32, we can assert that the strong singular V/t=2 character of the Jastrow potential is able to empty the 0 condensate, whereas a less singular v 1/q2 would not 0 π/4 π/2 3π/4 π q leadton 0inthethermodynamicli∼mit. Remarkably, |q| 0 → even though v 1/q 3, the structure factor in the in- q sulator has the c∼orrec|t|Nq q2 behavior, see Fig. 15. In FIG.16: Jastrowparametersvq multipliedby|q|3/2asafunc- turn,thisimpliesthatEq.(∼15)doesnothold. Inorderto tion of |q| for the potential of Eq. (3) with different V/t and prove more firmly that v β /q 3 can indeed lead to 20×20 (circles), 26×26 (squares), and 30×30 (triangles) q 3D ∼ | | clusters. Inset: The many-body variational parameter g N q2, we have calculated the density structure factor MB q ∼ as a function of V/t. with a non-optimized wave function of the form (5) with v β /q 3 and for different values of β . As shown q 3D 3D in F∼ig. 15,|f|or small β we obtain N q 3, implying 3D q ∼ | | that Eq. (15) is qualitatively correct in this case. How- 0.08 V/t=4 ever, above a critical β∗ , the small-q behavior of the 3D density structure factor turns into N q2, signaling a q ∼ remarkablebreakdownofEq.(15). FromFig.14itturns 0.06 2 out that the optimal value of β that we get variation- 3/ V/t=6 3D q| aβl∗ly,actotnhfiermsuinpgerflthuaidt-iNnsulatqo2r tinranthsietiionnsuislaltainrggerphtahsaen. / | 0.04 3D q ∼ Nq Most importantly, the change of behavior as a function V/t=8 of β is consistent with the binding-unbinding phase- 3D 0.02 transitionrecentlyuncoveredinaclassical3Dgaswitha V/t=9,10 1/q 3potential.33 Onceagain,asin1Dand2D,theMott | | insulating wave function in 3D is found to be closely re- 0 0 π/4 π/2 3π/4 π latedtothelow-temperatureconfinedphaseofaclassical model, where opposite charges tightly bound. |q| FIG. 17: Variational results for the density structure factor N(q)dividedby|q|3/2 forthepotentialofEq.(3)withdiffer- VI. 2D SYSTEMS WITH LONG-RANGE ent V/t and 20×20 (circles), 26×26 (squares), and 30×30 INTERACTION (triangles) clusters. In the previous paragraphs, we have shown that, in the Hubbard model (1), the superfluid regime can be be described by the variationalwave function (8) with a described by a long-range Jastrow wave function with Jastrow potential v β/q2. q ∼ v α/q . By increasing the on-site interaction U, For that purpose, we consider the more general Hub- q ∼ | | our variational approach describes a continuous transi- bardmodelofEq.(2)withalong-rangeinteractionΩ(r). tion to an insulating phase that corresponds to the con- LetusstartbyconsideringtherealisticcaseofaCoulomb fined phase of a classical model of interacting particles potential,namelyΩ(r) 1/rthatin2DleadstoEq.(3). ∼ with opposite charge. In particular, we have found ev- Then, we can vary its strength V to drive the system idences that the 2D Mott insulating wave function cor- across a superfluid-insulator transition. The small-q be- responds to a classical model with log(q)/q2 potential havior of the optimized Jastrow parameters is shown in − rather than to a 2D CGM with potential β/q2 in the Fig. 16. For V/t . 8, we obtain that v 1/q 3/2. q confined phase (β > β∗ 4π). This result, as we dis- On the contrary, for larger values of the inter∼actio|ns|, v q ≥ cussed, has a physical importance since the confined 2D becomes moresingularandinparticularis best fittedby CGM has power-law correlations,not possible in a Mott log(q)/q2,justlikewhatwefoundforshort-rangeinter- − insulatingphase. Nevertheless,itwouldbeinterestingto action. The structure factor in the weak-coupling phase search for bosonic Hamiltonians whose ground state can behaves as N q 3/2, see Fig. 17, compatible with a q ∼ | | 10 β β∗ at the transition. This indicates that the CGM 25 → 0 thatcorrespondstothevariationalwavefunctionisinthe 20 V/t=20 MB-0.4 pcolarsrmelaatpiohnasfeu,nwctitiohnesx.poInndeneetida,ltdheecasytirnugctduernesiftayc-dtoernsbitey- g V/t=18 -0.8 haveslike N q2 for allthe coupling strengths V/t, see q ∼ 2q| 15 4 8 12 16 20 Fig. 19. We believe that the weak-coupling phase has to | V/t=16 V/t beidentifiedwiththealgebraiclong-rangeorderedphase x found at high density in the continuum model by Magro q 10 v V/t=14 and Ceperley.21 This phase is characterized by absence V/t=10 of condensate but by a power-law decay of the single- 5 particle density matrix. On the contrary, the strong- V/t=6 coupling phase with Jastrow potential v log(q)/q2 q ∼ − 0 mustcorrespondtoagenuineMottinsulatorwithallcor- 0 π/4 π/2 3π/4 π relation functions decaying exponentially. |q| FIG.18: Thesame asin Fig. 16,for thepotentialof Eq.(4). VII. CONCLUSIONS Wehaveshownthatthelong-rangeJastrowwavefunc- tion gives a consistent picture of the superfluid-insulator V/t=6 transition of the bosonic Hubbard model in all spatial 0.06 dimensions and also in the presence of long-range interaction. V/t=10 2 Inonedimensionthevariationalresultsarecompatible q| 0.04 with a Berezinskii-Kosterlitz-Thouless phase transition /|q V/t=14 between the quasi-long-range ordered gapless phase and N the Mott insulator. From the point of view of the vari- 0.02 V/t=16 ational wave function, the gapless phase is characterized by a Jastrow potential v α/q , while the insulating q V/t=18 one by v 1/q2. ∼ | | q 0 V/t=20 In two ∼dimensions we have evidences of a second- 0 π/4 π/2 3π/4 π order phase transition between a superfluid phase and |q| aMottinsulator. Here,thesuperfluidwavefunctionstill has v α/q , compatible with the existence of sound q FIG.19: Thesame asin Fig. 17,for thepotentialof Eq.(4). modes,∼while| t|he insulating wave function is characterized by a more singular v log(q)/q2. This singu- q ∼ − lar behavior in the Mott phase does not change even if superfuid phase with 2D plasmons. We mention that long-range interaction is considered. For instance, for a a similar behavior has been found in continuum models CoulombinteractionΩ(r) 1/r,theJastrowpotentialin at high densities both analytically37 and numerically.38 thesuperfluidphasechang∼esintov 1/q 3/2,compati- q ∼ | | In the strong-coupling regime, the insulating behavior ble withthe existence of2Dplasmons,yetthe insulating N q2 is recovered. These results are confirmed by wave function has still v log(q)/q2. For an interac- q q ∼ ∼− GFMC (not shown), though the critical V/t is slightly tion Ω(r) log(r), we observe a transition between an decreased, i.e., V /t 7. In this case, as before with algebraic l∼ong-range ordered phase,21 characterized by c ∼ short-rangeinteraction, the optimized Jastrow never be- v 1/q2, and a Mott insulating phase, once again with q haves as the potential of a CGM. v ∼ log(q)/q2. q ∼− Therefore, let us turn to the more singular interaction Inthreedimensions,theMotttransitionasrevealedby given by Eq. (4), leading to Ω(r) log(r). In this thebehavioroftheJastrowpotentialbecomesmuchmore ∼ − case the potential in q-space is given by Ω(q) 1/q2 evident. Asusual,thesuperfluidphasehasv α/q . In q ∼ ∼ | | and we expect, similarly to what happens in the contin- thiscase,however,theMottinsulatingwavefunctionhas uumforhighdensity,21 thatintheweak-couplingregime a much more singular v 1/q 3, still with a structure q also v β/q2. Indeed, as shown in Fig. 18, this is the factor that correctly beha∼ves a|s|N q2. q ∼ q ∼ case for V/t . 16. Above this value, v becomes once Thisworkwasmainlyfocusedonbosons,butasimilar q again of the form log(q)/q2, just like in all previous variationalwavefunctioncanbeappliedalsotofermionic − examples. We note that the values of β extracted in the models, with the same key ingredient, i.e., a long-range weak coupling phase seem to be all smaller than β∗ of Jastrowfactorthatdrivesthe metal-insulatortransition, the 2D Berezinskii-Kosterlitz-Thouless phase transition, in spite of the uncorrelated wave function being metal- although we can not exclude by the numerical data that lic.15,17 According to this picture, the Mott insulating