Delocalization of Weakly Interacting Bosons in a 1D Quasiperiodic Potential

V. P. Michal,1 B. L. Altshuler,2 and G. V. Shlyapnikov1,3,41Laboratoire de Physique Théorique et Modèles Statistiques, Université Paris Sud, CNRS, 91405 Orsay, France

2Physics Department, Columbia University, 538 West 120th Street, New York, New York 10027, USA3Van der Waals-Zeeman Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

4Russian Quantum Center, Novaya street 100, Skolkovo, Moscow region 143025, Russia(Received 7 March 2014; revised manuscript received 23 June 2014; published 23 July 2014)

We consider weakly interacting bosons in a 1D quasiperiodic potential (Aubry-Azbel-Harper model) inthe regime where all single-particle states are localized. We show that the interparticle interaction may leadto the many-body delocalization and we obtain the finite-temperature phase diagram. Counterintuitively,in a wide range of parameters the delocalization requires stronger coupling as the temperature increases.This means that the system of bosons can undergo a transition from a fluid to insulator (glass) stateunder heating.

DOI: 10.1103/PhysRevLett.113.045304 PACS numbers: 67.85.-d, 03.75.Lm, 05.60.-k

Quantum mechanics of a single particle in a quasiperi-odic potential is a standard although quite nontrivialtheoretical problem [1–6] intensively studied for severaldecades. It turned out that for one-dimensional (1D)quasiperiodic potentials (superposition of two incommen-surate periodic potentials) the eigenfunctions can be eitherextended or localized, depending on the parameters of thepotential. The phenomenon of quantum localization wasobserved for 1D quantum gases in both random [7] andquasiperiodic potentials [8]. The experiments conductedin the regime of negligible interaction between the atomsof expanding Bose gas demonstrated good qualitativeagreement with the single-particle theories of localization[3,9]. The same applies to spreading of wave packets oflight in a quasiperiodic photonic lattice [10]. Recently[11,12], a Feschbach resonance was used to study bosonswith sizable and fine-tunable interaction in the quasiperi-odic potential.Theoretical description of many-body effects in disor-

dered fermionic [13,14] and bosonic [15] systems is basedon the idea of the localization of many-body wavefunctions in the Hilbert space—many-body localization.The quasiperiodic potential represents an intermediate casebetween periodic and disordered systems. Zero temperaturephase diagram for 1D bosons in the quasiperiodic potentialhas been previously discussed and calculated numerically[16–18], and the case of an infinite temperature for(nearest-neighbor) interacting spinless fermions has beenstudied in Ref. [19]. The problem of localized and extendedstates of two interacting particles in the 1D quasiperiodicpotential has been discussed in Refs. [20,21].In this Letter we study finite-temperature transport proper-

ties of interacting bosons in the 1D quasiperiodic potentialand predict the physical behavior which differs drasticallyfrom the many-body localization of bosons caused byrandom potentials. We show that, counterintuitively, in a

broad temperature range an increase in temperature induces atransition from fluid to glass.The standard model of a 1D quasiperiodic potential is the

Aubry-Azbel-Harper (AAH) model [1–3]—a tight-bindingHamiltonian with hopping amplitude J and periodicallymodulated on-site energies, at a period incommensuratewith respect to the primary lattice. The eigenstate ψα

j atenergy εα is determined by the equation

Jðψαjþ1 þ ψα

j−1Þ þ V cosð2πκjÞψαj ¼ εαψ

αj : ð1Þ

Here V is the modulation amplitude, and κ is an irrationalnumber. In 1D random potentials all single-particle statesare localized [22,23]. On the contrary, in the AAH modelall states are extended unless the amplitude of the modu-lation exceeds a critical value 2J. Then all states arelocalized, and the localization length in units of the latticeconstant is given by [3]

ζ ¼ 1

lnðV=2JÞ ⇒ ζ ≃ VV − 2J

≫ 1 for V − 2J ≪ V: ð2Þ

