Relaxation of a one-dimensional Mott insulator after an interaction quench

# Relaxation of a one-dimensional Mott insulator after an interaction quench

Marcus Kollar and Martin Eckstein Theoretical Physics III, Center for Electronic Correlations and Magnetism, Institute of Physics, University of Augsburg, 86135 Augsburg, Germany
April 28, 2008
###### Abstract

We obtain the exact time evolution for the one-dimensional integrable fermionic Hubbard model after a sudden change of its interaction parameter, starting from either a metallic or a Mott-insulating eigenstate. In all cases the system relaxes to a new steady state, showing that the presence of the Mott gap does not inhibit relaxation. The properties of the final state are described by a generalized Gibbs ensemble. We discuss under which conditions such ensembles provide the correct statistical description of isolated integrable systems in general. We find that generalized Gibbs ensembles do predict the properties of the steady state correctly, provided that the observables or initial states are sufficiently uncorrelated in terms of the constants of motion.

###### pacs:
03.75.Ss, 05.30.Fk, 71.27.+a, 02.30.Ik

## I Introduction

Recent experiments with ultracold atomic gases Bloch2007pre (); Greiner2002b (); Kinoshita2006 (); Hofferberth07a () have made it possible to study the time evolution of a tunable quantum many-body system that is kept in excellent isolation from the environment. For example, such a quantum system can be forced out of equilibrium by suddenly changing a parameter in the Hamiltonian. Then the system may or may not relax to a new steady state, which is not necessarily the thermal state predicted by statistical mechanics. After such a “quantum quench” the system evolves according to Schrödinger’s equation

 |Ψ(t)⟩=exp(−iHt/ℏ)|Ψ(0)⟩, (1)

where is the prepared initial state and is the new Hamiltonian for times . This situation has recently been studied by a variety of numerical and analytical techniques Altman2002a (); Polkovnikov2002a (); Sengupta04 (); Cazalilla06 (); Rigol06+07 (); Kollath07 (); Manmana07 (); Cramer2008 (); Eckstein07 (); Rigol2008a (); Gangardt2008 (); Barthel2008 (); Moeckel08a (); Anders2005 ().

Due to the unitary time evolution the wave function of an isolated system remains pure for all times and does not converge for . Only the state of a finite subsystem, for which the rest of the system effectively acts as a reservoir, can become stationary Cramer2008 (); Barthel2008 (). Nevertheless, also for the entire system we expect relaxation of the expectation value of an observable to a stationary value for large times. However, this global relaxation can happen only for (i) sufficiently large systems, (ii) sufficiently simple observables, and (iii) sufficiently complicated Hamiltonians, for the following reasons. First of all, (i) many degrees of freedom are needed, so that the thermodynamic limit may be taken, otherwise one expects finite recurrence times Montroll1961 () (see Rigol2005a () for a recent example). Furthermore, (ii) the expectation value of a complicated observable need not relax; for example, the expectation value of , involving the projectors onto two eigenstates of with different energies, oscillates for all times. Usually such projectors are highly nonlocal and their expectation values correspond to correlation functions of very high order. On the other hand, local and few-particle observables are usually simple enough to relax to new stationary values. Finally, (iii) the Hamiltonian that governs the dynamics must also be sufficiently complicated. For example, the magnetization of an Ising chain in a transverse magnetic field relaxes for long-range interactions, but keeps oscillating when only next neighbors are coupled Heims65 (). Similar “collapse-and-revival” oscillations of a many-body system were recently observed in experiments with ultracold atoms by Greiner et al. Greiner2002b (). In their experiments, a Bose condensate was prepared in the potential of an optical lattice which was suddenly steepened Greiner2002b (), then the bosons are essentially only subject to the Hubbard interaction but no hopping between lattice sites occurs. Since in this case has only integer eigenvalues, the wave function will oscillate for all times with period .

For small hopping between lattice sites it follows from perturbation theory that expectation values keep oscillating for short times. But it is not clear what will happen for long observation times. Is relaxation possible for a bosonic or fermionic Hubbard model if the spectrum has a Mott gap, or if the initial state is a Mott insulator? Or is it prevented by the Mott gap in the energy spectrum? The answer is no, not necessarily: in Sec. II we provide an example, the fermionic Hubbard model Gebhard92 (), which shows that relaxation in the presence of a Mott gap is indeed possible. Note that the formation of a fermionic Mott insulator was recently observed with ultracold atoms Jordens2008pre ().

