Threshold for Chaos and Thermalization in One-Dimensional Mean-Field Bose-Hubbard model Amy C. Cassidy,1,2 Douglas Mason,3 Vanja Dunjko,2 and Maxim Olshanii2 1Department of Physics and Astronomy, University of Southern California, Los Angeles, CA 90089, USA 2Department of Physics, University of Massachusetts Boston, Boston MA 02125, USA 3Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA (Dated: January 23, 2009) We study the threshold for chaos and its relation to thermalization in the 1D mean-field Bose- 9 Hubbard model, which in particular describes atoms in optical lattices. We identify the threshold 0 for chaos, which is finite in the thermodynamic limit, and show that it is indeed a precursor of 0 thermalization. Farabovethethreshold,thestateofthesystem afterrelaxation isgovernedbythe 2 usual laws of statistical mechanics. n PACSnumbers: 67.85.-d,03.75.Lm,05.45.Jn a J 3 Introduction.– We study the threshold for chaos and chaos threshold is governed only by the parameters and 2 the ability to thermalize of the 1D mean-field Bose- observables that are finite in the thermodynamic limit, Hubbard model (BHM) [1]. The study of thermaliza- and as a result it remains finite in that limit. ] r tion in non-linear systems dates back to the early work Ourprogramissimilartoacomprehensivecomparison e h of Fermi, Pasta, and Ulam (FPU) [2] on a non-linear between FPU and φ4 [12], where, however,the existence t string,modeledbyanharmonicallycoupledoscillators. It of the thermalization threshold in FPU is denied. o was expected that for a large number of degrees of free- In this paper we observe thermal behavior in time- . t a dom,evensmallnonlinearitieswouldcausethesystemto averagedmean-fieldquantities. Note thatinrecentwork m thermalize, resulting in energy equipartition. However, on thermalization in quantum systems, thermal proper- - equipartitionwasnotobserved. The absence ofthermal- ties emerge from individual quantum stationary states d ization was eventually explained in two complementary [13,14]. Studiesonthesemi-classicalregimesuggestthat n ways: one in terms of closeness to an integrable system, thetwoarerelated,althoughopenquestionsremain[13]. o the Korteweg-de Vries model [3], and another in terms Empirically, our system describes the motion of c [ of a chaos threshold given by the theory of overlapping bosonic atoms in a one-dimensional tight-binding opti- resonances put forth by Chirikov and Israilev [4, 5]. cal lattice [1, 15]. 3 v Since then further studies on thermalization and ap- Systemofinterest.–Westudythemean-fielddynamics 8 proach to equilibrium have been carried out in several of an interacting one-dimensional Bose gas on a lattice 8 classical field theories, including recent studies on the (1D Bose-Hubbard model (BHM)) with periodic bound- 3 classical φ4 model [6], Nonlinear Klein-Gordon equation ary conditions. The Hamiltonian in the momentum rep- 3 (NLKG) [7], Non-Linear Schr¨odinger equation (NLSE) resentation is . 5 [8,9],DiscreteNon-LinearSchr¨odingerequation(DNLS) µ 80 [9, 10] equivalent to BHM, and Integrable Discrete Non- H =X(cid:16)¯hωn|ψn|2− 20|ψn|4(cid:17)+ (1) Linear Schr¨odinger equation (IDNLS)[9]. n 0 µ : No conventional thermalization is expected in the 0 ψ∗ψ∗ψ ψ , v 2 X i j n i+j−n i NLSEandIDNLS,whicharebothintegrable. InNLKG, i,j,n X like in FPU, the ability of the system to reach thermal n6=i,j r equilibrium in the course of time evolution emerges only a where the indices span the range n, i, j = when the degree of nonlinearity exceeds a certain criti- 0, 1, 2, ..., Ns−1 (N is supposed to be odd). cal value (see [5, 11] for the thermalization threshold in ± ± ± 2 s Throughout the text the wavefunction ψ is normalized FPU). On the contrary,the φ4 model eventually reaches n to unity: ψ 2 = 1. The bare frequency of each equilibrium regardless of how small the nonlinearity is. Pn| n| momentum mode is given by ω = 2Jcos 2πn , and In our paper we show that the BHM (along with the n − (cid:16) Ns (cid:17) equivalent DNLS) belongs to the former class. the coupling constant is µ0 = UNa/Ns. Here J and Furthermore,wehavecomparedtwoquantitativemea- U are the nearest-neighbor site-hopping and on-site sures of thermalizability: maximal Lyapunov exponent repulsion constants of the standard Bose-Hubbard (whose positivity is a signature for chaos) and spectral model, respectively, and Na is the number of atoms. erinutmro)p.yB(wothhichmeparsouvridesesshaowdisatanshcaerptotthhreersmhoalldeaqsuiolinbe- Teqhueatciaonnosnoicfaml potaiiorns aarreeQginve=n bψyn,∂∂Ptψnn==i¯h−ψh¯n∗i,∂∂ψaHn∗n.d Wthee varies the nonlinearity strength, and the two thresholds define the dimensionless non-linearity parameter, κ, are undeniably close. Furthermore, we assert that the to be the ratio between the typical interaction energy 2 per site, U(Na/Ns)2, and the hopping energy per site, 1 κ=0.36 0.5 JCeoxfNhpaaoaon/rseNegnsict:orsni.–teinTrihopnehκassta≡aennsdµdpJa0acrch≡deasioUsisgt(nhNatahJtatru/etrNshehseoo)flsdetph.aefrrocahmtaioonLtiycbaenptawutneu(eo2rnve) ηSpectral entropy, 0000....2468 η 0000 ....02468 0 0κκκ.===00004λ...579 [42J ]0.08 0.12 ηλ 0000....1234 λapunov exponent, [J] y initially closetrajectoriesgrowsexponentiallywithtime, L for typical trajectories, as captured by a positive max- 0 0 imal Lyapunov exponent (MLE). In regular regions the 0 0.5 1 1.5 2 2.5 Nonlinearity parameter, κ separationgrowslinearly[16], resulting in zeroMLE.As we increase κ in our system, we expect the phase space FIG. 1: (color online). Averaged finite-time maximal Lya- to change from being dominated by regular regions for punov exponent (FTMLE), λ, and normalized spectral en- small κ to being dominated by chaotic regions for large tropy, η, as functions of the nonlinearity, κ. Ns=21. Inset: κ. In the present section, we use the MLEs to quantify Normalized spectral entropyof finaltime-averaged state ver- this transition to chaos,which, as we will see in the sub- susFTMLEforeachofthe100initialconditionusedtocom- sequentsection, coincides with a relatively broadchange putetheaveraged value for κ=0.36, 0.54, 0.72, 0.9. from unthermalizability to complete thermalizability. Consider two trajectories x(t) and x(t) with initial points x and x , respectively. The separation δx(t) = Thermalizability threshold from spectral entropy.– For 0 0 e x(t) x(t) inietially satisfies a linear differential equa- coupledanharmonicoscillators,asinthe FPUstudy,en- − ergy equipartition among the normal momentum modes teion, and the duration of this linear regime grows with- out bound as the initial separationx x goes to zero. signified thermalization. In the BHM, the additional 0 0 − conservation of the norm modifies the quantity that is The finite-time maximal Lyapunoveexponent (FTMLE) corresponding to the phase-space point x [17] is equipartitioned. To determine the best measure for the 0 equipartitionweusethevariationalHartree-FockHamil- λtfin(x0)=xe0li→mx0 tfi1n lnkxe(tkfixn0)−−xx(0tkfin)k. (3) tHtoioannritaarnlefie[-e2Fld0o].c,kTHehniHseFrpgri=oecse{Pd¯huωnrneH¯hFgω}inHvwFes|aψh¯snωr|e2Hg,Fawr=dheh¯edrωeas+tthh2eeµsveNatriao-f The limit t gives the MLE,eλ (x ). The FTM- n n 0 a− fin →∞ ∞ 0 µ, where µ is the chemical potential. LEsarethemselvesofintrinsicinterestandinthechaotic The Hartree-Fockapproximationis knownto overesti- regime the average over the FTMLE converges to the matetheinteractionenergyintheregimeofstronginter- standardMLE[17,18]. Wechoseaconvenientquantum actions. For this reason, we determine the temperature mechanicalmetric, x x 2 = ψ ψ 2 (see [19]). k − k Pn| n− n| T and the chemical potential using the time-averaged Initially, we studyethe FTMLE onea 21-site lattice for numerical kinetic energy (along with the norm) instead a class of initial conditions where only the k = 0, 1 of the total energy. The temperature and the chemical ± modes are occupied. In this subspace we sample uni- potentialwerecomputedindividuallyforeachinitialcon- formly from the intersection of the microcanonicalshells dition used. in energy and norm; the energy is chosen to be the infi- The new quantity to be equipartitioned is the nite temperature energyof the subsystem, andthe norm distribution of the Hartree-Fock energy, q (t) = n iwse1s.etFoars tehaechinviatilaulepoofinκt,swxe.saTmopeleac1h00inpitoiainltps,oiwnthiwche |ψn(t)|2¯hωnHF/Pn′|ψn′(t)|2¯hωnH′F. A quantitative mea- 0 sure of the distance from thermodynamic equilibrium is addasmallrandomvector,aslittle asmachineprecision the spectral entropy S(t) = q (t)lnq (t), or more allows, to obtain the corresponding x ’s. Each pair we −Pn n n 0 conveniently the normalized spectral entropy [11], propagateforatimet ,shortenoughtoensurelinearity fin e of the evolution of δx(t) but long enough to be able to S S(t) max η(t)= − , (4) clearly distinguish chaotic trajectories from regular ones S S(0) max onaplotoflnδx(t)versust: theformerincreaselinearly, − and the latter, logarithmically [18]. We also verify that where S = lnN is the maximum entropy, which oc- max s the average of the FTMLE’s over the ensemble of initial curs for complete equipartition of q . In Fig. 1 the spec- n conditionsdoesnotdependont aslongasbothcriteria tral entropy of the final time-averaged state, also aver- fin above are satisfied. In Fig. 1 the averaged FTMLEs are aged over 100 initial states (drawn from the same en- plottedasafunctionoftheinteractionstrength. Thereis semble that was used for the Lyapunov exponent calcu- a distinct regime with zeroLyapunovexponentfor small lation) is plotted for each value of κ. For large nonlin- κ < 0.5 and a strongly chaotic regime for κ > 1 where earities, κ > 1, the normalized spectral entropy goes to all∼initial conditions have positive exponent. ∼ zero, indica∼ting remarkable agreement between the final 3 ] 1 J 00..68 κεT==00.0.195 J thienfriimtniaaalll (cid:1)ponent, [ 00..43 NNNsss===421111 0.4 ex 0.2 2ψ(n)| 0. 20 punov 0.1 bution, | 00..68 κεT==00.4.451 J thienfriimtniaaalll (a) Lya 0.00.0 0.N5onli1n.0earity1 .p5aram2.e0ter, (cid:0)2.5 3.0 distri 0.4 1.0 m 0.48 u 0.2 ] Moment 0. 80 κε==11.8.4 J initial [Je, (cid:2)T0.8 00..3462 T final cl 0.6 thermal rti 0.6 0.30 a ] 0.4 er p 0.24 (cid:1) [J p 0.2 y 0.4 g 0.18 0 er -8 -6 -4 -2 0 2 4 6 8 en 0.12 Momentum index, n otal 0.2 0.06 T 0.0 0.00 FIG.2: (coloronline). Initial,finalandHartree-Fockthermal 0.0 0.5 1.0 1.5 2.0 2.5 3.0 momentum distributionsfor κ=0.09,0.45,1.8, startingfrom (b) Nonlinearity parameter, (cid:0) the same initial state. N=21. The initial state is a represen- tative state and the final state is time averaged. ǫ is the FIG. 3: (color online). (a) Averaged Finite Time Lya- T total energy per particle. punov exponent, λ/J, for three different system sizes, Ns = 11,21,41. For each κ, the same energy-per-particle was used for each lattice size. (b) Contour lines of the Lyapunov ex- state and the thermal predictions. Note that this corre- ponent versus the nonlinearity, κ and energy-per-particle, sponds to the chaos threshold observed previously. Fur- ǫT = (H −H0)/Na, where H is the Hamiltonian (1), and thermore,weverifiedthat forlargeκ,the fluctuations in H0=−2J+(1/2)µ0 isthegroundstatevalueofH. Ns =11. The first contour line corresponds to λ = 0.02. The circles c kineticenergyscaleas√N ,confirmingtheirthermalna- s and dotted line give the total energies (per particle) used in ture. Forκ<.5thenormalizedspectralentropyisabove thecalculation for (a). .5 signifying∼that during the time evolution the state of the system remains close to the initial state. As seen in the inset of Fig. 1, an individual initial state with larger excitedmodeinthecaseoftheHamiltonian(1). Assum- FTMLEtendstohavelowerspectralentropy,i.e. torelax ing that the shape of the momentum distribution ψ 2 n | | to a state which is closer to the thermal one. Beginning as a function of n/N should be fixed in the thermody- s at κ 0.5, where the averaged FTMLE is substantially namic limit, the left-hand-side of the above relationship ≈ non-zero,someoftheinitialstatesthermalizecompletely. also remains finite. These observations lead to a conjec- In Fig. 2, the initial and time-averaged momentum ture that the chaos criterion involves only the intensive distributions of a representative state are plotted for parameters and observables, i.e. those that are finite in κ = 0.09,0.36 and 0.9, along with the thermal Hartree- the thermodynamic limit. Fock predictions, ψ 2 =(T/N )/(h¯ω +2µ N µ). Our test for the above conjecture is based on the fact n a n 0 a h| | i − Chaos Threshold for Different Lattice Sizes.– Let us thatforachaoticmotionthemajorityofthe trajectories startfromthenotionthatthe parameterκintroducedin cover the whole available phase space, and as a result (2) is the only dimensionless combination of the param- the MLE becomes, for a given set of parameters, a func- eters of the problem that remain finite in the thermody- tion of just the conserved quantities: energy and norm. namic limit, N , N /N = const,J = const,U = Thisimpliesthatforthesameenergy-per-particle,norm, s a s → ∞ const. Curiously, the chaos threshold for N = 21 is at and nonlinearity parameter κ, the Lyapunov exponents s κ .5, i.e. κ 1. Another observationcomes from a re- for different lattice sizes should be similar. In Fig. 3a ≈ ∼ latedwork[8]onchaosthresholdinNLSEwithhard-wall the averaged FTMLE is plotted for three different lat- boundary conditions. The authors find that the bound- tices, N =11, 21, and 41. For each κ, the same energy- s ary between regular and chaotic motions of momentum per-particle (in units of J) is used for all three lattices. mode, n, is givenby (µ ψ 2)/(h¯ω n) 1, where ¯hω is The corresponding energies are shown by the solid line 0 n 1 1 | | ∼ the lowest excitation energy, e.g. the energy of the first in Fig. 3b. From the plot it is indeed evident that the 4 averagedFTMLE is universal with respect to the size of zero as N−2, i.e. towards the origin in Fig. 3(b), where s the lattice and that the values for N = 11 already give the motion has no dynamical instabilities. s a very good estimate of both the value of the averaged FTMLE and the threshold. We are grateful to Bala Sundaram, Kurt Ja- Two Parametric Theory of the Chaos Threshold.– The cobs, Isabelle Bouchoule, Boris Svistunov, and Anatoli universality observed above suggests the most relevant Polkovnikov for enlightening discussions on the subject. pairofvariablesformappingthechaosthreshold,namely This work was supported by grants from the Office of κ and the total energy-per-particle, ǫ /J. (In FPU one Naval Research (N00014-06-1-0455) and the National T variable is sufficient, ultimately because there is one less Science Foundation (PHY-0621703 and PHY-0754942). conserved quantity.) In Fig. 3b contour lines of the av- eraged FTMLE for N = 11 are plotted versus the non- s linearity parameter and energy-per-particle. We use two setsofinitialconditionswithn=0, 1andn=0, 1, 2 ± ± ± [1] D. Jaksch et al.,Phys. Rev.Lett. 81, 3108 (1998). momentum modes occupied. [2] E. Fermi, J. Pasta, and S. Ulam, Los Alamos Report One can observe a plateau in the averaged FTMLE LA-1940 (1955). for λ < λc = 0.02, given by the solid line. After [3] N. J. Zabusky and M. D. Kruskal, Phys. Rev. Lett. 15, crossing∼the critical line the averaged FTMLE increases 240 (1965). with uniform slope. The critical line resembles a hyper- [4] B.V.Chirikov,At.Energ.6,630(1959)[J.Nucl.Energy, Part C 1, 253 (1960)]. bola with the point of closest approach to the origin at [5] F. M. Izrailev and B. V. Chirikov, Sov. Phys. Dokl. 11, (κ, ǫ ) (0.5, 0.2J), so that the hopping parameter J T ∼ 30 (1966). appears to be a relevant energy scale. This is probably [6] D.Boyanovsky,C.Destri,andH.J.deVega,Phys.Rev. not an accident: for ǫT J the dispersion law ωn be- D 69, 045003 (2004). ≫ ginstodeviatefromthe(quadratic)dispersionlawofthe [7] G. J. L. Gerhardt et al., Chaos Solit. Fract. 13, 1269 integrable NLSE with periodic boundary conditions. (2002). Summary and outlook.– In this paper we consider the [8] P. Villain and M. Lewenstein, Phys. Rev. A 62, 043601 (2000). dynamics ofatomsinanopticallattice fromthe pointof [9] B. M. Herbst and M. J. Ablowitz, Phys. Rev. Lett. 62, viewofchaostheory. Weidentify thethresholdforchaos 2065 (1989). and show that it corresponds to the onset of thermal- [10] M. J. Ablowitz, C. Schober, and B. M. Herbst, Phys. ization. Far above the threshold, the final state of the Rev. Lett.71, 2683 (1993). system is governedby the usual statistical mechanics. [11] R. Livi et al.,Phys. Rev.A 31, 1039 (1985). We seetwopotentialapplicationsofourresults. First, [12] M. Pettini and M. Cerruti-Sola, Phys. Rev. A 44, 975 in quantum nonequilibrium dynamics, our results can (1991). [13] M. Srednicki,Phys. Rev.E. 50, 888 (1994). serve as a guide for identifying the dominant effects [14] M. Horoi, V. Zelevinsky and B.A. Brown, Phys. Rev. preventing thermalization in optical lattices. Based on Lett.74,5194(1995);P.JacquodandD.L.Shepelyansky, the studies of the validity of the classical field the- Phys. Rev. Lett. 79, 1837 (1997); V.V. Flambaum and ory for Bose condensates [21] our results will apply F.M. Izrailev, Phys. Rev. E 55, R13 (1997); M. Rigol, for the lattice site occupations satisfying Na/Ns V. Dunjko and M. Olshanii, Nature (London) 452, 854 ≫ max(κ, 1)max((∆n/N )−1, 1), where ∆n is the typical (2008). s width of the momentum distribution. We note that the [15] M.Greineretal.,Phys.Rev.Lett.87,160405(2001);T. St¨oferle et al., Phys.Rev.Lett. 92, 130403 (2004). Mott regime, N = integer N , ∆n = N , U/J a s s × ≥ [16] B. V.Chirikov, Phys. Rep. 52, 263 (1979). 2.2N /N [22], lies well outside of the above criteria. a s [17] B. Eckhardt and D. Yao, Physica D 65, 100 (1993); N. Second, in chip-based atom interferometry with dense VoglisandG.J.Contopoulos,J.Phys.A27,4899(1994). Bose condensates[23], our results illustrate the fact that [18] G.Contopoulos,L.Galgani,andA.Giorgilli, Phys.Rev. nonlinear instabilities cannot affect the performance of A 18, 1183 (1978). interferometric schemes. Recall that the force fields [19] ThismetricbecomesEuclidianunderthecanonicaltrans- used in interferometry are usually periodic with a pe- formation Qn = (2h¯)1/2Re(ψn), Pn = (2h¯)1/2Im(ψn): riod L = λ/2, where λ = 2π/k, and k is the wavevector kxe−xk2 = 21h¯{Pn+=(N−s(−N1s)−/12)/2[(Q′n−Qn)2+(Pn′−Pn)2]}. of light used to generate the interferometric elements. [20] Y. Castin, in Coherent Atomic Matter Waves, edited by For spatially uniform initial conditions, the time evolu- R.Kaiser,C.Westbrook,andF.David(Springer,Berlin, tioncanbedescribedbyaNLSEwithperiodicboundary 2001), p.32; P.OhbergandS.Stenholm,J.Phys.B 30, 2749 (1997). conditions. Inturn,theNLSEconstitutesthecontinuum [21] Y. Castin, J. Phys.IV France 116, 89 (2004); Y.Castin limitofourmodel,N ,wherewekeepconstantthe s →∞ et al., J. Mod. Opt. 47, 2671 (2000); Y. Kagan and B. ground-statechemicalpotentialµ0, the size ofsystemL, V. Svistunov,Phys.Rev.Lett. 79, 3331 (1997). andtheratiobetweentheenergy-per-particleET andthe [22] C.J. Hamer and J.B. Kogut, Phys. Rev. B 20, 3859 so-calledrecoilenergyE ¯h2k2/2m=π2J/N . Inthis (1979). R s limit both the parameter≡κ and the ǫ /J ratio tend to [23] Y. Wang et al.,Phys. Rev.Lett. 94, 090405 (2005). T