Below we assume ζ ≫ 1. Our analytical considerationbased on the semiclassical approach [2,5,24,25] to thesingle-particle problem, is valid when the period of themodulation is much larger than the period of the lattice,κ ≪ 1. We supplement this analysis by numerical calcu-lations for κ ∼ 1, in particular for κ equal to the goldenratio ð ffiffiffi

5p

− 1Þ=2.The semiclassical one-particle spectrum is organized

according to the continued fraction decomposition [2,5]of the irrational parameter κ ¼ 1=½n1 þ 1=ðn2 þ…Þ�. Forn1; n2;… ≫ 1 the spectrum has a hierarchical structure: itconsists of n1 narrow first-order bands (FOBs)—clusters ofL=n1 energy levels (L is the size of the system). Each FOBcontains n2 second-order bands, so that there are ∼n1n2second-order bands in total, etc. All eigenstates are located

PRL 113, 045304 (2014) P HY S I CA L R EV I EW LE T T ER Sweek ending25 JULY 2014

0031-9007=14=113(4)=045304(5) 045304-1 © 2014 American Physical Society

in the energy interval −V − 2J < ε < V þ 2J. The spacingω between FOBs is the frequency of the classical periodicmotion [26] in a single potential well of the size 1=κ ≈ n1[27], which includes n1 levels,

ω ¼ 2πHdx=v

≃ 2π2κVlnð64V2=jε2 − V2=ζ2jÞ ; ð3Þ

where v ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi4J2 − ½ε − V cosð2πκxÞ�2

pis the classical

velocity of a particle. The widths of FOBs are determinedby the tunneling between neighboring wells [5],

Γs ¼4ω

πexp

�−Z

dxjpj�: ð4Þ

The integral is taken over the classically forbidden region,and jpj ¼ arccosh½ðV cosð2πκxÞ − εsÞ=2J�. The index slabels FOBs centered at the energies εs. The action canbe approximated as

Rdxjpj ≈ jεsj=4κJ, which yields an

exponential dependence of the bandwidth on energy

Γs ≈32κJπ

expð−jεsj=4κJÞ: ð5Þ

Can the localization of the bosons be destroyed by theinteraction? It is known that the interaction can delocalizefermions [14] and 1D bosons [15] in the case of a randompotential. Experiments with interacting bosons in 1Dquasiperiodic potentials [11,12] indicated an interaction-induced localization-delocalization transition. It is alsoworth noting that experiments in quasiperiodic photoniclattices [10] have found that nonlinearity (interactions)increases the width of localized wave packets of light.Here we consider the AAH model with a weak on-siteinteraction,

Hint ¼U2

Xj

a†ja†jajaj; U ≪ J; ð6Þ

with aj being the bosonic field operators. In order toestimate the critical coupling constant Uc correspondingto the many-body localization-delocalization transition(MBLDT), we use the method developed in [14,15], whichis similar to the original estimation for the single-particleAnderson localization [9]. One has to consider the localizedone-particle states jαi and analyze how different two-particle states jα; βi hybridize due to the interaction. Thecriterion of MBLDT is

Pα ∼ 1; ð7Þ

where Pα is the probability that for a given one-particlestate jαi there exist three other states jβi; jγi; jδi, such thatthe two-particle states jα; βi and jγ; δi are in resonance; i.e.,the matrix element hγ; δjHintjα; βi≡Mγδ

αβ exceeds the

energy mismatch Δγδαβ ≡ jεα þ εβ − εγ − εδj where εα, εβ,

εγ and εδ are one-particle energies.For large occupation numbers Nβ, Nγ, and Nδ of the

states jβi, jγi and jδi the fluctuations are small. Selecting agiven single-particle state α and taking into account bothdirect and inverse processes we find [28]

Mγδαβ ¼

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffijNβð1þNγÞð1þNδÞ−NγNδð1þNβÞj

qUγδ

αβ; ð8Þ

with

Uγδαβ ¼ U

Xj

ψδ�j ψγ�

j ψβjψ

αj : ð9Þ

As discussed in Refs. [14,15], the matrix elements of theinteraction are small unless the energies εα; εβ; εγ; εδ arealmost equal pairwise, e.g., εα ≈ εγ and εβ ≈ εδ, while εαand εβ can differ substantially. Accordingly,M