Another central question is whether the steady state of a quenched isolated system can be described by an effective density matrix , such that yields the correct expectation value for any observable which relaxes. Statistical mechanics can be used to make an approximate but usually accurate prediction for this steady-state density matrix. For example, the microcanonical prediction is that const for states with energy close to , and zero otherwise. If indeed agrees with the long-time limit of , we say that the system thermalizes. Clearly, thermalization can be expected only for sufficiently coarse-grained observables; it is always possible to construct a complicated correlation function that depends on the details of the initial conditions and is not described by . As for classical gases, thermalization is generally expected for isolated interacting quantum systems. Indeed, for one-dimensional atomic Bose gases, the dynamics leading to the thermal state were recently observed by Hofferberth et al. Hofferberth07a ().

Thermalization of an isolated system is impossible in certain cases, usually because the system is integrable in the sense that there are infinitely many constants of motion. Simple theoretical examples are one-dimensional integrable models, such as the spin chain, for which the magnetization does not thermalize after a quench of the longitudinal magnetic field Girardeau1969 (); Barouch1969 (); Girardeau1970 (); Barouch1971 (), or the Ising chain in a transverse field, whose correlation functions do not thermalize after a field quench Sengupta04 (). Experimentally, the lack of thermalization was recently observed by Kinoshita et al. Kinoshita2006 () for bosonic atoms confined to one dimension, whose nonthermal stationary momentum distributions were attributed to integrability.

The fundamental assumption of statistical mechanics is that the equilibrium state is characterized only by a few thermodynamic variables such as internal energy and particle number. However if the system is integrable, its infinitely many constants of motion lead to a rather detailed memory of the initial state, because much fewer states are accessible to the dynamics. Nevertheless, even in the absence of thermalization, a statistical prediction for the steady state can be made with a generalized Gibbs ensemble (GGE), as discussed in the recent work of Rigol et al. Rigol06+07 (). For example, if the time evolution of an integrable system is determined by the effective Hamiltonian

 Heff =∑αϵαIα, (2)

where the operators commute, , then the standard choice for the statistical operator of the GGE is constructed from these constants of motion according to Rigol06+07 ()

 ρG =e−∑αλαIαZG, ZG =Tr[e−∑αλαIα], (3)

This choice maximizes the entropy ( ) for fixed expectation values , which are set to their initial-state expectation values by an appropriate choice for the Lagrange multipliers  Balian91a (). The GGE prediction for the steady-state expectation value of an observable is then . GGEs successfully predict some properties of the nonthermal steady states occurring after quenches in integrable or highly constrained systems. For example, they yield the correct nonthermal momentum distribution of one-dimensional hard-core bosons Rigol06+07 () (experimentally realized in Ref. Kinoshita2006, ), for the one-dimensional Tomonaga-Luttinger model Cazalilla06 (), and for the Falicov-Kimball model in infinite dimensions Eckstein07 (). The GGE also yields the correct double occupation for the fermionic Hubbard model, as shown in Sec. II below.

However, in some cases, the stationary values of some observables differ from the statistical predictions of the GGE. For one-dimensional hard-core bosons the unit-cell averaged one-particle correlation function is not described by the GGE Rigol06+07 (). Furthermore Gangardt and Pustilnik Gangardt2008 () pointed out that the GGE (3) may not capture correlations between the conserved quantities . As a consequence the merit of generalized Gibbs ensembles is currently somewhat controversial. In this situation rather general criteria for the validity of GGEs should be useful, which we derive in Sec. III. Our approach is complementary to the work of Barthel and Schollwöck Barthel2008 (), who recently showed that for finite subsystems the reduced density matrix converges to the GGE under certain mathematical conditions on the initial state and Hamiltonian.

It should be noted that an important ambiguity lingers in the construction of the GGE for the Hamiltonian (2). While all the are conserved, so are all their combinations, i.e., all products of the form , , …, leading to the question of whether all such products should also be included in the exponent of the density matrix (3). In Sec. III we show that for observables and initial states which involve little or no correlation, it suffices to fix the constraints , as in Eq. (3), but not their products. We also provide an example for which different choices of the GGE lead to different predictions.

