Phase diagram and thermodynamics of the three-dimensional Bose-Hubbard model

0.27 MB·English

by  B. Capogrosso-Sansone

Phase diagram and thermodynamics of the three-dimensional Bose-Hubbard model
B. Capogrosso-Sansone,1 N.V. Prokof’ev,1,2,3 and B.V. Svistunov1,3
7
0 1Department of Physics, University of Massachusetts, Amherst, MA 01003, USA
0 2BEC-INFM, Dipartimento di Fisica, Universita´ di Trento, Via Sommarive 14, I-38050 Povo, Italy
2 3Russian Research Center “Kurchatov Institute”, 123182 Moscow, Russia

n We report results of quantum Monte Carlo simulations of the Bose-Hubbard model in three
a
dimensions. Criticalparametersforthesuperfluid-to-Mott-insulatortransitionaredeterminedwith
J
significantly higher accuracy than it has been done in the past. In particular, the position of the

9 critical pointat fillingfactor n=1is foundtobeat (U/t) =29.34(2), andtheinsulatinggap ∆ is
c
measured with accuracy of a few percent of thehopping amplitudet. Weobtain theeffective mass
]
ofparticleandholeexcitationsintheinsulatingstate—withexplicitdemonstrationoftheemerging
r
e particle-hole symmetry and relativistic dispersion law at the transition tip—along with the sound
h velocityinthestronglycorrelatedsuperfluidphase. Theseparametersarethenecessaryingredients
t to perform analytic estimates of the low temperature (T ≪ ∆) thermodynamics in macroscopic
o
. samples. Wepresentaccuratethermodynamiccurves,includingtheseforspecificheatand entropy,
at for typical insulating (U/t = 40) and superfluid (t/U = 0.0385) phases. Our data can serve as a
m basis for accurate experimental thermometry, and a guide for appropriate initial conditions if one
attemptsto use interacting bosons in quantuminformation processing.
-
d
n PACSnumbers: 03.75.Hh,03.75.Lm,75.40.Mg
o
c
[ I. INTRODUCTION mensurate filling factor undergoes a superfluid-to-Mott

insulator(SF-MI)quantumphasetransition. Theground
1
stateofMIcanbe usedinquantuminformationprocess-
v Inthepastdecade,stronglycorrelatedlatticequan-
ing to initialize alargesetofqubits (the mainremaining
8 tum systems have been attracting a lot of interest and
7 effort. Remarkably, simple yet nontrivial models which challengeis in addressingsingle atoms to build quantum
1 gates, see Ref. [1] and references therein). Atomic sys-
contain most of the important many-body physics and
1 tems in optical lattices have the advantage of being well
known in the theory community for many years can be
0 isolatedfromtheenvironment. Thisresultsinarelatively
now realized and studied experimentally. For the first
7
long decoherence time of the order of seconds [1] and
0 time theoretical predictions and experimental data for
therefore the possibility of building long-lived entangled
/ strongly correlated states can be directly tested against
t many body states. These properties make MI ground-
a each other in the ideal setup when all model ingredients
m are known and controlled. states good candidates for building blocks of a quantum
computer. Another possible application is in interfero-
- Experimentally,latticesystemsarerealizedbytrap-
d metricmeasurements[11]. Ithasbeenargued[12,13,14]
ping atoms in an optical lattice, a periodic array of po-
n that using the superfluid-to-Mott-insulatorphase transi-
tential wells resulting from the dipole coupling of the
o tion to entangle and disentangle atomic Bose-Einstein
c atoms to the electric field of the standing electromag-
condensateonecangobeyondtheHeisenberg-limitedin-
: netic wave produced by a laser. Optical lattices are a
v terferometry.
very powerful and versatile tool. By changing the laser
i A system of bosons with short-range repulsive pair
X parameters and configuration, the properties and geom-
interaction trapped in an optical lattice is described by
r etry of the optical lattice can be finely tuned [1]. Ul-
a the Bose-Hubbard Hamiltonian:
timately, this results in the possibility of controlling the
Hamiltonianparametersandexploringvariousregimesof H = −t b†b + U n (n −1)− µ n , (1)
interest. In particular, ultra-cold Bose atoms trapped in i j 2 i i i i
<Xij> Xi Xi
an optical lattice are an experimental realization of the
Bose-Hubbard model. The model has been studied in whereb† andb arethebosoniccreationandannihilation
i i
the seminal paper by Fisher, Weichman, Grinstein, and operators on the site i, t is the hopping matrix element,
Fisher,Ref.2,anditsphysicalrealizationwithultra-cold U istheon-siterepulsionandµ =µ−V(i)isthesumof
i
atoms trapped in an optical lattice has been envisioned thechemicalpotentialµandtheconfiningpotentialV(i).
in Ref. [3]. Few years later, the Bose-Hubbard system In what follows, we consider bosons in the simple cubic
wasproducedinthe laboratory [4]. Since then, the field lattice. Atzerotemperatureandintegerfillingfactor,the
remains very active [5, 6, 7, 8, 9, 10], not only because competitionbetweenkineticenergyandon-siterepulsion
theoretical predictions and experimental techniques still inducestheMI-SFtransition. Whentheon-siterepulsion
have to be substantially improved to claim the quanti- isdominating,t/U ≪1theatomsaretightlylocalizedin
tative agreement, but also because of the new physical the MI ground state which is well approximated by the
applications. product of local (on-site) Fock states. The Mott state
At zero temperature, a systemof bosons with com- is characterized by zero compressibility originating from
2
an energy gap for particle and hole excitations. When and hole excitations inside the insulating phase. Close
the hopping amplitude is increased up to a certain criti- to the MI lobe tip, the data for the dispersion of the el-
calvalue(t/U) ,particledelocalizationbecomesenerget- ementary excitations are fitted by the relativistic law, in
c
ically more favorableand the system Bose condenses. In agreementwiththe theory(this alsoallowsus to extract
the chemicalpotentialvs. hopping matrixelement plane the value of the sound velocity in the critical region).
(energies are scaled by U), the T =0 phase diagramhas In the Mott state, the knowledge of gaps and effective
a characteristic lobe shape [2], see also Fig. 3 below, masses is sufficient to calculate the partition function in
with the MI phase being inside the lobe (there is one the low temperature limit analytically and to make reli-
lobe for each integer filling factor). The most interesting able predictions for the system entropy.
regioninthephasediagramisthevicinityofthelobetip, For such applications of the system as quantum in-
(µ = µc, t = tc), corresponding to the MI-SF transition formation processing and interferometry, controlling the
in the commensurate system. For other values of µ or temperature is of crucial importance. Most applications
t, the SF-MI criticality is trivial and corresponds to the are based on the key property of the good insulating
weaklyinteractingBosegasatvanishingparticle density state, which is small density fluctuations in the ground
[2]. Itisstraightforwardlydescribedprovidedtheparticle state. At zero temperature fluctuations are of quan-
(hole) effective mass is known. If, however, one crosses tum nature and can be efficiently controlled externally
the MI-SF boundary at constant commensurate density through the t/U ratio. At finite temperature, fluctua-
(this is equivalent to going through the tip of the lobe tions are enhanced by thermally activated particle-hole
atafixedchemicalpotential)thelong-waveactionofthe excitations. Only when the temperature is much smaller
system becomes relativistic and particle-hole symmetric. than the energy gap, the number of excitations is expo-
Nowthe phasetransitionis inthe four-dimensionalU(1) nentially small. Up to date, there are no available ex-
universality class [2]. It is worth emphasizing that here perimental techniques to measure the temperature of a
we havea unique opportunity ofa laboratoryrealization strongly interacting system. For weakly interacting sys-
of the non-trivial relativistic vacuum, a sort of a “hy- tems, the temperature can be extracted in a number of
drogen atom” of strongly-interacting relativistic quan- ways, e.g. from the interference pattern of matter waves
tum fields. Approaching the critical point from the MI [17] or the condensate fraction observed after the trap
side, one deals with the vacuum that supports massive is released and the gas expands freely [18] —these prop-
bosonic particlesand anti-particles(particlesand holes). erties are directly related to the momentum distribution
On the other side of the transition, the SF vacuum sup- function n(k). For strongly interacting systems, both
ports massless bosons (phonons) that do not have an temperatureandinteractionareresponsibleforfillingthe
anti-particle analog. In principle, one cansystematically highermomentumstates,whichmakesithardtoextract
studyuniversalmultiparticlescatteringamplitudesofthe temperature using absorption imaging techniques.
relativisticquantumfieldtheoryintheultra-cold”super-
The results presented in this paper can be used
collider”!
to perform accurate thermometry. Typically, the ini-
tial temperature, T(in), (before the optical lattice is adi-
The present study is focused on the three-
abatically loaded) is known. By entropy matching one
dimensional (3D) system. To the best of our knowledge,
can easily deduce the final temperature of the MI state,
previoussystematicstudiesofthe3Dcasewerelimitedto
T = T(fin), provided the entropy of the MI phase is
the mean-field (MF) [2] and perturbative methods [15].
known. To this end we have calculated the energy, spe-
In Ref. [15], the authors utilized the strong-coupling ex-
cific heat and entropy of the system in several impor-
pansion to establish boundaries of the phase diagram in
tant regimes which include MI and strongly correlated
the(µ/U, t/U)plane. Thisapproach,basedonthesmall
SF phases. These data can be used to suggest appropri-
ratio zt/U ≪1, where z =6 is the coordination number
ate initial conditions which make the Bose-Hubbardsys-
for the simple cubic lattice, works well only in the MI
tem suitable for physical applications, such as the ones
phase in the region far from the tip of the lobe, where
described above.
the insulating gap is larger then hopping, ∆/zt > 1.
Close to the critical region, where ∆/zt∼1, the strong- Anotherinterestingquestionconcernsthe natureof
coupling expansion is no longer valid. We present the inhomogeneous states in confined systems when the MI
results of large-scale Monte Carlo (MC) simulations of phase is formed in the trap center. The confining po-
the model (1) by worm algorithm [16]. With precise tential provides a scan in the chemical potential of the
data for the single-particle Green function, we are able phasediagramatfixedt/U [3]. As one movesawayfrom
to carefully trace the critical and close-to-criticalbehav- the trap center the system changes its local state. At
ior of the system, and, in particular, produce an ac- zero temperature, the density profile of the system can
curate phase diagram in the region of small insulating be read (up to finite-size effects) from the ground state
gaps∆≪t. Thoughthe correspondingparameterrange phase diagram. At finite temperature, this is no longer
is quite narrow, it is crucial to clearly resolve it to re- possible. In particular, the liquid regions outside of the
veal the emerging particle-hole symmetry and relativis- MI lobes could be normal or superfluid, depending on
tic long-wave physics at the tip of the MI-SF transition. temperature.
We also present data for the effective masses of particle So far experimental results have been interpreted
3
by assuming that liquid regions are superfluid, but there
T/t
were no direct measurements or calculations to prove NL
that this was the case. [Part of the problem is that 6
5.591
absorption imaging is sensitive only to n(k), which is
the Fourier transform of the single-particle density ma-
4
trix in the relative coordinate. All parts of the system SF
contribute to n(k) and it is hard to discriminate where
the dominantcontributioncomes from.] Itis almostcer- 2
tain that T(fin)/T(fin) of the strongly correlated system
c
is higher then T(in)/T(in). Indeed, since the entropy of MI
c 0
MI at ∆ ≫ T is exponentially small, most entropy will 0 10 20 30 U40/t
be concentrated in the liquid regions. At this point we
(U/t) =29.36(2)
notice that the transition temperature in the liquid is c
suppressedrelativeto thenon-interactingBosegasvalue
T(0) ≈ 3.313n2/3/mby both (i) effective massenhance- FIG. 1: (Color online). Finite-temperature phase diagram
c
ment in the optical lattice, m → 1/2ta2 (here a is the at filling factor n = 1. Solid circles are simulation results
latticeconstant),and(ii)strongrepulsiveinteractionsin (the line is a guidance for the eye), error bars are plotted.
the vicinity of the Mott phase, in fact, T → 0 at the T = 5.591t is the critical temperature of the ideal Bose gas
c
with the tight binding dispersion relation. At finite, but low
SF-MIboundary. ItseemsplausiblethattheMIphaseis
enough temperature,theMI domain isloosely definedasthe
alwayssurroundedbyabroadnormalliquid(NL)region.
part of the phase diagram to the right of the gray line. The
Itmayalsohappenthatsuperfluidityiscompletelyelim-
rest of the non-superfluid domain is referred to as normal
inated in the entire sample in the final state. [Strictly
liquid (NL).
speaking, at T 6= 0 the MI and NL phases are identical
in terms of their symmetries and are distinguished only
quantitatively in the density of particle-hole excitations, sion, work well in the region where zt/U ≪ 1 and the
i.e. in the Hamiltonian (1) the finite-temperature MI is system is deep in the MI phase. Under these conditions
continuously connected without phase transition to NL, thekineticenergytermintheHamiltoniancanbetreated
see Fig. 1. For definiteness, we will call NL a normal perturbatively and the unperturbed ground state is a
finite-T state which is superfluid at T = 0 for the same productoflocalFockstates. InRef.[15]theauthorscar-
set of the Hamiltonian parameters.] Fig. 1 shows the ried out an expansion, up to the third order in zt/U, for
finite-temperature phase diagram for filling factor n=1 the SF-MI boundariesandestimatedpositions ofcritical
(we will discuss how we determine the critical temper- points at the tips of the MI lobes (by extrapolating re-
ature in Sec. III). The critical temperature goes to zero sultstothe infinite expansionorder). Theirresultsagree
sharply,whileapproachingthecriticalpoint. Inthelimit with the mean field solution calculated in Ref. [2], when
of U → 0 the critical temperature is slightly above the the latter is expanded up to the third order in zt/U and
ideal-gasprediction(T =5.591twascalculatedusingthe the dimension of space goes to infinity. As already men-
tight binding dispersion relation), as expected (see, e.g. tioned, this approach starts failing when ∆ ∼ zt. Using
Ref. [19]). MC techniques we were able to calculate critical param-
The paper is organized as follows: in Sec. II we eters and predict the position of the diagram tip with
present results for the ground state phase diagram and much higher accuracy: with the worm algorithm (WA)
effective mass of particle (hole) excitations, at integer approach the energy gaps can be measured with preci-
filling factor n = 1. In Sec. III we investigate the ther- sion of the order of 10−2t [20]. The simulation itself is
modynamic properties of the system. We present data basedontheconfigurationspaceoftheMatsubaraGreen
for energy, specific heat and entropy and calculate the function
final temperature of the uniform and harmonically con-
fined system in the limit of large gaps. For the case of G(i,τ) = hT b†(τ)b (0)i, (2)
τ i 0
trapped system, we also determine the state of the liq-
uid at the perimeter of the trap. Brief conclusions are which is thus directly available. We utilize the Green
presented in Sec. IV. functiontodeterminedispersionrelationsforparticleand
hole excitations at smallmomenta [fromthe exponential
decayofG(p,τ) withthe imaginarytime] whichdirectly
give us the energy gap and effective masses.
II. GROUND STATE PROPERTIES
RecallthatinthemomentumspacetheGreenfunc-
tion of a finite size system G(p,τ) is different from zero
This section deals with the results of large-scale only for p = p = 2π(m /L , m /L , m /L ), where
m x x y y z z
Monte Carlo simulations for the ground state phase di- L is the linear system size in direction α (we per-
α=x,y,z
agram of the Bose-Hubbard system in three dimensions. formedallsimulationsinthecubicsystemwithL =L),
α
Analytical approaches, e.g. the strong coupling expan- and m = (m ,m ,m ) is an integer vector. Using
x y z
4
µ/U
0 50.6
∆(t) 0.45
50.4 0.8 0.40
-2
0.35
50.2
)) 0.00 0.05 1/L 0.10 0.6 0.30
τ 0.032 0.036
0, -4
=
p
( 0.4
G
(
n
l -6
0.2 Freericks & Monien
-8
0.0
26 28 30 32 34 36 38 40 0.00 0.01 0.02 0.03 0.04 0.05
τ(t) t/U
FIG. 3: (Color online). Phase diagram of the first MI-SF
FIG. 2: (Color online). Zero-momentum Green function in
lobe. Numerical data are shown by open circles. The error
the Mott phase with the chemical potential µ/U = 0.809,
bars are shown but are barely visible even in the inset. The
slightly below the upper phase boundary. Here we show
data for the system with N = 103 bosons at U/t = 70 and dashed lines and the square represent results of Ref. [15] for
the strong coupling expansion and the extrapolated position
T/t=0.025. Intheinset weplot theenergygap ∆ forlinear
of the diagram tip, respectively.
system sizes L=10 and L=20. Finite-size errors are within
thestatistical error bars.
Upto values oft/U ∼0.031no size effects weredetected
Lehman expansion and extrapolation to the τ → ±∞
within the error bars. [Here and throughout the paper
limit one readily finds that
errorbarsareoftwostandarddeviations]. Inthe critical
region the finite-size effects were eliminated using stan-
G(p,τ) → Z+e−ǫ+(p)τ , τ → +∞, (3) dard scaling techniques (see below). The dashed lines
(cid:26)Z−eǫ−(p)τ , τ → −∞. are the prediction of Ref. [15] based on the third-order
expansion in t/U. It becomes inaccurate quite far from
The two limits describe single-particle/hole excitations the tip when the insulating gap is about ∼ 6t. On the
in the MI phase. Here Z± and ǫ± are the particle/hole other hand, the value of the tip position extrapolated
spectral weight (or Z-factors) and energy, respectively. to the infinite order is right on target, within the error
In the grand canonical ensemble, excitation energies are bar of order 3t for the chemical potential and on-site
measured relative to the chemical potential. With this repulsion. In all simulations (performed at t/T = 40)
in mind, calculating the phase diagram of the system is the finite-temperature effects are negligible—the system
rather straightforward. At a fixed number of particles is essentially in its ground state.
N = L3 and t/U ratio one determines chemical poten- To eliminate finite-size effects in the critical region
tials µ± for which the energy gap for creating the parti- and pinpoint the position of the lobe tip, we employed
cle/hole excitation with p = 0 vanishes. The insulating standard scaling techniques based on the universality
gap is given then by ∆ = µ+ −µ−. For high precision considerations.
simulations of the gap one has to choose the value of µ First, let us briefly review the universal proper-
veryclosetoµ± andconsiderfinite,butzeroforallprac- ties of the insulator-to-superfluid transition (see Ref. [2]
ticalpurposes,valueoftemperaturesothatthefollowing for more details). There exist two types of transi-
two conditions are satisfied: tions: the “generic” transition, when the phase bound-
ary is crossed at fixed t/U, and a special transition
|µ−µ |≪t, |µ−µ |≫T . (4)
± ±
at fixed integer density, when the SF-MI boundary is
This is exactly how we proceed. By plotting ln[G(p,τ)] crossedat fixedµ/U. The generictransitionis drivenby
vs. τ we deduceǫ (p) fromthe exponentialdecayofthe the addition/subtraction of a small number of particles,
±
Greenfunction. AtypicalexampleisshowninFig.2. We and is fully characterized by the physics of the weakly-
use the valuesofthe hoppingamplitudet andthe lattice interactingBosegasformedbythesmallincommensurate
constant a as units of energy and distance, respectively. densitycomponentn−n ,wheren isthenearestinteger
0 0
InFig.3 we presentaccurateresultsfor the bound- to n. In particular, if δ is the deviation from the generic
aries of the first Mott insulator lobe. We have done cal- critical point in the chemical potential or t/U ratio then
culations for systems with linear sizes L=5, 10, 15, 20. |n−n |∼δ and T (δ)∼δ2/3 in the SF phase.
0 c
5
The special transition at the tip of the lobe hap-
∆L/t
pens at fixed integer density. It is driven by delocalizing
quantum fluctuations which for large values of t/U en-
able bosons to overcome the on-site repulsion and hop 20
within the lattice. As explained in Ref. [2], the effective
action for the special transition belongs to the (d+1)-
dimensional XY universality class which implies emer- 18
gent relativistic invariance (rotational invariance in the
imaginary-time–space,whichisequivalenttotheLorentz
invariance in real-time–space), and, in particular, an
16
emergent particle-hole symmetry. The upper critical di-
(t/U) =0.03408(2)
mension for this transition is (d + 1) = 4, so that for c
d > 3 the critical exponents for the order parameter, β,
14
andthecorrelationlength,ν,areofmean-fieldcharacter: 0.03400 0.03405 0.03410
t/U
β =ν =1/2 (with logarithmic corrections for d=3). In
thisstudy,wewerenotabletoresolvelogarithmicrenor-
FIG.4: (Coloronline). Finitesizescalingoftheenergygapat
malizations for realistic 3D systems and proceed below
the tip of the lobe. ∆L/t vs. t/U for system size L=5 (solid
with the analysis which assumes mean-fieldscaling laws. squares), L=10 (open circles), L=15 (solid circles), L=20
Denotingthedistancefromthecriticalpointasγ = (open squares). Lines represent linear fits used to extract
[(t/U)c−t/U],forasystemoflinearsizeLonecanwrite thecritical point.
∆(γ,L) = ξ−1f(ξ/L) = L−1g(γL2), (5)
-2
where∆istheparticle-holeexcitationgap,ξisthecorre- 10 3·. 4t/0U8
lationlength,andf(x)andg(x)aretheuniversalscaling
functions. In the last expression we have used the rela-
tion ξ ∝ γ−1/2. At the critical point, the product L∆ 3.406
does not depend on the system size. Therefore, by plot-
ting L∆ as a function of t/U one determines the critical
3.404
point from the intersection of curves referring to differ-
ent values of L, as shown in Fig. 4. This analysis yields
(Fig. 5 explains how finite-size effect in the position of 3.402
the crossingpoint originatingfromcorrectionsto scaling
was eliminated)
3.400
(t/U) = 0.03408(2) (n = 1). (6)
c 0.0 0.2 0.4 0.6 0.8 1.0 -2
(10·L)
Our final results are summarized in Fig. 3. We find
that the size of the critical region where 4D XY scal-
FIG. 5: Extrapolation to the thermodynamic limit. We
ing laws apply is narrow and restricted to small gaps of
showtheintersections (triangles) ofthecurves(L=5,L=10),
the order of ∆ ≤ t (inside the vertical error bar on the (L=10, L=15), (L=15, L=20), vs. L−2 . The fit (solid line)
strong-coupling expansion result in Fig. 3). It appears max
yields (t/U) =0.03408(2).
c
that resolving this limit experimentally would be very
demanding.
To perform analytic estimates of the MI state en- is a reasonable approximation for all values of k in the
ergy (and entropy) at low temperature T ≪ ∆ one has Brillouin zone. Note that, if one is to use the local den-
to know effective masses of particle and hole excitations, sity approximation (LDA) for the energy/entropy esti-
m±. For example, the particle/hole contributions to en- mates of trapped systems, then calculations have to be
ergy in the grand canonical ensemble are given by the performed in the grand canonical ensemble.
sums Todetermine effectivemasseswecomputed G(p,τ)
L 3 intheinsulatingstateanddeducedǫ±(p)forseverallow-
E± = ǫ±(k)nǫ ≈ dkǫ±(k)e−ǫ±(k)/T , est momenta from the exponential decay of the Green
(cid:18)2π(cid:19) Z
Xk function on large time scales. Dispersion laws were then
(7)
fitted by a parabola, with the exception for the diagram
where n is the Bose function and
ǫ tip, where the dispersion relation is relativistic. The re-
ǫ±(k) ≈ ±(µ±−µ)+k2/2m± (8) sult for m± is shown in Fig. 6. When t/U → 0 one can
calculate effective masses perturbatively in t/U to get
For large gaps the tight binding approximation
tm = 0.25−3t/U , tm = 0.5−12t/U . (10)
1−cosk + −
ǫ (k) ≈ ±(µ −µ) + α (9)
± ±
α=Xx,y,z m± Clearly, our data are converging to the analytical result
6
0.5
t·m E/Nt
±
10 400 nexc
0.4
300
8 200
0.3
100
6
0
0 10 20 30T/t
0.2
4
0.1 2
0
0.0
0.000 0.010 0.020 0.030 t/U 0 5 10 15 20 25 30 T3/5t
FIG. 6: (Color online). Effective mass for hole (solid cir-
cles) and particle (open circles) excitations as a function of
t/U. The exact results at t/U = 0 are m = 0.25/t and FIG. 7: (Color online). Energy per particle at t/U = 0.005,
+
m− =0.5/t. Bydashedlinesweshowthelowestorderint/U unity filling factor and linear system size L = 20. Solid cir-
correction to the effective masses. Close to the critical point cles are prime data (error bars within symbol size), the solid
thetwocurvesoverlap,directlydemonstratingtheemergence line is the analytical prediction from Eq. 7, where, at each
of the particle-hole symmetry. At t/U = 0.034, the sound temperature,thechemicalpotentialhasbeenfixedbyimpos-
velocity is c/t=6.3±0.4. ing equal number of particle and hole excitations. The inset
shows the total numberof excitations present in thesystem.
as t/U → 0. On approach to the critical point the ef-
A we probe the limits of applicability of semi-analytic
fective mass curves become identical for particles and
predictions in the Mott state. In Subsection B we cal-
holes indicating that there is an emergent particle-hole
culate the entropy of the Bose-Hubbardmodel, compare
symmetry at the diagram tip. In agreement with the
it to the initial entropy (i.e. before the optical lattice
theoretical prediction, the data taken at t/U = 0.034
is turned on) and estimate the final temperature. We
are fitted best with the relativistic dispersion relation
ε(p) = c m2c2+p2, where c is the sound velocity and consider both a homogeneous system in the MI and SF
the effecptive∗mass is defined as m∗ = ∆/2c2. At this states and a system of N ∼3·104 particles in a trap.
value of t/U we have found that c/t = 6.3±0.4 and
tm = 0.010±0.004.
∗
A. Comparison with low T semi-analytic
predictions
III. FINITE TEMPERATURE ANALYSIS
Awayfromthetipofthelobe,intheMIstate,semi-
Controllingthetemperatureisanimportantexperi- analytic predictions are reliable provided the tempera-
mentalissue,crucialformanyapplicationsofcoldatomic ture is low enough, i.e. T ≪ ∆. In other words, there
systems and studies of quantum phase transitions. In exist a range of temperatures , defined as T .const·∆,
this section,we discuss thermodynamicproperties ofthe where the quasi-particle excitations can be successfully
Bose-Hubbard model. We present data for energy, spe- describedasanon-interactingclassicalgas(seeEq. (7)).
cific heat and entropy, for some specific cases. In partic- The value of the constant depends on ∆, as the two fol-
ular, we focus on the most important hni ≈ 1 situation. lowing examples demonstrate.
Ourdatacanbeusedintwoways: (i)tounderstandlim- Let us first consider larger gaps, e.g. ∆∼200t, for
its of applicability of the semi-analytic approach (with which we have found that the low temperature analytic
calculatedeffectiveparameters)discussedabove,and(ii) predictions reproduce numerical data very well. In Fig.
to have reference first-principle curves for more refined 7 we plot the energy per particle in the low temperature
numericalanalysis. Unfortunately, a directsimulationof regime for ∆=181.6t. The analytic prediction from Eq.
a realistic case in the trap, i.e. with similar number of (7), where, at any given temperature, the chemical po-
particles asin experiments,is still a challengingproblem tentialhasbeenchosenbysettingequalthetotalnumber
thoughsimulationsofabout 105 particlesor moreatlow ofparticleandholeexcitations (asit isdone for intrinsic
temperature seem feasible in near future. semiconductors), is reliable up to temperatures T .35t,
The results are organized as follows: in Subsection i.e. T . 0.15∆. In the inset we plot the average num-
7
1.0
28 E cV S(kB)
Nt 0.6 N N
0.8
26
0.6
0.4
24
0.4
22 0.2
0.2
20
0.0 0.0
0 5 10 15 T/t 20 0 5 10 15 T/t 20 0 5 10 15 T/t20
FIG. 8: (Color online). Energy (left), specific heat (center) and entropy (right) per particle at t/U =0.025 (MI ground state)
and unity filling factor. On the left, solid circles refer to prime data (error bars within symbol size). the data were taken for
linearsystemsizesL=10andL=20. Withintheerrorbars,wearenotabletoresolveanyfinitesizeeffect. Solidlinesinallplots
are obtained from spline-interpolated datafor energy, with subsequent analytic differentiation/integration of theinterpolation
curve.
ber of particle-hole excitations. This number increases the chemical potential
rapidly with temperature justifying the grand-canonical
calculationforthe quasi-particlegas(atfixedtotalnum- µ(ieff) = µ − V(i). (11)
ber of particles). The quasi-particle number density is
In strongly interacting regimes with a short healing and
∼ 5% at T ∼ 35t. Apparently, for higher temperatures
correlation length, the LDA approach can be easily jus-
the ideal gas picture is no longer valid as it crosses over
tified in most cases (critical regions of phase transitions
tothatofthe stronglycorrelatednormalliquid. We con-
excluded). In Ref. [21], the authors directly compare
clude that for large enough gaps and T . 0.15∆, one
simulation results for 1D and 2D harmonically trapped
can rely on low temperature analytical predictions to do
systems with LDA predictions based on known homoge-
thermometry.
neous system phase diagram. As expected, the density
For smaller gaps,instead, we do not find any inter- profilesdifferonlyattheMI-SFinterface,andwebelieve
esting region(i.e. where temperature effects are visible), that the same will be true for the 3D case which is more
for which the classical description is valid. In Fig. 8 we “mean-field-like”.
showtheenergyperparticleasafunctionoftemperature Whenthesemi-analyticpredictionsarereliable(see
for t/U = 0.025 (the groundstate is MI with the energy Fig. 7), one can use numerical results for the effective
gap ∆ = 18.35t). To get the specific heat and entropy, massesand gaps to calculate the entropy of the homoge-
wefirstusesplineinterpolationoftheenergydatapoints neousquasi-particlegaswiththetight-bindingdispersion
toobtainasmoothcurveE(T). Thespecificheatisthen relation. The entropy is given by:
obtained by differentiating the spline. The maximum in
the specific heat is reached when temperature is about S = − V d3k∂[Ω+(k)+Ω−(k)] , (12)
half the energy gap. The entropy has been calculated (2π)3 Z ∂T
by numerical integration of c /T. In order to see any
V where
temperature effect one has to go as high as T ∼ 2.5t; at
thesetemperaturestheclassicaldescriptionisalreadyno Ω (k) = T ln 1−expǫ±(k) . (13)
longer applicable and one has to rely on numerical data ± (cid:18) T (cid:19)
to do thermometry.
As anexample,considera uniformweaklyinteract-
ing Bose gas (WIBG) of 87Rb with the gas parameter
na3 ∼ 10−6, which is loaded into an optical lattice with
s
λ = 840nm and t/U = 0.005. At low enough temper-
B. Loading the optical lattice: estimate of T(fin)
ature, T . 0.3T , one can calculate the initial (prior to
from entropy matching c
imposingthelattice)entropyofthesystemusingtheBo-
goliubovspectrum. InFig.9weplottheentropyperunit
The standard approach to convert results obtained volume before and after the optical lattice is imposed.
forahomogeneoussystemintopredictionsforsystemsin The bottom x axis is temperature in units of the critical
external fields is the so-called local density approxima- temperatureofthe WIBG,the topx axis is temperature
tion (LDA), which is actually a local chemical potential in units of t, the hopping matrix element. The dashed
approximation when the density at the site i is identi- line is the entropy of the WIBG, the solid and dashed-
fied with the density of the homogeneous system with dotted lines represent the entropy of the Bose-Hubbard
8
0 5 10 15 20 T/t wheret is the reducedtemperature T/Tc. At verylowT
(T < µ), Eq. (14) misses the contribution coming from
-5 collective excitations. We are interested in initial tem-
2.5 S/V·10 2.5
peratures T ∼ 0.2−0.3T , which are feasible in current
c
experiments [25]. Starting from Eq. (14), we calculate
2.0 2.0 theentropyoftheBECinitiallypreparedinthemagnetic
trap. After the optical lattice is adiabatically turned on,
-50 1.5 1.5 x10tphoetemntaiganleotficthpeotheonmtioaglepnreoovuidsessyastsecman[soeveeErqth.e(1c1h)e]m. ical
x1 -5 Adirectcomparisonwithexperimentsatfixednum-
1.0 1.0 ber of particles would require to calculateµ(T) from the
normalization condition. At low temperatures, one ex-
pects the dependence of the chemical potential on tem-
0.5 0.5
perature to be weak (this will be confirmed by direct
simulations, see below). For simplicity, we fix µ at a
0.0 0.0 value corresponding to N = 303 trapped atoms in the
0.00 0.05 0.10 0.15 0.20 0.25
T/T
C firstMott lobe, at zero temperature. Fromthis point we
proceed in two directions. On one hand, we analytically
calculatethelowtemperaturecontributiontoenergyand
FIG. 9: (Color online). Entropy per particle at t/U =0.005,
entropy arising from particle and hole excitations in the
unityfillingfactorandlinearsystemsizeL=20. Thedashed
line is the entropy of the uniform WIBG. The solid line is trapped MI state. On the other hand, we directly sim-
the result of analytical derivation/integration of numerical ulate the thermodynamics of the inhomogeneous system
data for energy and the dashed-dotted line (barely visible) at a fixed chemical potential. The results are shown in
is the analytical prediction, Eq. (12), where, at each temper- Fig. 10, where we plot the energy per particle, counted
ature, thechemical potentialhas been fixedbythecondition from the ground state. The solid circles are data from
ofhavingequalnumberofparticleandholeexcitations. Ifthe the simulations (error bars are plotted), the solid line is
system was initially cooled down to Tin = 0.25T , the final
c the (analytically calculated) contribution of the particle
temperature is Tfin =22t and nearly a hundred of thermally
and hole excitations. The inset shows the low tempera-
activatedparticle-holeexcitationsarepresentinthefinalstate
ture region. A large mismatch between the two results
(see inset in Fig. 7).
indicates that the main contribution to energy is given
by the liquid at the perimeter of the trap. At zero tem-
perature, there are about N ∼ 29000 particles in the
model calculated starting from the numerical results of
trap, 7% of which are not in the MI state (recall that
Fig.7andanalyticalpredictionsofEq.(12),respectively.
µ has been determined by placing 303 particles in the
If the system was preparedat T ∼0.25T , the final tem-
c
MI state). Simulation results show that, in the range
perature would be T ∼ 22t. Fig. 7 shows that at this
of temperatures considered, the total number of parti-
temperaturethesystemisquitefarfromitsgroundstate
cles increases by 0.7% at most, which confirms the weak
and the number density of thermally activated particle-
temperature dependence of the chemical potential. In
hole excitations is ∼1%. The circumstances of this kind
addition, for T = 8t, we performed a simulation with N
becomecrucialifoneistousethesysteminquantumin-
fixed at the groundstate value. The energies per particle
formationprocessing. This exampleis alsoillustrativeof
in the canonical and grand canonical simulations differ
how numerical data can be used to suggest appropriate
by 0.3% only, therefore we proceed calculating the en-
initial conditions.
tropyin the canonicalensemble andcompare it with the
Now we turn to a more realistic case of confined
initial entropy of the system, at a fixed particle number.
system and use LDA to convert results for the uniform
system into predictions for the inhomogeneous one. In Weareinapositiontoaddressthequestionofwhat
what follows, we consider a gas of N ∼ 3 · 104 87Rb is the final temperature of the system after the optical
atoms, magnetically trapped in isotropic harmonic po- lattice is turned on and the final state is MI with the
tential of frequency 2π60 Hz. Experiments with such exceptionofasmallshellatthetrapperimeter. InFig.11
number of particles were recently performed [9]. With we plot the entropy of the trapped WIBG with (solid
this geometry, the parameter η = 1.57(N1/6a /a )2/5 line) and without (dashed line) the optical potential. If
s ho
(seeRef.[22])is∼0.33,whichisatypicalvalueincurrent theinitialsystemiscooleddowntotemperaturesT(in) ∼
experiments. For temperaturesin the rangeµ<T <Tc, 0.25Tc, see, e.g., Ref. [25], the final temperature will be
where T is the critical temperature of the harmonically T(fin) =(2.35±0.30)t.
c
trapped ideal gas, one can accurately calculate energy Withtheinitialconditionsconsideredinthisexam-
using the Hartree-Fock [23] mean field approach [24]: ple, what is the final state of the liquid at the perimeter
of the trap? Before answering this question, we would
E 3ζ(4) 1 like to recall that, along the MI-SF transition lines, the
= t4+ η(1−t)2/5(5+16t3), (14)
NT ζ(3) 7 critical temperature for the normal-to-superfluid transi-
c
9
0.0 0.5 1.0 1.5 2.0 2.5 T/t
E/Nt
S(k )/N
0.15 B
1.5
0.10 0.15 0.15
0.05
1.0 0.00
0.5 1.0 1.5 2.0 2.5 0.10 0.10
0.5 0.05 0.05
0.0 0.00 0.00
0 2 4 6 8 T/t 0.05 0.10 0.15 0.20 0.25 T/T
C
FIG. 10: (Color online). Energy per particle at t/U =0.005,
FIG. 11: (Color online). Entropy per particle for the same
µ = 116.5t and trap frequency ω = 2π60Hz. Solid circles
system as in Fig.10. The dashed line is the entropy of the
are numerical data, the solid line is the energy of particle
trapped WIBG, before the optical lattice is loaded; the solid
and hole excitations in MI deduced from Eq. (7). The inset
line is the final entropy. If the system was initially cooled
shows azoom of thelow temperaturerange. AtT ∼2.3t the
down to T(in) = 0.25T , the final temperature is T(fin) =
contribution of MI excitations to energy is ∼10% only. c
(2.35±0.3)t and the liquid at the perimeter of the trap is
normal (see text).
tionis zero. The transitiontemperature increasesas one
moves away from the border of the Mott lobes (lower- of the lowest mode is finite: E = cp , with c ≈ 6t
min min
ing µ at fixed t/U in our case) and the quasi-particle and p = 2π/L). The specific heat and entropy are
min
density increases until it reaches its maximum at about calculated as described for Fig. 8. We were not able to
n ≈ 1/2 and then decreases. The maximum Tc can resolvetheSF-NLtransitiontemperaturefromthissetof
be estimated from the ideal Bose gas relation T(0) = dataalone: numericaldatacorrespondingtosystemsizes
c
3.313n2/3/m∗ = 4.174t, with n = 0.5 and m∗ = 1/2t, L=10 and L=20 overlapwithin error bars and we did
but interaction effects are likely to reduce this value. not see any feature at the critical temperature. This is
We have performed simulations at half filling factor and not surprising since the specific heat critical exponent α
fixed t/U =0.005,and found the criticaltemperature to is very small, and it is thus very difficult to resolve the
be T (n = 1/2) = 2.09(1)t. As a consequence, for the singularcontributionandfinite sizeeffects inenergyand
c
chosen initial conditions, we can conclude that the final specific heat. For a system of linear size L, the singular
state of the liquid at the perimeter of the trap is normal part of the specific heat, c(s), can be written as
V
and it gives the main contribution to the entropy. For
such low final temperature, the contribution to the en- c(s)(t,L) = ξα/νf (ξ/L) = Lα/νg (tL1/ν), (15)
tropy per particle due to thermally activated excitations V c c
in the MI state is only 10%. The largest chemical po-
wheret=(T−T )/T ,α≈−0.01,ν =(2−α)/3,ξ isthe
c c
tential is, in fact, deep in the first Mott lobe, and the
correlation length and f (x) and g (x) are the universal
c c
energy required to introduce an extra particle or hole is
scaling functions. At the critical point, finite size effects
much larger than T. Most excitations are located at the
for the two system sizes considered are ∼ 1%, within
perimeter ofthe Mott statein anarrowshellofradiusR
error bars.
and width ∼0.05R.
The critical temperature was extracted from data
Retrievingthesameinformationforexperimentsus- for the superfluid stiffness. The scaling of the superfluid
ing a larger number of particles, e.g. 105 ÷106, by di- stiffness at the critical temperature is n ∝ |t|ν. This
s
rect simulation is still computationally challenging. In allowsonetoaccuratelyestimatethecriticaltemperature
order to use LDA, one should study the uniform sys- from
tem, scanning through the chemical potential. As our
last example, we consider a uniform system which is in n (t,L)=ξ−1f (ξ/L) = L−1g (tL1/ν), (16)
s s s
the correlated SF ground state. Fig. 12 shows data for
E, c , and S for t/U = 0.0385 and unity filling factor, by plotting n L vs. T as shown in Fig. 13. From the
V s
close to the MI-SF transition. The system stays in its data taken for system sizes L = 5, L = 10, L = 20, we
ground state for T ≪ 2t; in finite systems, the energy estimate the critical temperature to be T = 3.25(1)t.
c
10
20 NEt 0.6 cNV 1.0 S(NkB)
18 0.8
0.4
0.6
16
0.4
0.2
14
0.2
12 0.0 0.0
0 5 10 15 T/t 20 0 5 10 15 T/t20 0 5 10 15 T/t 20
FIG.12: (Coloronline). Energy(left),specificheat(center)andentropy(right)perparticleatt/U =0.0385(SFgroundstate)
andunityfillingfactor. Ontheleft,solidcirclesrefertoprimedata(errorbarswithinsymbolsize). Dataweretakenforsystem
size L=10 and L=20. Within error bars, we are not able to resolve any finite size effect. Solid lines in all plots are obtained
from spline-interpolated data for energy, with subsequent analytic differentiation/integration of theinterpolation curve.
region in the vicinity of the diagram tip where universal
3
n ·L propertiesof the relativisticeffective theorycan be seen.
S
Comparison with the strong-coupling expansion shows
that the latter workswell only for sufficiently largeinsu-
lating gaps ∆>6t outside of the fluctuation region. We
2 havestudied the effective masses of particle andhole ex-
citationsalongtheMI-SFboundary. Ourresultsdirectly
demonstrate the emergence of the particle-hole symme-
try at the diagram tip, and provide base for accurate
theoretical estimates of the MI thermodynamics at low
1 and intermediate temperatures.
We have studied thermodynamic properties of the
(T/t) = 3.25(1) superfluid and insulating phases at fixed particles num-
C
berfortheuniformcase. Thesedatacanbeusedtomake
predictions for the inhomogeneous system using the lo-
0
2.8 3.0 3.2 3.4 3.6 3.8 4.0 T/t caldensityapproximation. We haveshownthatforlarge
enough gaps the low temperature analytical predictions
agree with numerical data. By entropy matching, we
FIG. 13: (Color online). Finite size scaling of the super-
have calculated the final temperature of the system (af-
fluid stiffness. nsL vs. T/t for system size L=5 (circles),
L=10 (squares), L=20 (triangles). We estimate the critical tertheopticallattice isadiabaticallyloaded),intheuni-
temperaturetobeT =3.25(1)t,alreadynearlyhalfthenon- form and magnetically trapped system, at t/U = 0.005.
c
interacting gas value. We have performed direct simulations of a trapped sys-
tem, using typical experimental values for the magnetic
potential and number of particles. For the initial con-
At this temperature, the entropy per unit particle is ditions considered, we found the final temperature and
∼ 0.195, or, translating to entropy density in physical demonstrated that the main contribution to the entropy
units,3.6·10−5JK−1m−3,whichcorrespondstoaninitial comes from the liquid at the perimeter of the trap. We
temperature ∼0.35T(in). Thereforeit seems plausible to have calculated the normal-to-superfluid transition tem-
c
reach T experimentally. perature at the half filing and concluded that the liquid
c
at the perimeter is normal.
IV. CONCLUSIONS
V. ACKNOWLEDGEMENTS
We have performed quantum Monte Carlo sim-
ulations of the three dimensional homogeneous Bose-
Hubbard model. We were able to establish the phase WearegratefultoMatthiasTroyerforusefuldiscus-
diagram of the MI-SF transition with the record accu- sions. This work was supported by the National Science
racy ∼ 0.1% and determine the size of the fluctuation Foundation Grant No. PHY-0426881.