γδαβ ≈ NβU

γδαβ.

The approximation (8) for the matrix elements remainsvalid for small occupation numbers.If α and γ (as well as β and δ) are nearest neighbors in

energy the energy mismatch is

Δγδαβ ¼ δα þ δβ; ð10Þ

where δα ¼ jεα − εγj is a typical spacing between the stateson the length scale ζ at energy close to εα. We estimate thematrix element [29]:

Mγδαβ ≈UNβ=ζ: ð11Þ

According to Eqs. (10) and (11) the probability Pγδαβ of

having Mγδαβ ≳ Δγδ

αβ is

Pγδαβ ≈UNβ=ζðδα þ δβÞ: ð12Þ

The probability Pα which enters the criterion (7) ofMBLDT is the sum of Pγδ

αβ over all single-particle statesjβi; jγi; jδi. Since for given α and β the number of relevantpairs of states jγi; jδi is of order unity, only the summationover β is important,

Pα ¼Xβ;γ;δ

Pγδαβ ≈

X0

β

UNβ=ζðδα þ δβÞ; ð13Þ

whereP0

means that the summation is over the eigenstateson the length scale of ζ. Substitution of Eq. (13) into Eq. (7)leads to the criterion of MBLDT:

X0

β

UcNβ=ζðδα þ δβÞ ¼ 1: ð14Þ

The critical coupling strengthUc in Eq. (14) depends on thechoice of the state α through the quantity δα. One has to

PRL 113, 045304 (2014) P HY S I CA L R EV I EW LE T T ER Sweek ending25 JULY 2014

045304-2

choose α which minimizes Uc. Besides, Eq. (14) expressesthe critical couplingUc in terms of the occupation numbersNβ which are determined by the chemical potential. Inorder to formulate the MBLDT criterion in terms of theexperimentally controllable filling factor ν one has tocomplement Eq. (14) with the number equation, whichrelates ν to the occupation numbers:

ν ¼Xα

Nα=L: ð15Þ

We are now ready to apply the MBLDT criterion tointeracting bosons in the AAH model. On the insulatorside, the occupation numbers are given by

Nα ¼ fexp½ðεα − μþUNα=ζÞ=T� − 1g−1; ð16Þ

where μ is the chemical potential shifted by the interactionenergy of a particle in the state jαiwith particles in all otherstates. Equations (14), (15), and (16) yield the criticalcoupling Uc as a function of temperature.On the length scale ζ there are ∼ζ states with a

significant amplitude of the wave function. Thus, thereare at most ζ states contributing to the sum in Eq. (14). Thecalculation simplifies in the limit n1 ≪ ζ ≪ n1n2 whereeach FOB contains ∼ζ=n1 ≈ κζ overlapping states. Therelated example is shown in Fig. 1: κ is close to 1=8, thelocalization length is ζ ≃ 40, and n2 ¼ 7.At temperatures much smaller than the spacing between

the FOBs, i.e., T ≪ ω, only single-particle states from thelowest energy FOB participate in MBLDT. We assume thatthe spacing between the states in this band is approximatelyconstant and thus equal to δβ ≈ Γ0=κζ. Using the fact that

the sumP0βNβ=ζ over the states on the length scale ζ in

Eq. (14) is equal to the sumP

αNα=L in Eq. (15) over allstates we obtain [30]:

νUc ≈ 2Γ0=κζ; T ≪ ω: ð17ÞAt temperatures T ≫ ω many FOBs are occupied and

particles in these bands participate in MBLDT. The bandsare so narrow that all levels in the sth FOB have the sameoccupation, so that the corresponding level spacing isδβ ≈ Γs=κζ. Then, summing over β within each FOBone can rewrite Eq. (14) as

Xn1−1s¼0

ζUcNsκ2=ðΓ0 þ ΓsÞ ≈ 1 ð18Þ