Finally, we note that thermalization might also be prevented in nonintegrable systems due to many-body effects, e.g., the presence of a Mott gap. Nonthermal steady states in nonintegrable systems were observed and argued for in recent numerical studies for finite one-dimensional soft-core bosons Kollath07 () and spinless fermions Manmana07 (). By contrast, thermalization was observed for hard-core bosons on a two-dimensional lattice Rigol2008a (). Fast relaxation to a nonthermal quasisteady state, so-called prethermalization Berges2004a (), was recently observed for the fermionic Hubbard model in high dimensions Moeckel08a (). Further studies of relaxation in nonintegrable many-body systems are therefore desirable, but will not be the subject of this paper.

Our goals in the present paper are thus two-fold. (i) On the one hand, we provide an explicit example of relaxation in a fermionic Mott insulator: In Sec. II we obtain the exact time evolution of the one-dimensional fermionic Hubbard model with hopping and repulsive interaction  Gebhard92 (), starting from the metallic ground state at or the insulating ground state at . We find that the expectation value of the double occupation relaxes with algebraically damped oscillations to a new stationary value for all , i.e., relaxation is not inhibited by the presence of a Mott gap. The long-time limit differs from the thermal value, as expected for an integrable system, but is described by an appropriate GGE. (ii) On the other hand we discuss under which circumstances GGEs describe the steady state of integrable systems after relaxation. We show in Sec. III that their validity depends on the observable, on correlations in the initial state, and possibly the system size.

## Ii Interaction quench in a fermionic Hubbard model

### 1/r Hubbard chain

We consider sudden changes in the interaction parameter of the one-dimensional fermionic Hubbard model,

 H =H0+H1,      H1=U∑ini↑ni↓, (4a) H0 =∑i,j=1..Lσ=↑,↓tijc†iσcjσ=∑|k|<πσ=↑,↓ϵkc†kσckσ, (4b)

with repulsive on-site interaction , bandwidth , and dispersion =, which corresponds to hopping amplitudes that decay proportionally to inverse distance. This lattice model was introduced by Gebhard and Ruckenstein Gebhard92 () as a parent system of the Haldane-Shastry Heisenberg chain Haldane88 (); Shastry88 (), to which it reduces in the limit of large for a half-filled band with density . We consider only ; larger densities can be treated by means of a particle-hole transformation,  Gebhard97 (). For and any number of lattice sites , the model (4) is represented by an effective noninteracting bosonic Hamiltonian, from which ground-state and thermodynamic properties can be obtained analytically Gebhard92 (); Gebhard94b (); Gebhard97 (). For the ground state of (4) is of course the Fermi sea,

 |ψ0⟩ =∏k

with particle density . For , on the other hand, the ground state is Gebhard92 (); Gebhard94a ()

 |ψ∞⟩ =∏i(1−ni↑ni↓)|ψ0⟩, (6)

i.e., the Fermi sea with all doubly occupied sites projected out. At half-filling a Mott-Hubbard metal-insulator transition occurs at interaction strength , with the Mott gap given by for  Gebhard92 (). This metal-insulator transition is also captured by correlated variational wave functions Gebhard94a (); Dzierzawa95 ().

### Interaction quenches

We now consider the following nonequilibrium situation. For times the system is prepared in the ground state for interaction parameter or , i.e.,

 |Ψ(0)⟩ =|ψ0⟩√⟨ψ0|ψ0⟩ metallic state, or (7a) |Ψ(0)⟩ =|ψ∞⟩√⟨ψ∞|ψ∞⟩ Mott insulator. (7b)

Then at time the interaction is suddenly switched to a new value , so that the time evolution for is governed by the Hamiltonian (4), i.e., the system evolves according to Eq. (1). We refer to these two types of quenches as    (starting from the metallic state (7a)) and    (starting from the Mott-insulating state (7b)), respectively. More general initial states corresponding to intermediate values of can also be used; they lead to similar results which are omitted here.

### Bosonic representation

In the bosonic representation of Ref. Gebhard94a, the initial states (7) factorize. They can be written in terms of hard-core bosons () in the form Gebhard92 (); Gebhard94a ()

 |ψ0⟩ (8a) ≡∏K′|ψ0K,K+ΔK⟩, (8b) |ψ∞⟩ =∏K′(1−DK,K+ΔK) |ψ0K,K+ΔK⟩, (8c)

where , , and the prime indicates that only every other bosonic pseudomomentum appears. The Hamiltonian (4) only acts separately on each space spanned by the bracketed configurations and for neighboring pseudomomenta, and  Gebhard92 (); Gebhard94a (),

 Heff =∑K

