Bose-glass to Superfluid transition in the three-dimensional Bose-Hubbard Model

0.27 MB·English

by  Peter Hitchcock

Bose-glass to Superfluid transition in the three-dimensional Bose-Hubbard Model
Peter Hitchcock∗ and Erik S. Sørensen†
Department of Physics and Astronomy, McMaster University, Hamilton, ON, L8S 4M1 Canada
(Dated: February 5, 2008)
WepresentaMonteCarlostudyoftheBose-glasstosuperfluidtransitioninthethree-dimensional
Bose-Hubbard model. Simulations are performed on the classical (3 + 1) dimensional link-current
6 representationusingthegeometricalwormalgorithm. Finite-sizescalinganalysis(onlatticesaslarge
0 as 16x16x16x512 sites) of the superfluid stiffness and the compressibility is consistent with a value
0 of the dynamical critical exponent z = 3, in agreement with existing scaling and renormalization
2 group arguments that z = d. We find also a value of ν = 0.70(12) for the correlation length
exponent, satisfying the relation ν ≥ 2/d. However, a detailed study of the correlation functions,
n
a C(r,τ), at the quantum critical point are not consistent with this value of z. We speculate that
J thisdiscrepancy could beduetothefact that thecorrelation functions havenot reached their true
asymptoticbehaviorbecauseoftherelativelysmallspatialextentofthelatticesusedinthepresent
0
study.
1

] PACSnumbers: 67.40.Db,67.90.+z
l
e
-
r I. INTRODUCTION fluidscalingarguments7,basedongeneralizedJosephson
t
s relations, then show that:
.
at Quantum phase transitions as they occur in models κ δν(d−z), (1)
m comprised of bosons have been the focus of considerable ∼
interest lately. Most notably, experiments by Greiner et
- where δ describes the distance to the critical point in
d al.1 haveobservedthetransitionbetweenaMottinsulat-
terms of the control parameter. Since both the BG and
n ing (MI) phase and a superfluid phase in an optical lat-
SF phases are compressible it then seems natural to ex-
o tice loaded with 87Rb atoms. The observed phase tran-
c pect the compressibility to be finite also at the critical
sition between the superfluid and the insulating phase
[ point and one must then conclude that
is thought to share the universal properties of a variety

1 ofphysicalsystems,including He4 inporoussubstrates2,
z =d. (2)
v Josephson-junctionarrays3aswellasthinsuperconduct-
2
8 ing films4,5,6. While disorder can be neglected in the It can in fact be shown that assuming the compressibil-
1 experiments by Greiner et al.1, it clearly plays a central ity either diverges or tends to zero at the BG-SF critical
1 roleinseveraloftheexperimentsjustmentioned. Asout- pointleadstoimplausiblebehavior. InRef.7therelation
0 linedin the seminalpaperby Fisher et al.7 the statics as z = d was therefore argued to hold not only below the
6 well as the dynamics of the quantum phase transitions upper critical dimension but for any dimension d 1.
0 occurring in the Bose-Hubbardmodel are strongly influ- This result implies that there is not an upper critica≥l di-
/
t encedby disorderand anew insulating glassyphase,the mension for this transition in any conventionalsense. In
a
m Bose-glass phase, should occur. The transition of inter- d = 1 analytical work9 strongly supports the conclusion
est in this case is directly between the superfluid (SF) that z = d = 1 and numerical work11,12,18 in d = 2 also
d- and the Bose-glass (BG). Recent experiments showing find z = d = 2. Long-range interactions are likely to
n that this disordered transition can also be studied using yield a different z7, but we shall not be concerned with
o optical lattices8 has therefore created considerable ex- that case here. Based on Dorogovtsev’s19 double-ǫ ex-
c citement. Previous theoretical studies9,10,11,12,13,14 have pansion, an alternative scenario has been proposed20,21
v: mainly focused on the transition as it occurs in one or in which the critical exponents jump discontinuously to
i two dimensions and very recent theoretical works15,16, their mean field values at the critical dimension d = 4.
X c
addressing directly the experimental situation relevant In this scenario, for d > d two stable fixed points are
c
ar for optical lattices, have also mainly analyzed this case. present: the gaussian fixed point is stable at weak dis-
In light of the experiments by Lye et al.8,17, we consider order with the random fixed point remaining stable at
in the present study the transition as it occurs in three strong disorder. Below the critical dimension, d < d ,
c
dimensionsinthepresenceofshort-rangeinteractions. In onlytherandomfixedpointisstable. Itshouldbe noted
particularwedeterminethedynamicalcriticalexponent, thatinordertoobtainsensibleanswersfromthedouble-ǫ
z, characterizing this transition in d=3. expansionanultravioletfrequencycutoffω hastobein-
Λ
Onehall-markoftheMottinsulatingphaseisthatitis troduced20—a procedure which seems difficult to justify.
incompressible. In contrast to this, both the superfluid This was later remedied upon21 still within the frame-
and the Bose-glass phase are compressible and have a work of the double-ǫ expansion. It is also possible to
finitenon-zerocompressibility,κ=0. Closetothequan- questionthe validity of anyǫ-expansionapproachon the
6
tum criticalpoint between the Bose-glassandthe super- groundsthattheexistenceofanuppercriticaldimension
2
for the disordered transition is implicitly assumed. One The paper is outlined as follows: the Bose-Hubbard
mightthenaskifboundsontheuppercriticaldimensions modelisintroducedbelowinfurtherdetail,includingthe
exist. Indeed, using the exact inequality22 mappingtoa(d+1)dimensionalclassicalmodelonwhich
we perform our study. Section IIIA outlines the scaling
ν 2/d, (3) theories necessaryto extractthe critical exponents. Sec-
≥
tion IIIB details the numerical techniques and discusses
valid for stable fixed points in the presence of correlated
the difficulties we encountered in properly equilibrating
disorder,oneimmediatelyseesthattherequirementthat
oursimulations. Finally,wepresentourresultsinsection
ν =1/2attheuppercriticaldimension,d ,yieldsd 4.
c c ≥ IV.
Since the existence of an upper critical dimension is de-
batable it is then natural to instead focus on the lower-
critical dimension, which is well established. Indeed, it
II. MODEL
is also possible to develop an expansion away from the
lower-critical dimension (d = 1)23,24. Using this ap-
l WebeginwiththeBose-HubbardHamiltonian,includ-
proach the dynamical critical exponent is exactly z = d
ing on-site disorder in the chemical potential:7
and the correlation length exponent ν can be calculated
to second order in √d 1 yielding good agreement with U t
experimental and num−erical results in d = 2. The onset HBH =X(cid:16)2nˆr(nˆr−1)−µrnˆr(cid:17)−2 X(Φˆ†rΦˆr′+H.c.).
ofmean fieldbehavior,if any,is therefore atbest uncon- r hr,r′i
ventional in this model and still a matter of debate. In (4)
particular, it would be valuable to know if the relation Here Φˆ†r and Φˆr are boson creationand annihilation op-
z =dcontinues toholdin dimensionshigher thand=2. eratorsatsiter,andnˆr =Φˆ†rΦˆr isthenumberoperator.
In the present work, we present a Monte Carlo study The on-site repulsive interaction U localizes the bosons
of the three-dimensional Bose-glass to superfluid transi- and competes with the delocalizing effects of the tunnel-
tion at strong disorder. Our focus is reliable estimates ing coefficient t. The random chemical potential µr is
of the dynamical critical exponent z in order to test distributed uniformly on [µ ∆,µ+∆]. As noted by
−
the relation z = d = 3. As outlined above, much of Damski et al.16, the introduction of a random potential
the numerical work to date has focused on the one- and in anoptical lattice will also generaterandomness in the
two-dimensionalcases,whicharelessdemandingcompu- hopping(tunneling)term,t. However,this typeofdisor-
tationally. However, recently a very efficient geometric der can be ignored16. The phase diagram7 for the pure
wormalgorithm18,25 has been developedfor the study of modelconsistsofasuperfluidphaseathight/U whichis
bosonicphasetransitionsmakingsignificantlylargersys- unstable to a series of Mott-insulating regimes centered
tem sizes available. Even though the geometrical worm at commensurate densities at low t/U. In the presence
algorithmwe usehas provento be veryefficient fordeal- of disorder, the Bose-glass phase stabilizes between the
ingwithlargelatticesintwodimensions(alsointhepres- insulating and superfluid phases.
ence of disorder),we were only able to study a relatively Following standard methods11,31, by integrating out
limited number of system sizes in d = 3 that were large amplitude fluctuations of the Bose field to second order,
enough to avoid finite-size effects but small enough to we transform the Bose-Hubbard Hamiltonian into an ef-
properly equilibrate with a feasible amount of computa- fective classical Hamiltonian in (d + 1) dimensions well
tional effort. In the present study therefore, we focus suited for Monte Carlo study:
exclusively on the three dimensional case. We leave for
1 1
afubtluerientsetruedsytftohretchaeseapopfrdoa=ch4ttohamtewaonufildeldbesuogfgceosntseiddebry- H = K Xh2J2(r,τ)−µrJ(τr,τ)i. (5)
(r,τ)
Weichman and collaborators20,21.
It has been suggested that in the vicinity of commen- TheintegercurrentsJ(r,τ)aredefinedonthebondsofthe
suratevaluesofthechemicalpotentialdisorderisweakly lattice and obey a divergenceless constraint J(r,τ) =0.
irrelevant26,27,28,29. InthiscaseadirectMI-SFtransition Theresultingcurrentloopsareinterpretedin∇thiscontext
could occur even in the presence of weak disorder. How- as world lines of bosons11,31 and represent fluctuations
ever, recent high precision numerical work30 reachedthe from an average, non-zero density. The coupling K acts
opposite conclusion for d = 2, although very large sys- asaneffectiveclassicaltemperatureanddrivesthephase
temsizeswereneededinordertoshowthis. Sinceweare transition: atlowK therepulsiveshortrangeinteraction
mainly interested in the BG-SF transition, we minimize U dominates and the system is insulating, while at high
any cross-over effects by performing all of our calcula- K the hopping t dominates andthe systemis superfluid.
tions at a value of the chemical potential where in the In this transformation amplitude fluctuations of the bo-
absence of any disorder the system is in the superfluid son fields are integrated out to second order. However,
phase for any non-zero hopping and there is no MI-SF thesefluctuationsarenotexpectedtoaffecttheuniversal
transition7. The BG-SF transition we observe is there- details of the transition.
fore induced by the disorder and corresponds in the RG TheadvantageoftheformulationoftheBose-Hubbard
sense to a random fixed point. model in terms of the link-currents, Eq. (5), is that an
3
extremely efficient wormalgorithm18 is available for this latticesizeswhosetemporalsizesscalewiththeexponent
model. This wormalgorithmcanalsobe directed25 even z; that is to work at a fixed aspect ratio:
in the presence of disorder. However, the memory over-
head associated with using a directed algorithm is pro- α=Lτ/Lz. (11)
hibitive for the present study due to the presence of dis-
order andthe high spatialdimension. We havetherefore The quantity:
exclusively used the more standard formulation of the
algorithm18. Ld+z−2ρ(L1/νδ,α= Lτ) (12)
Lz
should then be a universal function of α at K . We can
III. METHOD c
hold the aspect ratio constant by working with systems
of dimension Ld αLz. If an initial estimate for the dy-
A. Scaling ×
namicalcriticalexponentz isavailable,the criticalpoint
can be located by plotting Ld+z−2ρ(L1/νδ,α) versus K
The BG-SF transition is a continuous quantum phase for several different linear system sizes L. If the correct
transition, and is characterized by two diverging length value of z is used, these curves will all intersect at the
scales: a spatial (ξ) and a temporal (ξτ) correlation criticalpoint. Thisunfortunatelyrequiresthattheinitial
length, related by the dynamical critical exponent z: estimateforthevalueofzbemadebeforesimulationsare
run.
ξ ξz (δ−ν)z, (6)
τ ∼ ∼ Alternatively, with an estimate of Kc, the first argu-
mentofEq. 12canbe heldconstant. CurvesofρLd+z−2
with δ = (K K )/K . These form the basis of the
scaling theory−usedcto acnalyze our data. plotted for different L against the ratio Lτ/Lz should
then collapse for the correct value of z32,33. We use
The twoobservablesofprimaryinterestare the super-
bothapproaches,withthe hopethataconsistentpicture
fluid stiffness ρ and the compressibility κ. The super-
emerges.
fluid stiffness ρ (proportional to the superfluid density)
OnceK hasbeenlocatedbytheabovemethodwecan
is defined in terms of the change in free energy associ- c
proceed to study the behavior of the compressibility at
ated with a twist in the spatialboundary conditions. As
thecriticalpoint. Ifindeedzisequaltod,asdescribedin
for the compressibility,the criticalbehaviourofρ canbe
the introduction, the compressibility should be roughly
derived as a generalization of the Josephson scaling re-
constant as K varies through K , and should not show
lations for the classical transition7,11,31, and scales with c
any dependence on L. In particular κ should neither
the correlation length as
diverge nor go to zero at K .
c
ρ ξ−(d+z−2). (7) We alsoconsider the particle-particlecorrelationfunc-
∼ tion C(r,τ), which is expected7 to decay asymptotically
The compressibility can similarly be found as a re- as
sponse to a twist applied to the temporal boundary con-
ditions, and is found to scale as7,11,31 C(r,τ =0) r−yr, C(r =0,τ) τ−yτ (13)
∼ ∼
κ ξ−(d−z), (8) with exponents given by
∼
leading to the relation z =d previously mentioned. y = d+z 2+η,
r
−
The first step in the study of the critical properties of y = (d+z 2+η)/z. (14)
τ
theBose-Hubbardmodelisaprecisedeterminationofthe −
location of the critical point through finite-size scaling Inadditiontodefiningtheexponentη,withreliableesti-
analysis. The presence of twocorrelationlengths implies matesofthe correlationfunctions this wouldinprinciple
thatfinite-sizescalingfunctionswillhavetwoarguments: provide for an independent calculation of the dynamical
f(L/ξ,Lτ/ξτ). Hence critical exponent z. Fisher et al.7 also derive bounds for
η:
ρ=ξ−(d+z−2)f(L/ξ,L /ξ ), (9)
τ τ
2 (d+z)<η 2 d. (15)
or,by appropriatelyscalingthe argumentsandrequiring − ≤ −
that the stiffness remain finite at the critical point for a The lower bound arises from the requirement that the
system of finite size, correlation functions decay, and the upper bound is ar-
guedfromthescalingofthedensityofstatesintheBose-
1
ρ= ρ¯(L1/νδ,L /Lz). (10) glass and superfluid phases. These can be simply stated
Ld+z−2 τ
as bounds on the exponents y and y :
r τ
This complicatesthe scalinganalysisaswe mustworkin
atwo-dimensionalspace. Thefirstapproachistoworkat 0<yτ 1, 0<yr z. (16)
≤ ≤
4
30
t
3x1s04 1
25 1x105
3x105
20 13xx110066 0.01 20 4ρL>) 1
4ρ<L>)15 13xx1100780.0001 HHαα,, αβ(0(t)ts) 4ρ<L>) P(<0.01
P( P(
10
10 101 103 t, t 105 107 0.001 0.01 <ρL04.>1 1
s L = 8
5 L = 12
0 0
0.05 0.1 0.15 0.2 0.25 0 0.1 0.2
<ρL4> <ρL4>
FIG.1: Mainpanel: histogramsofhρL4iforasetof1000dis- FIG.2: DistributionsofhρL4iforL=8and12atKc,plotted
orderrealizationsonasystemofsize8x8x8x64atourestimate on linear (main panel) and logarithmic (inset) axes. The av-
of Kc =0.190(1). Curves show the evolution of the distribu- erage of each distributions is close to [hρL4i]av =0.08, which
tion P(hρL4i) as a function of the number of Monte Carlo is significantly higher than the typical values of the distribu-
sweepst performedafterequilibrationoneachdisorderreal- tions, indicating that the broad tail of the distribution must
s
ization. The peak in the distributions at hρL4i = 0 persists be well sampled to obtain an accurate estimate of the true
for many times the relaxation time, as measured by the con- average.
vergence of the Hamming distances, while the broader peak
near hρL4i = 0.025 grows. Inset: Hamming distances calcu-
latedonthesamesetofdisorderrealizations. Hα,β isplotted and is defined in terms of the temporal winding num-
against the total number of sweeps performed, t. Hα,α0 is bers:
plotted against the number of sweeps performed after equili-
bration, t . Open symbols are calculated in terms of spatial L
s κ= τ n2 n n . (19)
currents. Solid symbols are calculated in terms of temporal Ld(cid:2)h τi−h τiαh τiβ(cid:3)av
currents.
Sincethisexpressioncontainsthedisorderaverageofthe
square of a thermal average, the systematic error is re-
The final exponent we calculate is the correlation
duced by calculating this average on two independent
length exponent ν. Taking the derivative of (10) with
lattices α and β with the same disorder realization34,35.
respect to K we see that
We are also interested in the derivative of ρ with re-
dρ spect to the coupling. This can be calculated thermody-
Ld+z−2 =L1/νρ¯′(L1/νδ,L /Lz). (17)
dK τ namically from the stiffness and total energy E:
Plotting this derivative against L at K should yield
c dρ 1
a power law with exponent 1/ν. The crossing data dK = K2(cid:2)hρEi−hρiαhEiβ(cid:3)av. (20)
ρLd+z−2 calculatedfordifferentsystemsizeswithafixed
aspect ratio, α, should also collapse to a single curve This wasfound to producebetter estimates thannumer-
when plotted against L1/νδ. ically differentiating ρ in the two-dimensionalcase36.
Finally, since the construction of each ‘worm’ is es-
sentially equivalent to propagating a boson through the
B. Numerical Method lattice,itispossibletocalculatetheparticle-particlecor-
relationfunctiondirectlyfromthebehaviourofthe algo-
Both the superfluid stiffness and the compressibility rithm. Details can be found in Ref. 25.
can be calculated in terms of the link-current wind- As always, the system must be be run for t Monte
0
ing numbers11,31 nγ = L−γ1Pr,τJrγ,τ in each direction Carlo sweeps at each disorder realization to ensure that
γ = x,y,z,τ (here L = L). The superfluid stiffness equilibriumhasbeenreachedbeforebeginningtosample
x,y,z
is related to a twist in the spatial boundary conditions the generated configurations. To confirm this, two sim-
and so can be calculated in terms of the spatial winding ulations are carried out simultaneously on lattices with
numbers: the same disorder realization but different initial config-
urations (we set all of the currents in eachdirection to a
1
ρ= n2 . (18) differentintegerconstantforαandβ). Theseinitialcon-
Ld−2L (cid:2)h γ=x,y,zi(cid:3)av
τ figurationsarefarfromequilibrium. Itisusefultodefine
We denote thermal averages by and disorder aver- ‘Hamming distances’ between different current configu-
hOi
ages by . Similarly, the compressibility is associ- rations in order to measure the relaxation time of the
(cid:2)O(cid:3)av
ated with a twist in the temporal boundary conditions algorithm (this is done in the spirit of Ref. 31, though
5
0.06 0.003 L = 8
α = 1/32 L = 8 L = 12 α = 1/16
0.25 0.04 L = 10 L = 16
L = 12 0.002 L = 20
0.02 L = 14
0.20 L = 16
00.0.0200 L = 20 0.001
α = 1/64
4 0.015
ρL 0.15 00..001005 α = 1/8 3ρL 0.000
0.10 0.000 0.188 0.190 0.192 0.06 α = 1/8
0.05
0.03
0.00
0.188 0.189 0.190 0.191 0.192
K 0.00
0.196 0.198 0.200 0.202 0.204
K
FIG.3: ThemainpanelshowsρLd+z−2=ρL4 plottedversus
K for lattices of size Lτ = 18L3 assuming d = z = 3. The FIG. 4: Search for a crossing in ρLd+z−2 for z = 2 using
lines are guides to the eye. The crossing gives an estimate of lattices of size L3 ×αLz=2, (see also Fig. 3). Results are
the critical point K = 0.190(1) and is consistent with z =
c shown for two different aspect ratios α = 1/8,1/16. The
3. The two insets show equivalent results with two different
crossings show significant drift between different lattice sizes
aspect ratios α=1/32 and α=1/64.
andbetween thetwo aspect ratios shown. Consequently,z6=
2.
the definitions used here are slightly different). We de-
fine the Hamming distance between the two lattices α ture of the winding numbers which are used to calculate
and β after performing a total of t Monte Carlo sweeps ρ and κ. Since most disorder realizations have a small
on their initial configurations: but finite superfluid stiffness at K , many independent
c
Hαν,=βx,τ(t)= Ld1L X(cid:2)Jαν,(r,τ)(t)−Jβν,(r,τ)(t)(cid:3)2. (21) rceoanlfiizgautrioantiotnosaocfhineγvemaursetlibaeblgeenesetriamteadtef.oFriegaucrhed1issohrodwesr
τ (r,τ) the distribution of ρL4 as a function of the number of
h i
sweeps t performed after equilibration on each realiza-
Similarly, we define the Hamming distance between the s
tion of the disorder. For t t , there are many realiza-
configurationoflattice αatthe sweept0 where webegin tionsforwhichnoconfigursat≈ionrisgeneratedwithafinite
tosamplethe generatedconfigurations(denotedα )and
0 winding number. Only after running for much longer at
the configuration of α after a further t =t t sweeps:
s − 0 each realization are the true features of the distribution
Hν=x,τ(t )= 1 Jν (t +t ) Jν (t ) 2. resolved. Under-samplingtheserealizationsrunstherisk
α,α0 s LdL X(cid:2) α,(r,τ) 0 s − α,(r,τ) 0 (cid:3) of underestimating the disorder average. This issue is
τ (r,τ) discussed at further length in Ref. 35.
In the present study, t is chosen prior to running the The computational demands of equilibration grow
0
simulations and is held constant. quickly with system size. Typical systems we studied,
For the sake of simplicity, we define a Monte Carlo of linear size L = 10 to L = 14, required up to 3 105
×
sweep to be the construction of a single worm; while sweepstorelaxanduptoafurther1 108sweepsforthe
not ideal for comparing its characteristics to other al- distributions of ρL4 to equilibrate.×For L=16, nearly
gorithms, this is sufficient for our discussion here. Since 1 106 sweepswhereriequiredtorelaxandwewereunable
×
the initialconfigurationofα andβ aredifferent, Hγ (t) to run simulations for sufficiently long to see the distri-
α,β
butions equilibrate. Some results for L = 16 are shown
willbelargeatthebeginning ofthe simulation(tsmall).
below as they demonstrate consistent scaling. Figure 2
Conversely, since shortly after t the configuration of α
0
will not have changed substantially, Hγ (t ) will be showsthe equilibrated distributions of ρL4 for two lat-
α,α0 s tice sizes, L=8 and12. They showthehbroaidtailwhich
small. If t has been chosen larger than the relaxation
0
necessitatesrunningonatleast103 disorderrealizations.
time of the algorithm, then for sufficiently large values
Moreover, at least for the system sizes we were able to
of t and t the configurations of α , α, and β will be
s 0
study, the distributions do not narrow as L increases.
independent, equilibrated states, and thus the two dis-
Hence, the system is likely not self-averaging37.
tances should converge. The inset of figure 1 shows
thisapproachtoequilibriumonan8x8x8x64latticewith
t = 3 107 at our estimate of K . The distances con-
0 c
verge a×fter t 105, indicating that the relaxation time IV. RESULTS
≈
for this systemis approximatelythirty thousandsweeps.
A further complication in the calculation of the dis- We begin by a discussion of our results for the scal-
order averages was encountered due to the discrete na- ing of the stiffness ρ. All our results are for the three-
6
0.15
L = 8 0.24 α = 1/8
L = 10
0.12 L = 12
L = 8
0.22 L = 10
L = 12
0.09 L = 14
4
L κ 0.2
ρ
0.06
0.18
0.03
0.16
0
0.01 0.1 1 0.188 0.189 0.19 0.191 0.192
L/L3 K
τ
FIG. 5: ρL4 versusα=L /Lz for system sizes L=8,10,12 FIG. 6: The compressibility, κ shown versus K near Kc =
τ
assuming z = 3 at K = 0.19. Lattices of size L3 × L 0.190(1). Results are shown for L = 8, 10, 12, and 14 with
τ α=1/8. Note the absence of features at K and the lack of
were used. The inflection points and curvature show a very c
dependenceon system size, indicating again that z =d=3.
goodcollapseindicatingthatthedynamicalcriticalexponent
is likely z = 3. The slight vertical drift with increasing L
likelyindicatesthatweareslightlyoffthetrueK =0.190(1),
c
within theboundsof ourerrorbars. Kc = 0.190(1) obtained with z = 3 is correct we can
now check these estimates self-consistently by plotting
ρLd+z−2=4 versus L /Lz=3 for a range of L at K =
τ c
dimensionalcase,hence,thescalingrelationEq.12states 0.190(1). For a detailed explanationof the procedure we
thatLz+1ρ(L1/νδ,α)shoulddisplayacrossingatthecrit- refer to Ref 32. Our results are shown in Fig. 5. The
icalpoint. Asalreadymentioned,wefocusonthevalueof curvesshow ρLd+z−2=4 plotted for systems oflinear size
the chemical potential µr distributed randomly between L = 8, 10, and 12 as a function of the aspect ratio α =
[µ ∆,µ+∆] with µ= 1 and ∆= 1. L /Lz=3. All the curves for different L and L collapse
− 2 2 τ τ
We begin with the ansatz z = d = 3, using lattices of onto a single curve for z = 3. This result is a rather
size L L L αLz=3. In Fig. 3 we show plots of ρL4 strong confirmation that the assumption z =3 indeed is
× × ×
versus K for linear system sizes ranging from L = 8 to correct. If one studies the curves in detail a very slight
16. A clean crossing is observed at K = 0.190(1) for verticaldriftwithincreasingLisnoticeable. Thisdriftis
c
three values of the aspect ratio, α = 1/8 (main panel), likelytheresultofasmalldeviation,withinourerrorbars,
α = 1/32, and α = 1/64 (insets). As can be seen from of the actual critical temperature from our estimate of
Fig. 3, our results for systems of size L=8 do not scale K = 0.190(1). It would have been quite interesting to
c
as well as larger system sizes. This effect is more pro- study systems of linear size L > 12 or systems with an
nounced for L < 8 (not shown). This is most likely due aspectratiolargelyexceedingL /Lz 1,however,given
τ
∼
to corrections to scaling which for these small system ourvalueofthedynamicalcriticalexponentofz =3this
sizes can not be neglected. We also note that an im- is computationally too demanding.
provement in the scaling behaviour of L = 8 with the We now discuss our results for the compressibility, κ.
aspect ratio is visible in Fig. 3. This is presumably due In Fig. 6 the compressibility is shown for a range of K
to the fact that the optimal value for α is given by the around the estimated critical point K = 0.190(1) for
c
relation (ξ/L)z ξ /L 35. severaldifferentsystemsizes,usingthesamelatticesizes
τ τ
∼
In order to test how sensitive we are to the value of z asforinFig. 3. Asexpectedifd=z,thecompressibility
we also tried z = 2, the value for the dynamical critical showsnodependenceonthelinearsizeofthesystemand
exponent quite well established in two dimensions. Our is approximately constant through the transition. Fur-
results for this value of z are shown in Fig. 4 for lattices thermore, at K the compressibility neither diverges nor
c
of dimension L3 αLz=2. For this value of z our results does it tend to zero. This is again a strong confirmation
×
show significant drift in the crossing for different system that z = d = 3. One should note that the results in
sizes and different aspect ratios as can be clearly seen in Fig.6areonlyreallyusefulfordetermining z ifthe criti-
Fig. 4. The apparent crossing between two system sizes calpoint, K is known(since one generallywouldexpect
c
areonlyslightlyhigherthantheclearcrossingseenusing κtobe independentofLfarawayfromthecriticalpoint
z =3in Fig.3 butonly twosystemsizescanbe made to where ξ L). In the absence of any knowledgeof K , κ
c
≪
crossatagivenK. FromtheresultsinFig.4weconclude wouldhave to be calculated for the entire range of phys-
that z =2. ically relevant K. It is therefore interesting to compare
6
AsexplainedinsectionIIIAwecannowholdconstant the results shown in Fig. 6 with the attempt of locating
thefirstargumentofthescalingfunction(12)byworking the critical point assuming z = 2 shown in Fig. 4. The
at our estimate of K . Assuming that our estimate of latter results show curves intersecting pairwise around
c
7
1000 K = 0.188
0.2 K = 0.190
4ρL KK == 00..119901, fixed µ 10-2
0.1
K = 0.192 τ)
K 0 0,
4ρdL/d 100 -0.05 δL01/ν 0.05 C(r, 0), C(10-4 LLL === 81102 α = 1/8, K = 0.19
L = 14
L = 16
L = 20 α = 1/200, K = 0.1905
10
10-6
8 10 12 14 16
L 1 10 r, τ 100
FIG. 7: The derivative of the stiffness, L4dρ/dK, as calcu-
FIG. 8: Spatial and temporal dependence of the particle-
lated directly during the simulation, plotted against the lin- particle correlation function, C(r,τ = 0) and C(r = 0,τ).
ear system size L. Data is shown for a range of couplings
Calculated at K = 0.19 and α = 1/8 for L = 8 to 16. Solid
near Kc. The solid lines indicate power-law fits to the data lines are a fit to C(r,τ = 0) = A(r−yr +(L−r)−yr) and
at K = 0.190, with exponent 1.11 and K = 0.191 with an C(r=0,τ)=A(τ−yτ+(L −τ)−yτ). Eachcorrelation func-
exponent of 1.69. Data is also shown for simulations run in τ
tionwasfittedwiththesamepower,butauniquecoefficient.
the micro-canonical ensemble of disorder at K = 0.190 (see
Fitted values are y = 4.4(4) and y = 1.06(7). Also shown
text). Thedashedlineindicatesapower-lawfittothemicro- arecurvescalculaterdforasystemofτsize203×40(α=1/200)
canonicaldatawithanexponentindistinguishablefromthefit
at K =0.1905. This correlation function showsa clear expo-
tothecanonicaldata. Theinsetshowsthebestcollapseofthe
nential dependence in the spatial direction, though the tem-
scaling curves ρL4 plotted against δL1/ν, with K = 0.1902
c poraldependenceisconsistentwiththesamepowerlawdecay
and ν =0.70. Combining these results we findν =0.70(12).
as for systems of larger aspect ratio.
K 0.2 with the intersections moving downwards with
L.∼If indeed the dynamical critical exponent was z = 2 rorbars of our estimate for Kc). Combining these two
and not z =3 as we show here, then it would seem very estimates of ν we conclude ν =0.70(12). This value of ν
unlikely that the compressibility shown in Fig. 6, just satisfies the “quantum” Harris criterion, ν 2/d, which
≥
below this range of K 0.2, could be so featureless. shouldholdforalldisorder-driventransitions22andcould
One would in that case∼have expected it to behave as be consistent with this inequality being satisfied as an
κ δν(d−z)=ν, resulting in the finite-size scaling form equality. Italsoinverygoodagreementwiththevalueof
κ=∼(1/Ld−z=1)κ¯(L/ξ,L /Lz). SincetheresultsinFig.6 ν =0.69obtainedbyHerbut24 inanexpansioninpowers
τ
areindependentofLtheythereforeclearlyexcludez =2. of √d 1 awayfrom the lower criticaldimension dl =1.
−
We now turn to our results for the correlation length It has been suggested26,38,39 that the inequality ν
≥
exponent ν. Figure 7 shows in the main panel a plot of 2/d can be violated and that in fact ν can be less than
the derivative of the stiffness (times L4) for the system 2/d. At the center of this debate is the way the average
sizes we studied. This quantity is calculated during the over the disorder is performed. If the disorder is cho-
simulation without the use of any numerical derivatives. sen from a random distribution without any constraints,
As explainedinsectionIIIAwe expectthis derivativeto calledthecanonicalensembleofdisorder,onemightaskif
behave as L1/ν at K . The results for K = 0.188 equivalentresultsareobtainedifconstraintsareimposed
c
∼
clearly deviate from a power-law consistent with this ontherandomdistribution,theso-calledmicro-canonical
value of K being below K . From the results shown in ensembleofdisorder. Forinstance,forthemodelconsid-
c
Fig. 3 we know that K =0.191 is likely above K , while eredhereonecouldconstraintherandomchemicalpoten-
c
K = 0.190 is likely very close to K . Hence we fit the tial to have exactly the same average for each generated
c
results for K = 0.190 and K = 0.191 to a power-law, sample. The proof22 of the “quantum” Harris criterion
L4dρ/dK =cL1/ν, finding anexponent of 1.11and 1.69, reliesonthephysicallymorerelevantcanonicalensemble
respectively (shown as the solid lines in Fig. 7). This ofdisorderbeingused. Subsequentwork40,41,42,43,44 have
allows us to bracket the estimate of ν to the interval shown that even though the two ensembles in principle
0.59 0.91. The rather broad range of this interval is yield the same results in the thermodynamic limit, their
−
due to the fact that the relatively large value of z makes finite-size properties can in some cases be different. All
thiswayofdeterminingν extremelysensitivetoaprecise our results have been obtained using the canonical en-
determinationofK . Wethereforecombinethisestimate semble of disorder, in that at each site a potential was
c
with a standard scaling plot of ρL4 versus δL1/ν shown drawn from the uniform distribution [0,1]. Hence, for a
in the inset of Fig. 7. The best scaling plot is obtained givensampletheaveragechemicalpotentialisnotexactly
with ν = 0.70 and K = 0.1902 (well within the er- 1/2. In light of the above discussion we have therefore
c
8
performed additional simulations directly at K but this largerlatticesizesforthespatialcorrelationfunctionsby
c
time imposing the constraint that the average chemical decreasing the aspect ratio dramatically. We have at-
potential must be exactly 1/2 for every disorder realiza- tempted this by performing calculations for a lattice of
tion. Our results for this micro-canonical ensemble of size 203 40 at K (also shown in Fig. 8), correspond-
c
×
disorderarealsoshowninFig7andareindistinguishable ing to an aspect ratio of α = 1/200. In this case the
from our results obtained with the canonical ensemble. spatial correlation functions show clear evidence for an
One should note that the uniform disorder distribution exponential decay probably due to a dimensional cross-
weemployislikelylesssensitivetothedifferencebetween overinducedbytheextremelysmallaspectratio. Infact,
the canonical and micro-canonical ensemble of disorder it would seem to be more reasonable to increase α to a
thanabimodaldistribution. Henceweconcludethatour moreoptimalaspectratiodeterminedbytherequirement
procedure for calculating ν is valid. that (ξ/L)z ξ /L . Due to the large value of z, this
τ τ
∼
We finally discuss our results for the correlation func- has proven impossible for the present study.
tions which are shown in Fig. 8 for a range of different V. CONCLUSION
system sizes calculated at K . To determine the power
c
law decay of the correlationfunctions we used lattices of
In the present paper we have shown strong numerical
size L3 αLz with α=1/8 and z =3. All the temporal
evidence that the dynamical critical exponent z at the
×
correlation functions can be fitted with the same expo-
BG-SF critical point is equal to the spatial dimension
nent y =1.06(7). We note that this value for y would
τ τ for the three-dimensional site-disordered Bose-Hubbard
appear to satisfy the equality y 1 as an equality. Fol-
τ model. The relation z = d therefore continues to hold
≤
lowingsectionIIIAweassumethaty =(d+z 2+η)/z
τ in d = 3. Results for the scaling of the stiffness versus
−
and using the above determined value for z =3 we then
K, as well as at K =0.190(1) versus L , are consistent
c τ
infer
with z = d = 3. The compressibility is almost constant
andindependentofsystemsizethroughthecriticalpoint
η 1. (22)
∼− consistent with this value for z. In addition we have ob-
tained values for the critical exponents ν =0.70(12)and
This would then satisfy Eq. (15) as an equality. Our
η 1, with the cautionary note that a reliable deter-
results for the temporal correlation functions have been ∼ −
minationofη isimpededbytheverysmallspatialextent
determinedouttolargelatticesizesL =512andappear
τ of the lattices available. In light of the recent work by
quite stable. We are therefore relatively confident that
Bernardet et al.39, it would have been very interesting
these results are trustworthy.
to analyze each disorder realization seperately, in order
Our results for the spatial correlation functions
C(r,τ = 0) appear less clear. The value we obtain to determine a characteristic Ki for each sample. The
scaling analysis should then be redone using K . Unfor-
for y = 4.4(4) is clearly not consistent with our other i
r tunately, we were not able to perform such an analysis
estimates for z, since the ratio y /y z then yields
r τ ≡ for this work.
z = 4.15. If one assumes that this is the correct value
for z it again follows that η 1. We strongly suspect
∼ −
that this value of y is due to the relativelysmall spatial
r
Acknowledgments
extent (L 16) of the lattices used. For such small lat-
≤
ticesizesthecorrelationfunctionshavelikelynotreached
their asymptotic behavior. We have varied the strength AllourcalculationshavebeenperformedattheSHAR-
of the disorder, and found the same critical behaviour Cnet computational facility and we thank NSERC and
at ∆ = 0.6 and ∆ = 1.0. One might attempt to reach CFI for financial support.
∗ Electronic address: [email protected] and A. M. Goldman, Phys.Rev.B 60, 4320 (1999).
† Electronic address: [email protected] 7 M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D.S.
1 M. Greiner, O. Mandel, T. Esslinger, T. W. H¨ansch, and Fisher, Phys.Rev.B. 40, 546 (1989).
I. Bloch, Nature415, 39 (2002). 8 J.E.Lye,L.Fallani,M.Modugno,D.S.Wiersma,C.Fort,
2 P.A.Crowell,F.W.VanKeuls,andJ.D.Reppy,Physical and M. Inguscio, Phys. Rev.Lett. 95, 070401 (2005).
Review B 55, 12620 (1997). 9 T. Giamarchi and H. J. Schulz, Phys. Rev. B 37, 325
3 H. S. J. van der Zant, F. C. Fritschy, W. J. Elion, L. J. (1988).
Geerligs, and J. E. Mooij, Phys. Rev. Lett. 69, 2971 10 W. Krauth,N.Trivedi, and D.Ceperley, Phys.Rev.Lett.
(1992). 67, 2307 (1991).
4 Y. Liu, D. B. Haviland, B. Nease, and A. M. Goldman, 11 E.S.Sørensen,M.Wallin,S.M.Girvin,andA.P.Young,
Phys. Rev.B 47, 5931 (1993). Phys. Rev.Lett. 69, 828 (1992).
5 N. Markovi´c, C. Christiansen, and A. M. Goldman, Phys. 12 G. G. Batrouni, B. Larson, R.T. Scalettar, J. Tobochnik,
Rev.Lett. 81, 5217 (1998). and J. Wang, Phys.Rev.B 48, 9628 (1993).
6 N. Markovi´c, C. Christiansen, A. M. Mack, W. H. Huber, 13 M.Makivic,N.Trivedi,andS.Ullah,Phys.Rev.Lett. 71,
9
2307 (1993). 30 N. Prokof’ev and B. Svistunov, Physical Review Letters
14 D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and 92, 015703 (2004).
P. Zoller, Phys. Rev.Lett. 81, 3108 (1998). 31 M.Wallin,E.Sørensen,S.Girvin,andA.P.Young,Phys.
15 R.RothandK.Burnett,Phys.Rev.A 68,023604(2003). Rev.B 49, 12115 (1994).
16 B. Damski, J. Zakrzewski, L. Santos, P. Zoller, and 32 M. Guo, R. N. Bhatt, and D. A. Huse, Phys. Rev. Lett.
M. Lewenstein, Phys. Rev.Lett. 91, 080403 (2003). 72, 4137 (1994).
17 C. Fort, L. Fallani, V. Guarrera, J. E. Lye, M. Modugno, 33 M. Guo, R. N. Bhatt, and D. A. Huse, Phys. Rev. B 54,
D. S. Wiersma, and M. Inguscio, Phys. Rev. Lett. 95, 3336 (1996).
170410 (2005). 34 R.N.BhattandA.P.Young,PhysicalReviewB37,5606
18 F. Alet and E. Sørensen, Phys. Rev. E 67, 015701(R) (1988).
(2003). 35 P. Hitchcock and E. S. Sørensen, preprint, cond-
19 S. N.Dorogovtsev, Phys.Lett. A 76, 169 (1980). mat/0601180 (2006).
20 P.B.WeichmanandK.Kim,PhysicalReviewB40,R813 36 F.Alet and E.Sørensen, Phys.Rev.B 70, 024513 (2004).
(1989). 37 A. Aharony and A. B. Harris, Phys. Rev. Lett. 77, 3700
21 R. Mukhopadhyay and P. B. Weichman, Physical Review (1996).
Letters 76, 2977 (1996). 38 F. P´azm´andi, R. T. Scalettar, and G. T. Zim´anyi, Phys.
22 J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, Rev.Lett. 79, 5130 (1997).
Phys. Rev.Lett. 57, 2999 (1986). 39 K. Bernardet, F. P´azm´andi, and G. G. Batrouni, Phys.
23 I. F. Herbut,Phys.Rev. B 58, 971 (1998). Rev.Lett. 84, 4477 (2000).
24 I. F. Herbut,Phys.Rev. B 61, 14723 (2000). 40 A. Aharony, A. B. Harris, and S. Wiseman, Phys. Rev.
25 F. Alet and E.Sørensen, Phys.Rev.E 68, 026702 (2003). Lett. 81, 252 (1998).
26 F. P´azm´andi and G. T. Zim´anyi, Phys. Rev. B 57, 5044 41 S. Wiseman and E. Domany, Phys. Rev. Lett. 81, 22
(1998). (1998).
27 L. Zhang and M. Ma, Phys. Rev.B 45, 4855 (1992). 42 F. Igl´oi and H. Rieger, Phys. Rev.B 57, 11404 (1998).
28 K. G. Singh and D. S. Rokhsar, Phys. Rev. B 46, 3002 43 A.DharandA.P.Young,Phys.Rev.B 68,134441(2003).
(1992). 44 C. Monthus, Phys. Rev.B 69, 054431 (2004).
29 J. KiskerandH.Rieger, Phys.Rev.B 55, R11981 (1997).