(we selected the state jαi to be in the lowest energy band,s ¼ 0). According to Eq. (5) the width Γs exponentiallyincreases with εs (i.e., s) in the interval 0 > εs > ε0≃−V − 2J, since jεsj decreases. Hence, at T ≪ 8J the sumover s in Eq. (18) is dominated by s ¼ 0, which leads to

Uc ≈ 2Γ0=N0κ2ζ: ð19Þ

For εs − μ < T we may use the occupation numberexpression [29]

Ns ≈ T=ðεs − μÞ; ð20Þand put Ns ≈ 0 for larger εs. In particular N0 ≈ T=ðε0 − μÞ.With the use of Eq. (20), one can rewrite Eq. (15) as

ν ≈XðTþμ−ε0Þ=ω

s¼0

κT=ðsωþ ε0 − μÞ: ð21Þ

In the temperature range ω ≪ T ≪ 8J we then find thechemical potential dependence on temperature [29]

μ ≈ ε0 −κTν

�1þ T

8νJln

�Tω

��; ð22Þ

and using equation (19) we obtain the critical coupling:

νUc ≈2Γ0

κζ

�1þ T

8νJln

�Tω

��: ð23Þ

Since T ≪ 8J, the second term in square brackets is asmall correction. Nevertheless, it is important. According toEq. (23) the temperature dependence of the critical cou-pling is anomalous: Uc increases with T; i.e., an increase intemperature favors the insulator state.This behavior originates from the cluster structure of the

spectrum, with exponentially increasing cluster width whengoing from the lowest (highest) cluster energy to the middleof the spectrum. Therefore, at T ≪ 8J (and ζ ≫ n1 ≈ κ−1)

FIG. 1 (color online). The critical coupling strength Uc versustemperature obtained by directly using Eqs. (14)–(16) and theone-particle spectrum computed by exact diagonalization, for κclose to 1=8, V ¼ 2.05J, and the filling factor ν ¼ 1; 2; 4. AtT ¼ 0 we recover universality in νUc0. The dashed line is theT → ∞ asymptotics. The inset shows the spectrum on the lengthscale of ζ (the number of states is ζ ≈ 40).

PRL 113, 045304 (2014) P HY S I CA L R EV I EW LE T T ER Sweek ending25 JULY 2014

045304-3

the localization-delocalization transition is provided onlyby the particle states in the lowest energy cluster. Thefraction of these particles decreases with increasing temper-ature, thus ensuring an increase in the critical couplingstrength Uc.For temperatures T → ∞ (T ≫ 8J; 8νJ) all eigenstates

are equally populated and Ns ¼ ν. Then the main con-tribution to the sum in Eq. (18) comes from s ¼ 0 ands ¼ n1 − 1 as Γn1−1 ≈ Γ0 ≈ ð32κJ=πÞ expð−1=κÞ. Havingin mind Eq. (17) we thus obtain

νUc∞ ≃ Γ0=κ2ζ ≃ νUc0=2κ ≫ νUc0; ð24Þ

where Uc0 is the zero temperature critical coupling atζ ≫ n1 ≈ κ−1. Therefore, under this condition we alwaysexpect the anomalous ”freezing with heating” behavior athigh enough temperatures.Our analytical results are confirmed by numerics using

the single-particle spectrum obtained by exact diagonaliza-tion. The results for κ ≈ 1=8 and ζ ≃ 40 are shown inFig. 1. The critical coupling turns out to be very small sinceit is proportional to the width Γ0, which is several orders ofmagnitude smaller than J. This justifies the validity of ourperturbative approach with respect to the interparticleinteraction.In the opposite limit 1 ≪ ζ ≲ κ−1, the situation changes.

Indeed in this case single-particle states participating inMBLDT belong to different FOBs (not more that one statefrom a given cluster). The characteristic spacing betweenthese states is ∼8J=ζ and the cluster structure of thespectrum is not important. The resulting critical couplingνUc is ∼J. For κ ≪ 1 one can use the quasiclassicalapproach with the density of states κ=ωðεÞ and ωðεÞ givenby Eq. (3). The results are consistent with our calculationsbased on exact diagonalization for the one-particle spec-trum and Eqs. (14), (15), and (16). They suggest a slowdecrease of Uc with increasing temperature and are