with . In this representation it is then straightforward, although tedious, to obtain the propagator , its action on , and the expectation value of the double occupation . We take the thermodynamic limit, with fixed density ; the sums over are then replaced by integrals which can be evaluated analytically.

### Results for the double occupation

Setting the bandwidth to (and also 1), our results for the interaction quenches    and    can be written as

 d(t)∣∣0→U∞ =c++f(t), (10a) d(t)∣∣∞→U0 =c−−f(t)U, (10b)

for the two types of quenches, with the abbreviations

 c± =n28∓Δ232U2[2nU+Ω2lnωΩ], (11) f(t) =g(Ω,t)−g(ω,t)−n8ωtsin(ωt)+3cos(ωt)Ut2. (12)

These expressions involve several energy scales, apart from the interaction and the bandwidth ( ), namely the Mott gap , the total bandwidth of the spectrum , and , a characteristic density-dependent energy scale appearing in the holon and spinon excitation energies Gebhard92 (); Gebhard94b (). As functions of the constants have the remarkable symmetry that both are invariant under the replacement (for all ). The function in Eq. (12) is given by

 g(η,t)=132U2[−Ω2Δ2Ci(ηt)+6−Δ2t2t3ηsin(ηt) +6−(Δ2+8U)t2t4cos(ηt)], (13)

where is the integral cosine. For the quench to ( ) this reduces to and , where

 f1(t) =−2t2+316t4+38t3sin(2t)−4t2−316t4cos(2t). (14)

We note that in all cases relaxes with damped oscillations from its initial value ( for the metallic state (7a) or 0 for the Mott insulator (7b)) to a new stationary value. This long-time limit, , always exists, even when quenching to , , or . For these final values of we find

 limt→∞d(t) =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩n24U: 0→∞n2(3−2n)6U: ∞→0n28U: 0→1 or ∞→1. (15)

For the quench to we note in particular that the stationary value of the double occupation differs from the thermal value, . This is discussed in detail at the end of Sec. III.

For general and the function (12) behaves asymptotically as

 f(t)=−n(1−n)2sin(ωt)ωt−cos(Ωt)4Ut2 +(1−3n)ω2+(1−n)Ω22ω2cos(ωt)4Ut2+O(1t3) (16)

for large times, which at half-filling reduces to perfect beating,

 c± =18∓[(1−U)216U+(1−U2)216U2ln∣∣1−U1+U∣∣], (17) f(t) =−cos(Ωt)+cos(Δt)4Ut2=−cos(Ut)cos(t)2Ut2. (18)

We conclude that the relaxation of involves the frequencies and and that it falls off rather slowly as for densities , or as for . This type of algebraic decay is typical for one-dimensional systems Calabrese06a (). Mathematically it can be traced to the Riemann-Lebesgue lemma, i.e., the decay of one-dimensional integrals over oscillating functions Barthel2008 (); Olver1997 ().

The results for for several quenches are shown in Fig. 1,

together with the long-time limit (dotted blue lines) and the thermodynamic prediction (solid red lines). The latter is determined from the exact grand-canonical potential  Gebhard92 (); Gebhard94b (); Gebhard97 (), using the temperature and chemical potential that correspond to the same internal energy and density as the final state. The density is given by , which yields the chemical potential by inversion, and the internal energy per site is . The temperature is obtained as the solution of , where the energy in the final state is for the metallic state (7a) and for the insulating state (7b). Finally, the thermal value of the double occupation per site is obtained as , evaluated at the determined temperature and chemical potential .

The long-time limit of clearly differs from the thermodynamic prediction (see Fig. 1). This is the expected behavior for an integrable system with an infinite number of constants of motion. They act as constraints on the accessible states and must be taken into account into the statistical description of the steady state. As discussed in the Introduction, this can be achieved by employing a GGE Rigol06+07 (). For the Hubbard chain (4) the effective bosonic Hamiltonian is indeed of the form (2), where labels the pairs of bosonic pseudomomenta (,) and the constants of motion are the projectors onto one of the two eigenstates of [Eq. 9]. For the initial states (7) of the form the statistical operator of the GGE (3) becomes

 ρG =∏K′2∑ν=1|⟨νK,K+Δ|ψK,K+Δ⟩|2|νK,K+Δ⟩⟨νK,K+Δ|, (19)

