Dynamics of thermalisation in small Hubbard-model systems S. Genway,1 A. F. Ho,2 and D. K. K. Lee1 1Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom 2Department of Physics, Royal Holloway University of London, Egham, Surrey TW20 0EX, United Kingdom (Dated: January 24, 2011) We study numerically the thermalisation and temporal evolution of a two-site subsystem of a fermionic Hubbard model prepared far from equilibrium at a definite energy. Even for very small systems near quantum degeneracy, the subsystem can reach a steady state resembling equilibrium. Thisoccursforanon-perturbativecouplingbetweenthesubsystemandtherestofthelatticewhere relaxation to equilibrium is Gaussian in time, in sharp contrast to perturbative results. We find similar results for random couplings, suggesting such behaviouris generic for small systems. 1 1 0 Understandingthe originofstatisticalmechanicsfrom s (b )ofenergyε (ǫ ),actssolelyonthesubsystem 2 | iS | iB s b a purely quantum-mechanical description is an interest- (bath). λV is the coupling between the subsystem and n ing areaof activeresearch[1–9]. Ofparticularinterest is the bath. For λ=0,the eigenstatesareproductsofsub- a the situation of an isolated quantum system partitioned system and bath eigenstates, denoted sb , with energies J | i into a subsystem and a bath. We ask the question: how E = ε + ǫ . The homogeneous case corresponds to 1 sb s b doobservablesonthesubsystemthermalisewhentheto- λ=1. Specifically,we choosea two-sitesubsysteminan 2 tal system is in a pure state? Seminal works [1, 2] have L-site Hubbard ring of fermions: ] demonstrated the concept of ‘canonical typicality’ that h H = J (c† c +h.c.)+U(n n +n n ), p most random pure states of well-defined energy for the S − X σ 1σ 2σ 1↑ 1↓ 2↑ 2↓ - totalsystemleadtothermalisedreduceddensitymatrices σ=↑,↓ nt (RDMs) for the small subsystem. Numerical works have L−1 L H = J (c† c +h.c.)+U n n , a demonstratedthermalisationinspinorbosonsystemsfor B −X X σ iσ i+1,σ X i↑ i↓ u variousobservablesofthesubsystem[3–5](andoftheen- i=3 σ=↑,↓ i=3 q tire closed system [6, 10]). Recent theoretical work [11] V = J (c† c +c† c )+h.c. , (1) [ has investigatedwhether thermalisation of small subsys- − X σh 2σ 3σ 1σ Lσ i σ=↑,↓ 2 tems, initially far from equilibrium, is generic. v where n =c† c is the number operatoron site i with 3 In this letter, we investigate the temporal relaxation iσ iσ iσ spin σ. The lattice is a ring with the subsystem sites 6 towardsasteadystate,focusingonthe regimewherethe at i = 1,2 and bath sites at i = 3 to L. The hop- 8 steady state appears thermalised. We consider a small 4 ping integrals are Jσ =J(1+ξ sgn(σ)), with a non-zero Hubbard ring of fermions away from half filling, with . ξ =0.05toremoveleveldegeneraciesduetospinrotation 5 two adjacent sites as a subsystem and the other sites as symmetry. We set the on-site repulsion U = J to give 0 abath. Wepreparethesysteminaproductstateofsub- 0 us a metallic system with interacting bath modes while systemandbathpure states ina narrowenergywindow. 1 avoiding the formation of strong features in the many- : Evenforsuchasmallsystem,wefindasteady-stateRDM body density of states at U J arising from Hubbard v close to a thermal state, down to quantum degenerate ≫ i interactions. We will let J be the unit of energy. This X temperatures. Moreover,wefindthattheRDMdiagonal Hamiltonianpreservesthe totalparticlenumber, N,and r elementsapproachasteadystateasanexponentialdecay spin component, Sz, but not the total spin S2. In this a forweaksubsystem-bathcoupling. ThisbecomesaGaus- work, we fill the L = 9 lattice with 8 fermions of total sian decay at a non-perturbative coupling, with a decay spin Sz = 0. The two-site subsystem has 16 eigenstates ratethatdepartssignificantlyfromtheFermigoldenrule. and the bath has 8281 eigenstates, while the compos- WenotethatthisisdistinctfromtheGaussianbehaviour ite system has a total of 15876 states with average level in driven systems that remain out of equilibrium [12], spacing ∆ 10−3. This is small enough to allow exact and in decoherence dynamics [13] of off-diagonal RDM ≃ diagonalisation, but large enough to provide a smooth elements in systems that cannot thermalise. density of states. The Model. Taking motivation from cold atoms in op- Considerasystempreparedinapurestateoftheform tical lattices [14, 15] where local addressing is possible, we study a local cluster in a generic (non-integrable) bu 1 Ψ(t=0),E = s b (2) interacting system with a quasi-continuous spectrum. | 0i X √B| i ii bi=bl We will examine how a local subsystem (S) thermalises with the rest of the system as a bath (B) via unitary where s is the initial subsystem state, e.g. , i S S | i | ↑ ↑i evolution of the whole system under the Hamiltonian withparallelspinsonthetwosites. Ψ(t=0) containsa | i H = H +H +λV where H (H ), with eigenstates linear combinationof B bath eigenstates b within an S B S B i B | i 2 energyshellofwidthδ ,chosensuchthat ΨH Ψ =E . value λ 1. The RDM becomes virtually diagonal — B 0 th h | | i ≪ Thewidthδ (=0.5inthiswork)issmallonthescalefor even the sum over the fluctuating off-diagonal elements, B variationsinthe densityofstates. The systemevolvesin [ r2 ]1/2, is 10−1 to 10−3 smaller than each diago- time: Ψ(t) = e−iHt Ψ(0) . The subsystem is described nPalse6=lse′msesn′t. Fig. 1 shows the steady-state values of the | i | i τ completelybytheRDMwhichtracesoverthebathstates diagonal elements of the RDM, r = lim ρ(t)dt/τ τ→∞R0 b : ρ(t)=Tr Ψ(t) Ψ(t). Thisisevaluatedusingthe as a function of the composite energy E for a coupling B B 0 | i | ih | eigenstates of H from exact diagonalisation. of λ = 0.5. For a variety of initial states, ρ(t) is only Equilibrium States. Before discussing relaxation dy- weakly dependent on the details of the initial state at namics during thermalisation, we identify first the pa- long times for 3 > E > 6, approaching the thermal 0 − rameter regime where ρ does relax to thermal equilib- form ω expected from the canonical ensemble. rium. We say that a subsystem thermalises if its RDM ρ(t) approaches the thermal RDM ω after some time t 1.2 20 (shorter than ∆−1.) The thermal RDM, ω, is diagonal 1 15 withelementsω = sω s N (E ε ,N n ,Sz serzsg),ywεsh,enres p|sairStsicissleasSasnhudb|ssyp|sitineSms∝zseaingBdenNstB0a−t(eEs,|snibS,−swzb)istihstehn−e- ∆, rω 00..68 Teff 1 50 number of bath states with energies in [E,E+δE] with σ 0 nb particlesandspinszb. Wehavetospecifyenergy,num- 0.4 0.1 λ 1 ber and spin because they are globally conserved by the 0.2 Hamiltonian H. We can define an effective temperature 0 T =[∂logN /∂E ]−1 providedthatthe systemisin 0.01 0.1 1 eff B |E0 λ a state with energy uncertainty δE ∆, the level spac- ≫ ing. (In the thermodynamic limit, ω takes the form of a FIG. 2. Dependenceon the initial state, ∆r (solid), and dis- Gibbs canonical distribution [10] — if particles are not tance from the thermal state, σω (hollow), as a function of exchanged with the bath, ωss ∝e−εs/T.) coupling λ, at composite energies E0 = −2 (circle) and 1.77 (square). Inset: effective temperature T for E =−2. eff 0 0.3 λ=0.5 0.25 ε Next, we establish the range of the coupling λ over 1 0.2 whichthesystemforgetsitsinitialstateandthermalises. We expect the system to retain memory of the initial rss 0.15 stateatweakcoupling(λ 1). Moreover,forλ 1,the ≪ ≫ eigenspectrum becomes significantly altered by the cou- 0.1 ε pling,splittingintoseveralbandsandweseeoscillations. 2 0.05 ε This is a feature of the projection of the initial state on 3 ε 4 the strongly-coupled link. Therefore, we expect that the 0 lossofmemoryoftheinitialstateandthermalisationare -6 -4 -2 0 2 4 6 8 10 E possible only in a range of intermediate couplings. To 0 quantify this, we calculated the root-mean-square varia- FIG. 1. (Colour online) Time average of diagonal RDM el- tionindiagonalRDMelementsduetousingdifferentini- eomf ceonmtsp,orssiste, efonrertghyeEn0s, f=or2in,itsizsal=sub0syssetcetmorsatsataesf|u↑n↓c,t↑ioiSn ti.a.l.suabvseyrastgeinmgsotvaetresa:ll∆1r6=ini12tiPalss[thartsess2iin−thhresssiu2b]12sy,swteitmh (solid line), |↑,↑iS (dashed) and |↑,↓iS (dotted). The four Fhocki basis (i.e., eigenstates at J =0). A small ∆r indi- elements are labelled by their energies, ascending from ε to 1 cates memory loss. We have also measured the closeness ε . Thick lines: the corresponding elements for the thermal st4ateω,foundbycountingbathstateswith spinSz−szs and tothe thermalstateω usingσω = 21Psh|rss−ωss|i. We N −ns particles within a Gaussian energy window of width see from Fig. 2 that memory loss and thermalisation oc- δE =0.5, centered on energy E0−εs. curinthe intermediaterangeλth >λ>3with crossover value λ 0.1 at E = 2 and 1.77. th 0 ≃ − We also find that the relative probabilities of different We now present our results for a system starting from states in the n =2, sz =0 sector fit a Boltzmann form: s s the initial states (2). We avoid the regime of very small logr = ε /T +const. For states near the centre of ss s eff − subsystem-bath coupling λ where the subsystem RDM, theeigenspectrum(E 1.77),theeffectivetemperature 0 ≃ ρ, is strongly dependent on the initial state even at long T is infinite. At E = 2, we find T 2 up to λ 2 eff 0 eff − ≃ ∼ timesduetofinite-sizeeffects. Nevertheless,wefindthat (Fig. 2 inset). We estimate the chemical potential to be even such a small system can reach a steady state for 2J 2 so that, unlike in previous work, we see thermal- ≃ couplings λ larger than a surprisingly small crossover isation at temperatures down to quantum degeneracy. 3 We note that these thermalised systems are surpris- 100 ingly small. Popescu et al. [2] give an estimate of the 10 number d of composite-system eigenstates spanned by R the initial state sufficient for thermalisation — if the 1 probability that σ > Y = 0.1 is at least as small as ω X = 0.01, then d > (9π3/2Y2)ln(2/X) 70000. This R ≃ Γ 0.1 isalmosttwoordersofmagnitudelargerthanthenumber of states (&δB/∆) spanned by our initial state which is 0.01 aslowas950. Moreover,wefindthermalisationatU =J for smaller systems than at U J. We believe strong 0.001 ≪ λH λR λH λR inelastic scattering in the interacting bath enables effi- exp exp Gauss Gauss 0.0001 cient thermalisationat U =J when system size is larger 0.01 0.1 1 10 100 than the inelastic scattering length ( J2/U2 for small λ ∝ U/J). Time Evolution. Having established the coupling FIG.4. (Colouronline)Ratesofdecayofρsisi forthesubsys- range for thermalisation for model (1), we will now dis- tem state si =|↑,↑i with initial state (2) at E0=−2 ( /+) and1.77((cid:4)/×)andrandommodelatE =−2(#/△). Weak cussourmainresultsforthetemporalrelaxationtowards 0 the steady state. Fig. 3 shows examples of the time evo- coupling/exponential decay: −dρsisi/dt|t=0 at short times foundfrom exponentialfits( ,(cid:4),#)agreewithFermigolden lutionofthediagonalRDMelementρ (t)withs=s for ss i ruleprediction,γ (solid lines, gradient2). Moderatecou- FGR two coupling strengths. The system is again preparedin pling/Gaussian decay: fit to Gaussian with rate Γ (×,+,△) the product state (2) with |siiS =|↑,↑iS. These results agreeswithΓ1 (dashedlines, gradient1). Verticallinesmark are computed for energy E0 = 2. We do not expect estimates of the crossover values λexp and λGauss (λH/R for − our results to depend strongly on E unless the system Hubbard/random models). Data in the crossover region are 0 is close to a strongly correlated ground state. rates obtained from attempted fitsto either form. 1 0.8 λ=0.1 λ=1 in λ, the RDM element corresponding to a subsystem ρ(t)ss 00..46 state s6=si is approximated by ρ˜(t): 0.2 t1 2 0 4λ2 bu sin[(E E )t] 0 80 1t60 240 0 0.5 t1 1.5 2 ρ˜ss(t)= B Xb (cid:12)(cid:12)(cid:12)(cid:12)bXi=bl Esbsb−−Essiibbii 2 hsb|V|sibii(cid:12)(cid:12)(cid:12)(cid:12) (3) (cid:12) (cid:12) FIG.3. TheRDMelementρssasafunctionoftimetwiths= after the composite system is prepared in the state (2). |↑,↑i for the initial state (2) at total energy E0 =−2 (Teff ≃ Theelementρ˜sisi ismostreadilyfoundbyusingTr(ρ˜)= 2) with si = s. Left: coupling λ = 0.1 with exponential fit 1 to give ρ˜ (t) = 1 ρ˜ (t). This perturbation (dashed). Right: λ=1 with Gaussian fit (dashed). sisi −Ps6=si ss theoryisvaliduntiltimet whenρ hasdroppedsignif- 2 sisi icantly below unity. For times between t and t , eq. (3) 1 2 We find qualitatively different relaxation behaviour follows the Fermi golden rule (FGR): ρsisi(t) decreases forperturbativeandnon-perturbativecouplings(Fig.3). linearly in time with dρsisi/dt∝−λ2. Beyond the FGR Whereas the RDM relaxes towards the steady state ex- regime,weexpectto seeexponentialdecay(see,e.g., ap- ponentially in time at weak coupling (λ < λ ), the proximateMarkovianschemes ofthe Lindbladtype [16]) exp relaxation follows a Gaussian form at larger coupling asisfoundinourdata(Fig.3left)forλ>λexp =0.1. In (λ > λ ). Interestingly, this Gaussian regime cov- our case, the initial state is not a bath eigenstate. This Gauss ers the coupling range where the system thermalises. givessmallfluctuations ontopofasimple linear-tdecay, We can understand our results at short times or duetointerferencebetweentermsintheinnersumin(3). weak coupling. At short times, we can approximate We check in Fig. 4 that the FGR prediction agrees |Ψ(t)i ≃ (1 − iHt)|Ψ(0)i. It can be shown (and quantitatively with dρsisi/dt|t=0 for E0 = −2 and 1.77, our numerics agree) that the element ρ (t) 1 found from the parameters obtained for the exponential Γ21t2 for t < t1 = 1/max(Esb − E0)s,iswi ith≃Γ1 −= fit to ρsisi for t > t1: ρsisi(t) ∼ Ae−(t−t0)/τ +(1−A). λfe[rPensc6=esib,be|thwsebe|Vn|sΨta(0te)si|2c]o1u/2p.ledThbeymhoapxpiminugm(Ve)neisrgoyfdthife- TtimheeFtGbeRtwreaetne,tγFaGnRd,tis: foγund bty=avte[rρ˜agin(gt′)(3)1o]dvte′r/ta. 1 2 − FGR R0 sisi − orderofthesingle-particlebandwidth4J andsot 1/4. This procedure is needed for a non-zero level spacing ∆. 1 ≃ Atweakcoupling,wecangobeyondt bytreatingthe We point out that the exponential fit fails at very weak 1 couplingλV asaperturbationtotheuncoupledHamilto- coupling (λ 10−2) when the system barely relaxes. ∼ nianH +H . It is readily shownthat, to leading order We now considerlargercouplings whereλ O(1). In- S B ∼ 4 steadofexponentialrelaxation,wefindagoodfit(Fig.3 subsystem and the bath, by inelastic collisions of the right)to aGaussiandecay: ρ (t) Ce−Γ2t2+(1 C). fermionswithinthistimescale. Thisshouldbeinsensitive sisi ∼ − This is seen for couplings λ ? λ = 1. The decay to system size for systems larger than the inelastic scat- Gauss rate Γ now increases linearly with λ and is as large as tering length. Therefore,we speculate that the observed the bandwidth scale 1/t 4. It appears insensitive to Gaussian relaxation should remain for large systems. 1 ≃ energy E . Interestingly, we see in Fig 4 that, in the 0 regime where Γt1 ∼ 1, the decay rate Γ is well approxi- expdoenceanytial crossover Gdauecssaiyan mated by Γ from the short-time expansion, which sug- λ 1 gestsΓ=Γ1/C21. (Inourdata,0.97<C21 <1.) Inother λth λexp λGauss memory thermalised words,perturbationtheorygivestheearly-timeprecursor to the full Gaussian form. This suggests the interpreta- FIG. 5. Schematic behaviour of the reduced density matrix tionthatthetimeintervalofvalidityoftheFermigolden ρ(t) as a function of subsystem-bath coupling λ. Top: relax- rule (t <t<t ) narrows and vanishes as λ increases to 1 2 ation to steady state. Bottom: steady state of ρ(t). Memory λGauss. For coupling range λexp < λ < λGauss, the be- of initial state also occurs at large λ in the Hubbardcase. haviourislessclearcut—thedecaystartsasaGaussian but becomes exponential at later times. The amplitude of this exponential tail decreases with increasing λ, be- We are grateful to Miguel Cazalilla for useful discus- coming negligible as λ reaches λ . sions. We wish to thank Imperial College HPC for com- Gauss Random Couplings. To verify that the two regimes of putingresourcesaswellasEPSRCforfinancialsupport. relaxation are not specific to our model Hamiltonian, we proceeded to study an alternative model where the subsystem-bathcouplingV isreplacedbyarandomHer- mitian matrix W which still respects the conservationof [1] S.Goldstein,J.L.Lebowitz,R.Tumulka,andN.Zangh`ı, theglobalparticlenumberN andspinSz. Eachnon-zero Phys. Rev.Lett. 96, 050403 (Feb 2006). matrix element of W is Gaussian distributed, with the [2] S.Popescu,A.J.Short,andA.Winter,NaturePhys.2, variance chosen such that Tr(W2) = Tr(V2). Thus, we 754 (2006). cancompareH =H +H +λV withH =H +H +λW [3] R.V.Jensen andR.Shankar,Phys.Rev. Lett.54, 1879 S B S B with similar decay rates. In this model, we expect 1/t (Apr 1985); K. Saito, S. Takesue, and S. Miyashita, J. 1 to be of the order of the full bandwidth 20 for N =8, Phys. Soc. Jpn. 65, 1243 (1996). Sz =0. Wefindexponentialrelaxationat∼weakcoupling, [4] M.J.Henrich,M.Michel,M.Hartmann,G.Mahler,and J. Gemmer, Phys. Rev.E 72, 026104 (Aug2005). λ>λ =0.8,andwe recoverGaussianrelaxationwith exp [5] S.Yuan,M.Katsnelson, andH.DeRaedt,J.Phys.Soc. a linear-λ decay rate for λ ? λGauss = 8 (Fig. 4, hol- Jpn. 78, 094003 (2009). low symbols). The crossovervalues, λexp andλGauss, oc- [6] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 cur at nominally higher couplings than for the Hubbard (2008); M. Rigol, V. Dunjko, V. Yurovsky, and M. Ol- ring (1). Theybecome closerto theHubbard-ringvalues shanii, J. Phys. A 40, F503 (2007). if we mimic the structure of V by restricting the states [7] H. Tasaki, Phys.Rev. Lett.80, 1373 (Feb 1998). coupledbyW: hs′b′|W|sbi6=0onlyif|Es′b′−Esb|<4J, [[89]] PM..RSeriemdnanicnk,i,PJh.yPs.hRyse.vA. L2e9t,t.L9795,(116909460)4. (Oct 2007). the single-particle bandwidth. [10] D. C. Brody, D. W. Hook, and L. P. Hughston, We summarise our results in Fig. 5. We have shown J. Phys.A:Math. Theor. 40, F503 (2007); B. V. Fine, that a two-site subsystem of the Hubbard model relaxes Phys. Rev.E 80, 051130 (Nov 2009). to steady states resembling canonical thermal states, [11] N. Linden, S. Popescu, A. J. Short, and A. Winter, even for systems with a handful of sites and at quantum Phys. Rev.E 79, 061103 (Jun 2009). degenerate energies. This occurs at a non-perturbative [12] D.Segal,D.R.Reichman,andA.J.Millis,Phys. Rev.B 76, 195316 (Nov 2007). coupling between the subsystem and bath, correspond- [13] D. Rossini, T. Calarco, V. Giovannetti, S. Montangero, ing to nearly homogeneous systems. In this regime, the and R. Fazio, Phys. Rev.A 75, 032333 (Mar 2007); reduceddensitymatrixρ(t)displaysGaussianrelaxation F. M. Cucchietti, S. Fernandez-Vidal, and J. P. Paz, to the thermal state, with a decay rate Γ linear in the ibid. 75, 032337 (Mar 2007); V. V. Dobrovitski, A. E. coupling λ. This contrasts sharply with the perturba- Feiguin,D.D.Awschalom,andR.Hanson,Phys. Rev.B tiveregimewhereρ(t)exhibitsanexponentialrelaxation 77, 245212 (Jun 2008); S. Yuan, M. I. Katsnelson, and with a λ2 decay rate. We believe that the Gaussian re- H. DeRaedt, ibid.77, 184301 (May 2008). [14] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, laxation to thermalisation is a generic feature of closed I. Bloch, and S. Kuhr,Nature467, 68 (2010). nanoscalesystems,asissupportedbyourresultsforran- [15] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. dom Hamiltonians. Gillen,S.Folling,L.Pollet,andM.Greiner,Science329, Finally, we note that it can be shown that Γ1t1 547 (2010). ∼ λJt λ irrespective of system size. The subsystem [16] M.-P. Breuer and F. Petruccione, The Theory of Open 1 ∼ thermalises on the time scale of a few hops between the Quantum Systems (Oxford University Press, 2006).