displayed in the inset of Fig. 2(a) for ζ and κ−1 both closeto 8. However, since our approach is based on theperturbative treatment of the interactions, its predictionsat νUc ≳ J at least require a large filling factor ν. A detailedanalysis of this question will be given elsewhere.For κ ∼ 1 the quasiclassical approach is no longer valid

and one has to rely only on the numerics based on exactdiagonalization for the one-particle problem and Eqs. (14),(15), and (16). The results for κ ¼ 1=ð1þ 1=ð1þ…ÞÞ ¼ð ffiffiffi

5p

− 1Þ=2 (golden ratio) and for κ close to 0.24 at ζ ≈ 7

are shown in Fig. 2. For the latter case our results at T ¼ 0

and ν ¼ 1 are consistent with the DMRG calculations ofRefs. [16] and [17] (using κ ¼ 0.77 which is equivalent toκ ¼ 0.23) extrapolated to V ¼ 2.3J.Reference [19] presented results of the numerical sim-

ulation for spinless fermions with nearest neighbor inter-action subject to a quasiperiodic potential at T ¼ ∞. Thisproblem (different from bosons with the onsite interaction),can also be attacked with our approach at any temperature.The results will be published elsewhere, but alreadynow we can say that at T → ∞, they agree fairly wellwith Ref. [19].Our results at finite temperatures indicate an anomalous

UcðTÞ dependence. The experiment [11] has been per-formed for κ ≈ 1.24, which according to Eq. (1) is equiv-alent to κ ¼ 0.24. The extrapolation of experimental resultsto V ≈ 2.3J gives νUc=J ∼ 0.3, which is consistent withour calculations.In conclusion, we have developed the many-body

localization theory of weakly interacting bosons in a 1Dquasiperiodic potential and obtained the phase diagram interms of temperature and interaction. The most unexpectedprediction based on our calculations is the transition fromfluid to insulator (glass) with heating.

We are grateful to I. L. Aleiner and G. Modugno forfruitful discussions and acknowledge support from Triangle

FIG. 2 (color online). The same as in Fig. 1 for V ¼ 2.3J (ζ ≈ 7). In (a) κ ≈ 0.24, and in (b) κ is equal to the golden ratio. The inset in(a) shows νUcðTÞ for κ ≈ 1=8 and V ¼ 2.25J (ζ ≈ 8).

PRL 113, 045304 (2014) P HY S I CA L R EV I EW LE T T ER Sweek ending25 JULY 2014

045304-4

de la Physique through the project DISQUANT, IFRAF,and the Dutch Foundation FOM. The research leading tothese results has received funding from the EuropeanResearch Council under European Community’s SeventhFramework Programme (FR7/2007-2013 Grant AgreementNo. 341197).

[1] P. G. Harper, Proc. Phys. Soc. London Sect. A 68, 874(1955).

[2] M. Ya. Azbel, Sov. Phys. JETP 17, 665 (1963): Sov. Phys.JETP 19, 634 (1964); Phys. Rev. Lett. 43, 1954 (1979).

[3] S. Aubry and G. André, Ann. Isr. Phys. Soc. 3, 133(1980).

[4] J. B. Sokoloff, Phys. Rev. B 22, 5823 (1980); 23, 6422(1981); J. B. SokoloffPhys. Rep. 126, 189 (1985).

[5] I. M. Suslov, Sov. Phys. JETP 56, 612 (1982). In this paperthe notation Γ s is used for the halfwidth.

[6] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. denNijs, Phys. Rev. Lett. 49, 405 (1982).

[7] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P.Lugan, D. Clement, L. Sanchez-Palencia, P. Bouyer, and A.Aspect, Nature (London) 453, 891 (2008).

[8] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M.Zaccanti, G. Modugno, M. Modugno, and M. Inguscio,Nature (London) 453, 895 (2008).

[9] P. W. Anderson, Phys. Rev. 109, 1492 (1958).[10] Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, N.

Davidson, and Y. Silberberg, Phys. Rev. Lett. 103, 013901(2009).