from which we obtain as and for the two types of quenches. The double occupation that is predicted by the GGE thus agrees precisely with the long-time limit [Eq. (11)], i.e., . In the next section we discuss general rules when GGEs are valid, which will explain in particular why the GGE gives the correct stationary value of the double occupation in the Hubbard model.

## Iii Validity of generalized Gibbs ensembles

### Long-time average and diagonal ensemble

In order to address the validity of GGEs we first need to determine the long-time limit of a quantum system after an arbitrary quench. Suppose that an isolated system is prepared at time in an initial state which is described by the density matrix , while the time evolution for is governed by an arbitrary time-independent Hamiltonian . For an initial pure state (e.g., as in Eq. (7)) the density matrix has the form , whereas for a statistical mixture of orthogonal states with probabilities (e.g., for an initial grand-canonical ensemble, as in Ref. Eckstein07, ).

For the time evolution of the density matrix is given by

 ρ(t) =eiHtρ0e−iHt, (20)

and the expectation value of an observable is

 ⟨O⟩t =Tr[Oρ(t)] =∑nn′gg′e−i(En−En′)t⟨ng|O|n′g′⟩⟨n′g′|ρ0|ng⟩, (21)

where are the eigenstates of with energies , and labels possible degeneracies. If the long-time limit exists, then it is necessarily equal to the long-time average ,

 ¯¯¯¯¯¯¯¯⟨O⟩ =limT→∞1TT∫0dt⟨O⟩t =∑ngg′⟨ng|O|ng′⟩⟨ng′|ρ0|ng⟩, (22)

assuming that the limit can be taken termwise (which is allowed for the large but finite systems that we have in mind). In a steady state the system is thus described by the “diagonal ensemble” vonNeumann1929 (); Deutsch1991a (); Brody2007 (); Rigol2008a (),

 ρdiag =∑ngg′|ng⟩⟨ng|ρ0|ng′⟩⟨ng′|=∑nPnρ0Pn. (23)

Here is the projector onto the subspace spanned by the eigenvectors corresponding to the energy eigenvalue . The statistical operator correctly describes the long-time limit, if it exists, of any observable , i.e., .

The diagonal ensemble correctly yields any stationary expectation value, regardless of the transient behavior. However, not a lot is gained by regarding as a “statistical” prediction for the steady state, because each energy eigenstate contributes according to the initial conditions given by . For a nonintegrable system one expects, by contrast, that only a few conserved quantities such as energy and particle number need to be fixed for a successful statistical description in terms of thermal Gibbs ensembles, which can emerge from the diagonal ensemble by means of “eigenstate thermalization” Deutsch1991a (); Srednicki1994a (); Rigol2008a (). Similarly one can ask when a GGE (3) for an integrable system yields the same prediction as the diagonal ensemble (23). This is discussed in the next subsection.

### Gibbs ensemble for an integrable system

We now consider an integrable system whose time evolution after the quench is governed by the effective Hamiltonian (2). For simplicity we consider two typical cases of Hamiltonians only. Either (a) the constants of motion have the eigenvalues 0 and 1 and can thus be represented by fermions or hard-core bosons, , with , ; or (b) the have the eigenvalues and can be represented by bosons, , with . Examples for case (a) are the effective Hamiltonians for hard-core bosons in one dimension Rigol06+07 (), free fermions with quenched disorder Eckstein07 (), and the fermionic Hubbard chain (Sec. II), whereas case (b) applies to the Luttinger model Cazalilla06 (). For both cases the Lagrange multipliers in Eq. (3) are then given by (a) and  (b) .

The Hamiltonian (2) has the eigenstates with occupation numbers and energy eigenvalues . For simplicity we assume that the degeneracy of energy eigenvalues is irrelevant, i.e., the observable or the initial-state density matrix are diagonal in the subspace of eigenvectors with the same energy. This assumption will be examined in detail at the end of this section. From Eq. (22) the diagonal ensemble and long-time average are then given by

 ¯¯¯¯¯¯¯¯⟨O⟩ =∑m⟨m|O|m⟩⟨m|ρ0|m⟩. (24)

Is this the steady-state value predicted by the GGE (3)? We answer this question for two types of observables: for case (a) we consider the observable

 A=∑α1⋯αmβ1⋯βmAα1⋯αmβ1⋯βma†α1⋯a†αmaβm⋯aβ1, (25)

while for case (b) we allow for powers of the bosonic operators and consider (for )

 B=∑α1⋯αm,r1⋯rmβ1⋯βm,s1⋯smBα1⋯αm,r1⋯rmβ1⋯βm,s1⋯sm(b†α1)r1⋯(b†αm)rm(bβm)sm⋯(bβ1)s1, (26)

We assume without loss of generality that vanishes whenever two indices or two indices are the same.

It is straightforward to obtain the long-time average (24) and GGE average (3) of the observables and by using the occupation number basis and the fixed GGE averages . In case (a) we find

 ¯¯¯¯¯¯¯¯⟨A⟩ =∑α1⋯αm˜Aα1⋯αm⟨m∏i=1Iαi⟩0, (27a) ⟨A⟩G =∑α1⋯αm˜Aα1⋯αmm∏i=1⟨Iαi⟩0, (27b)

where we used the identity

 ⟨m∏i=1Iαi⟩G=m∏i=1⟨Iαi⟩G=m∏i=1⟨Iαi⟩0 (28)

in the second line. In case (b) we have

 ¯¯¯¯¯¯¯¯⟨B⟩ =∑α1⋯αm˜Br1⋯rmα1⋯αm⟨m∏i=1ri−1∏k=0(Iαi−k)⟩0, (29a) ⟨B⟩G =∑α1⋯αm˜Br1⋯rmα1⋯αmm∏i=1[ri!(⟨Iαi⟩0)ri], (29b)

where we used the bosonic operator identity

 (b†αi)ri(bαi)ri=ri−1∏k=0(Iαi−k) (30)

in the first line, and the identity

 ⟨m∏i=1 (b†αi)ri(bαi)ri⟩G=m∏i=1⟨(b†αi)ri(bαi)ri⟩G =m∏i=1⟨ri−1∏k=0(Iαi−k)⟩G=m∏i=1[ri!(⟨Iαi⟩G)ri] =m∏i=1[ri!(⟨Iαi⟩0)ri] (31)

in the second line. Furthermore we defined the permutation-averaged matrix elements and .

From these results we obtain rather general sufficient conditions for the validity of the GGE predictions, namely the factorization of initial-state expectation values of (a) products or (b) polynomials of the constants of motion  as follows:

 If~{}⟨m∏i=1Iαi⟩0=m∏i=1⟨Iαi⟩0 then~{}~{}~{}¯¯¯¯¯¯¯¯⟨A⟩=⟨A⟩G. (32a) If~{}⟨m∏i=1ri−1∏k=0(Iαi−k)⟩0=\makebox[0.0pt][l]$m∏i=1[ri!(⟨Iαi⟩0)ri]$ then~{}~{}~{}¯¯¯¯¯¯¯¯⟨B⟩=⟨B⟩G. (32b)

The and in (32a) and (32b) are those for which and , respectively, are nonzero. The criteria (32) are the central result of this section, and we now discuss their implications in detail.

First we note that for simple observables, which involve at most one factor , the factorizations (32) occur trivially. This is the case, e.g., for the double occupancy in the Hubbard model, explaining why the GGE (19) works in this case. Typical observables of an interacting integrable system, however, are often rather complicated when expressed in terms of the constants of motion . Thus all correlations between the constants of motion must vanish in the initial state in order for (32) to be fulfilled. This seemingly restrictive condition can nevertheless be met, because often the initial state is not strongly correlated in terms of the . Moreover, some correlations among the are allowed provided that their contribution to the sum (27) or (29) is negligible, e.g., if it vanishes in the thermodynamic limit.

For example, one-dimensional hard-core bosons are represented by a free-fermion Hamiltonian. For an alternating potential (as studied in Refs. Rigol06+07, ) the initial state contains correlations only between the fermionic momentum number operators and . One can then show from (27) that the GGE (3) makes correct predictions, up to finite-size corrections, for observables that are restricted to a finite region of real space. Another example is the fermionic Luttinger model, which maps to a free-boson Hamiltonian. For an interaction quench (studied in Ref. Cazalilla06, ) the initial state is a product state with correlations only between the bosonic momentum occupation and .

On the other hand, correlations between constants of motion in the initial state, which remain for all times, cannot be described by the GGE (3), as noted in Ref. Gangardt2008, . However, in interacting integrable systems, such observables usually correspond to complicated many-particle operators in terms of the original microscopic degrees of freedom, and thus are not measurable in practice. As mentioned in the Introduction, the microcanonical Gibbs ensemble faces the same problem: one can always construct fine-grained observables that do not thermalize.