[11] L. Tanzi, E. Lucioni, S. Chaudhuri, L. Gori, A. Kumar, C.D’Errico, M. Inguscio, and G. Modugno, Phys. Rev. Lett.111, 115301 (2013).

[12] C. D’Errico, E. Lucioni, L. Tanzi, L. Gori, G. Roux, I. P.McCulloch, T. Giamarchi, M. Inguscio, and G. Modugno(to be published).

[13] B. L. Altshuler, Yu. Gefen, A. Kamenev, and L. S. Levitov,Phys. Rev. Lett. 78, 2803 (1997).

[14] D. Basko, I. L. Aleiner, and B. L. Altshuler, Ann. Phys.(Amsterdam) 321, 1126 (2006); Problems of CondensedMatter Physics, edited by (Oxford University, New York,2007), Chap. 5.

[15] I. L. Aleiner, B. L. Altshuler, and G. V. Shlyapnikov, Nat.Phys. 6, 900 (2010).

[16] G. Roux, T. Barthel, I. P. McCulloch, C. Kollath, U.Schollwöck, and T. Giamarchi, Phys. Rev. A 78, 023628(2008).

[17] X. Deng, R. Citro, A. Minguzzi and E. Orignac, Phys. Rev.A 78, 013625 (2008).

[18] T. Roscilde, Phys. Rev. A 77, 063605 (2008).[19] S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Phys.

Rev. B 87, 134202 (2013).[20] S. Flach, M. Ivanchenko, and R. Khomeriki, Europhys. Lett.

98, 66002 (2012).[21] G. Dufour and G. Orso, Phys. Rev. Lett. 109, 155306

(2012).[22] M. E. Gertsenshtein and V. B. Vasiliev, Theory Probab.

Appl. 4, 391 (1959).[23] N. F. Mott and W. D. Twose, Adv. Phys. 10, 107 (1961).[24] M. Wilkinson, Proc. R. Soc. A 391, 305 (1984).[25] M. Albert and P. Leboeuf, Phys. Rev. A 81, 013614 (2010).[26] L. D. Landau and E. M. Lifshitz, Quantum Mechanics

(Non-relativistic Theory), Course of Theoretical Physics,3rd ed. (Elsevier, New York, 1977), Vol. 3.

[27] Hereinafter considering the limitκ ≪ 1, we put everywheren1 ¼ κ−1.

[28] The matrix element hγ; δjHintjα; βi by itself is proportionalto

ffiffiffiffiffiffiNα

p. However, the related probability (proportional to

Nα) determines the time derivative of Nα, so that Nα dropsfrom the expression for the probability of hybridization of agiven state of the Nα manifold. This is equivalent toconsidering a single particle state jαi with Nα ¼ 1.

[29] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.113.045304 for detailson the approximate exponential energy dependence of theFOB widths, evaluation of the overlap sums, and calculationof the critical interaction temperature dependence.

[30] In fact, the level spacing δβ is minimal in the center of theFOB and somewhat increases closer to the band edges. Thisslightly increases the critical coupling compared to the resultof Eq. (17). However, the quantity νUcκζ=2Γ0 remainsequal to unity within a factor of 2 up to κ ≈ 0.25, where thecluster structure of the spectrum is still pronounced but thequasiclassical approach is strictly speaking not valid. Adecrease in the level spacing with increasing energy in thelow-energy part of the FOB can also influence the temper-ature dependence of Uc at T < ω. Indeed, as shown inFig. 1, the curveUcðTÞ has a dip at very low temperatures. IfT < ω and the occupation of the lowest energy FOB islarge, the critical coupling Uc first decreases with increasingtemperature (this is analogous to the behavior of UcðTÞdisplayed in the inset of Fig. 2, where not more than onestate of a given FOB participates in MBLDT). Once thetemperature exceeds the spacing between the FOBs theexponential dependence of the FOB widths dominates andthe temperature dependence of Uc becomes anomalous.

PRL 113, 045304 (2014) P HY S I CA L R EV I EW LE T T ER Sweek ending25 JULY 2014

045304-5