Can a GGE be improved if it does not yield the correct long-time average? The minimal necessary extension of the ensemble depends on the observable in question, but as one sees from Eq. (32), it always suffices to fix not only the expectation values of the constants of motion but also the expectation values of all of their products, i.e., by using

 ˜ρG ∝exp(−∑αλαIα−∑abλabIαIβ−⋯), (33)

where the Lagrange multipliers are chosen to fix all products, . Then it follows immediately that and for the observables considered above. Thus any steady state can be described by a sufficiently extended GGE, i.e., by fixing sufficiently many products of constants of motion.

While fixing all products as in Eq. (33) yields the exact long-time average, this extension of the GGE can hardly be regarded as a statistical description of the steady state (as noted in the previous subsection for the diagonal ensemble), because it uses almost the full information about the initial state. In fact, any nondegenerate Hamiltonian acting on a Hilbert space of dimension has pairwise commuting and linearly independent constants of motion Rigol2008a (); fixing all of them in a GGE recovers the diagonal ensemble Brody2007 (). For a nonintegrable system one can choose, e.g., the projectors onto the eigenstates of  Brody2007 (), or linearly independent integer powers of  Manmana07 (). In practice, however, these extensions of the GGE are as hard to calculate as the long-time average.

### Degenerate energy levels

In the previous subsection we assumed that the degeneracy of energy levels is irrelevant [as defined above Eq. (24)], which allowed us to move from Eqs. (22-23) to (24). Below we provide an example for which this assumption does not hold. In that case the expression (24) for the long-time average cannot be used, and thus neither (27a) nor (29a) are available.

We consider a quench to in a general fermionic Hubbard model. Fermions (with spin ) on a Bravais lattice (with lattice sites) are prepared in a correlated unpolarized initial state with fixed densities The time evolution is governed by the free Hamiltonian

 H =∑ijσVijc†iσcjσ=∑kσϵkc†kσckσ, (34)

where labels the crystal momentum (we suppress the vector notation); periodic boundary conditions are assumed for simplicity. This Hamiltonian of the form (2), with the number operators playing the role of the constants of motion . We are interested in the steady-state expectation value of the double occupation .

Assuming again that the degeneracy of energy levels is irrelevant, we obtain for the long-time average (24), using the basis ,

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨ni↑ni↓⟩ =∑m⟨m|ρ0|m⟩1L2∑kk′mk↑mk′↓ =Tr[ρ0L2∑ijni↑nj↓]=n↑n↓=n24. (35)

The same value is obtained from the canonical and grand-canonical ensemble, and also from the generalized Gibbs ensemble which uses the number operators as constants of motion:

 ⟨ni↑ni↓⟩G =1L2∑kk′⟨nk↑⟩G⟨nk′↓⟩G =1L2∑kk′⟨nk↑⟩0⟨nk′↓⟩0=n↑n↓=n24. (36)

Thus we conclude that the double occupation thermalizes to the value after a quench to in any Hubbard model, provided that the degeneracy of energy levels is indeed irrelevant.

Interestingly, this statement disagrees with our exact results for the Hubbard chain from Sec. II: when quenching from to we obtained the long-time limit as  [Eq. (15)]. This differs from the long-time average (35) because the degeneracy of energy levels is in fact relevant for this quench: the initial-state density matrix does not factorize in the free-fermion basis and the linear dispersion leads to a massive degeneracy for the free-fermion energy eigenstates . Therefore the long-time limit is not given by Eq. (24) and does not equal .

Furthermore this example shows that GGEs based on different representations of the constants of motion can yield different results. The free-fermion GGE (36) predicts the wrong value , whereas the GGE (19) which uses the effective bosonic representation gives the correct value . In fact the choice of constants of motion for the construction of a GGE (3) is always ambiguous, as discussed in the Introduction. Nevertheless it is possible to determine the correct GGE a priori, i.e., without knowing the real-time dynamics, by verifying that the degeneracy of energy levels is irrelevant and that the conditions (32) are fulfilled.

## Iv Conclusion

The exact real-time dynamics of the double occupation in the fermionic Hubbard chain shows that the presence of a Mott gap does not inhibit the relaxation after an interaction quench. Its steady-state properties are correctly predicted by generalized Gibbs ensembles. Furthermore we showed for a general class of integrable quantum systems that the GGE prediction equals the long-time average, provided that the observables or initial states are sufficiently uncorrelated in terms of the constants of motion.

## Acknowledgements

Useful discussions with Marcos Rigol, Stefan Kehrein, Corinna Kollath, Krzysztof Byczuk, and Dieter Vollhardt are gratefully acknowledged. This work was supported in part by the SFB 484 of the Deutsche Forschungsgemeinschaft.

## References

• (1) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 85 (2008).
• (2) M. Greiner, O. Mandel, T. W. Hänsch, and I. Bloch, Nature 419, 51 (2002).
• (3) T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
• (4) S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm, and J. Schmiedmayer, Nature 449, 324 (2007).
• (5) E. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 (2002).
• (6) A. Polkovnikov, S. Sachdev, and S. M. Girvin, Phys. Rev. A 66, 053607 (2002).
• (7) K. Sengupta, S. Powell, and S. Sachdev, Phys. Rev. A 69, 053616 (2004).
• (8) M. A. Cazalilla, Phys. Rev. Lett. 97, 156403 (2006).
• (9) M. Rigol, A. Muramatsu, and M. Olshanii, Phys. Rev. A 74, 053616 (2006); M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
• (10) C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
• (11) S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, Phys. Rev. Lett. 98, 210405 (2007).
• (12) M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev. Lett. 100, 030602 (2008).
• (13) M. Eckstein and M. Kollar, Phys. Rev. Lett. 100, 120404 (2008).
• (14) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
• (15) D. M. Gangardt and M. Pustilnik, Phys. Rev. A 77, 041604(R) (2008).
• (16) T. Barthel and U. Schollwöck, Phys. Rev. Lett. 100, 100601 (2008).
• (17) M. Moeckel and S. Kehrein, Phys. Rev. Lett. 100, 175702 (2008).
• (18) F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 (2005).
• (19) E. Montroll, in Lectures in Theoretical Physics, edited by W. E. Brittin, B. W. Downs, and J. Downs (Interscience, New York, 1961), Vol. III, p. 221.
• (20) M. Rigol, V. Rousseau, R. T. Scalettar, and R. R. P. Singh, Phys. Rev. Lett. 95, 110402 (2005).
• (21) S. P. Heims, Am. J. Phys. 33, 722 (1965).
• (22) F. Gebhard and A. E. Ruckenstein, Phys. Rev. Lett. 68, 244 (1992).
• (23) R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, eprint arXiv:0804.4009.
• (24) M. D. Girardeau, Physics Letters A 30, 442 (1969).
• (25) E. Barouch and M. Dresden, Physical Review Letters 23, 114 (1969).
• (26) M. D. Girardeau, Physics Letters A 32, 67 (1970).
• (27) E. Barouch and B. McCoy, Physical Review A 3, 786 (1971).
• (28) R. Balian, From Microphysics to Macrophysics: Methods and Applications of Statistical Physics, vol. 1 (Springer, Berlin, 1991).
• (29) J. Berges, S. Borsányi, and C. Wetterich, Phys. Rev. Lett. 93, 142002 (2004).
• (30) F. D. M. Haldane, Phys. Rev. Lett. 60, 635 (1988).
• (31) B. S. Shastry, Phys. Rev. Lett. 60, 639 (1988).
• (32) F. Gebhard, The Mott metal-insulator transition: models and methods (Springer, Berlin, 1997).
• (33) F. Gebhard, A. Girndt, and A. E. Ruckenstein, Phys. Rev. B 49, 10926 (1994).
• (34) F. Gebhard and A. Girndt, Z. Phys. B 93, 455 (1994).
• (35) M. Dzierzawa, D. Baeriswyl, and M. Di Stasio, Phys. Rev. B 51, 1993 (1995).
• (36) P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 (2006).
• (37) F. W. J. Olver, Asymptotics and Special Functions (AK Peters, Ltd., Wellesley, MA, USA, 1997).
• (38) J. von Neumann, Z. Phys. 57, 30 (1929).
• (39) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
• (40) D. C. Brody, D. W. Hook, and L. P. Hughston, J. Phys. A: Math. Theor. 40, F503 (2007).
• (41) M. Srednicki, Phys. Rev. E 50, 888 (1994).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters