Thermalization in Systems with Bipartite Eigenmode Entanglement

It is analytically shown that the asymptotic correlations in exactly solvable models following a quantum quench can behave essentially as thermal correlations provided the entanglement between two eigenmodes is sufficiently strong. We provide one example and one counter example of this observation. The example illustrates the fact that the thermal correlations arise from initial states where the entanglement between the eigenmodes stems from the existence of a large energy gap in the initial state. On the other hand, the counter-example shows that when the bi-partite entanglement of the eigenmodes stems from interactions that do not open a gap, the correlations at asymptotically long times are non-thermal. We also show that the thermal behavior concerns only the asymptotic correlation functions, as the difference with an actual thermal ensemble can be observed measuring the energy fluctuations of the system. The latter observation implies a breakdown of the fluctuation-dissipation theorem.


I. INTRODUCTION
Experiments with ultracold gases [1][2][3][4][5], and, in particular, the ongoing efforts to build quantum emulators using ultracold atoms loaded in optical lattices, have aroused much interest in understanding the thermalization mechanisms of integrable models [1,[6][7][8][9][10][11][12][13]. The latter can be used as simple systems to validate a quantum emulator by comparing the outcome of the experiment to the exact solution, prior to using the emulator to study other, more complex models, for which no exact solutions are known.
However, in order to understand the outcome of a simulation of an integrable or exactly solvable model it is important to understand the effect of the initial conditions. Since quantum emulators of ultracold atoms are largely isolated systems and its evolution is essentially unitary, it becomes necessary to understand the conditions under which the asymptotic state of the system can be described by a standard statistical (i.e. Gibbsian) ensemble, or, as it was pointed out recently by Rigol and coworkers [6], must be described by a generalized Gibbs ensemble (GGE). The latter captures the fact that the existence of a non-trivial set of integrals of motion strongly constraint the non-equilibrium dynamics of the system. In Ref. [12], we showed that the GGE can be analytically derived for exactly solvable models and a general class of initial states. In particular, we showed that dephasing between different modes makes equal time correlation functions of both local and non-local operators non-ergodic, in the sense that in the thermodynamic limit, their infinite time limit only depends on the occupation of the eigenmodes in the initial state only. Thus, the asymptotic values of those correlation functions can be equivalently obtained by assuming that the correlations with other eigenmodes produce an effective temperature. This yields a description of the asymptotic correlations that is entirely equivalent to the GGE.
Nevertheless, it was also noticed in Refs. [7,8] that certain kinds of initial states can lead to asymptotic values of the observables that are essentially indistinguishable from those computed with a standard thermal Gibbs ensemble. Other cases of (pre-) thermalization have been found in integrable [14] and non integrable [15][16][17] systems for particular classes of initial conditions. Recently, Mitra and Giamarchi [18] also showed that the adiabatic introduction of a non-linearity following a quantum quench in the Luttinger model, can lead to thermalization as described by the standard Gibbs ensemble. The authors of Ref. [18] emphasized the importance of "mode coupling" for thermalization. In this work, we show that for certain classes of quenches for which two sets of modes are strongly entangled in the initial state, the GGE ensemble can be arbitrarily close to the standard Gibbs ensemble. Using the methods of Ref. [12], we find a simple instance of the mechanism by which initial states can lead to correlations that essentially look thermal. As we show below, this happens when the eigenmodes of the Hamiltonian that describes the evolution of a system following a quantum quench have a certain kind of bi-partite quantum correlations (i.e. entanglement) in the initial state. As an application, we find an analytical explanation for the numerical observations first reported in Ref. [7]. This fact could be used as a simple protocol to produce asymptotic correlations that essentially look as (rather high temperature) thermal correlations in exactly solvable systems. Furthermore, we will show below that the effective temperature that characterizes the asymptotic thermal correlations can be related to the entanglement spectrum of a subset of entangled modes, which means that the latter is accessible experimentally by measuring the effective temperature that characterizes the correlations at long times after the quench.
Before illustrating the above point with one example and one counter example, let us first describe the general set up which will be addressed below. Consider a system containing two subsystems A and B that are initially coupled together. For times t ≤ 0, the system is described by a Hamiltonian of the form H 0 = H A + H B + H AB where H A ,H B and H AB are quadratic in some eigenmodes {a k , b k } which carry a quantum number k (which forms an continuum in the thermodynamic limit) and can be fermionic or bosonic, i.e.
The dispersion relations are assumed such that ǫ A (k) = ǫ B (k) for essentially all k, which is required (see below) for dephasing between the two subsystems to occur as t → +∞. We can assume that the system is prepared in an initial state in contact with a thermal reservoir at a temperature T , i.e. ρ 0 = Z −1 0 e −H0/T (such that Tr ρ 0 = 1). For t > 0, the coupling between the two subsystems H AB disappears, and the two subsystems evolve unitarity and uncoupled, according to the Hamiltonian: The existence of the coupling H AB for all t ≤ 0 implies that in the initial state, ρ 0 , there are correlations (i.e. bi-partite entanglement) between the eigenmodes, i.e. a † k b k = 0. According to the conjecture of Rigol et al. [6], the asymptotic state of the system can be described by a 'generalized' Gibbs ensemble (GGE) density matrix that is obtained using the maximum entropy principle taking into account that the system dynamics is constrained by the existence of the set of integrals of motion given by I a (k) = a † (k)a(k) and I b (k) = b † (k)b(k). The GGE density matrix thus obtained reads: where the Lagrange multipliers are determined by the initial conditions, i.e. α(k) = ln[(1 ± n a (k))/n a (k)] and α(k) = ln[(1 ± n b (k))/n b (k)], with n a (k) and n b (k) given by (9) and (10) (the + applies to bosonic and the − to fermonic modes). Alternatively, one can arrive at an equivalent result by a different route [12]. Let us first consider the expansion of a local operator in terms of the eigenmodes of H: (6) At asymptotically long times after the quantum quench, provided ǫ A (k) = ǫ B (k) and certain conditions of smoothness are met, dephasing renders the two-point correlation function O † (x, t)O(0, t) to the following form [12]: Thus, we see that the asymptotic correlations of O(x) depend only on the eigenmode occupation in the initial state, a behavior that has been termed non-ergodic in Ref. [12]. The above sum over k in Eq. (7) allows to define a mode-dependent temperature for each mode [12]. Indeed, this statement is equivalent to the GGE (cf. Eq. 5) for a broad class of (Gaussian) initial states (see Ref. [12] and below). Thus, it follows that: The above result, valid for local operators, can be combined with Wick's theorem to show that the asymptotic behavior of non-local operators is also described by the GGE [12]. Alternatively, when the correlations are bi-partite, we can regard the effective temperature for the modes in the subsystem A as due to their entanglement with the modes in the subsystem B (and viceversa). Thus, whenever we are dealing with a † k a k = Tr we can trace out one of the subsystems, and write: where ρ A = Tr B ρ 0 and ρ B = Tr A ρ 0 . Therefore, the GGE density matrix can be written as: We can regard the result in Eq. (11) as a way to relate the density matrix of the GGE ensemble to the reduced density matrices of the subsystems A and B. Furthermore, since both ρ A and ρ B are hermitian, it is possible to write these objects as follows [19,20]: where we have introduced the entanglement Hamiltonian of the A (B) subsystem H A(B) , which is also a hermitian operator. Thus, we see ρ GGE is determined by the total entanglement Hamiltonian, The reduced density matrix describing a subsystem A of either a pure state or a thermal mixed state plays an important role in quantum information theory applied to condensed matter systems [19]. For a pure state, the von Neumann entropy S A = − Tr ρ A log 2 ρ A measures the entanglement between two subsystems A and B. The latter can be expressed in terms of the entaglement spectrum of H A . Recently, the von Neuman entropy and entanglement spectrum have become an important tool, as they can be used to characterize topological quantum phases in various kinds of quantum systems, such like graphene [21], topological insulators [22], and quantum spin chains [23]. In this context, an important question that has been addressed in recent times are the conditions under the entanglement Hamiltonian H A can be proportional to the subsystem Hamiltonian H A . Some examples of this fact been discussed in the literature [22,[24][25][26][27]. As we show in the example below (see section II), when this happens to be the case in a system like the one described above, we can expect the asymptotic correlations after the quantum quench to become essentially thermal.
Using the methods of Refs. [28], the entanglement Hamiltonian H A(B) can be determined for a (Gaussian) initial state of the form . Thus [28], which, by comparison with Eq. (5), allows us to identify the Lagrange multipliers α(k) and β(k) of the GGE with the entanglement spectrum of the subsystems A and B. Thus, the entanglement spectrum determines the asymptotic state following a quantum quench. Similar ideas have been discussed by Qi et. al. for the particular case of two-coupled edge states using boundary conformal field theory [22]. Conversely, provided Lagrange multipliers α(k) and β(k) could determined experimentally, we would be able to access the entanglement spectrum and the von Neumann entropy of the subsystems A and B. However, in actual experiments it may be difficult to obtain the full functional dependence of α(k) and β(k). Thus, below we shall focus on two cases where the entanglement spectrum takes a simple form, which may be easier to measure experimentally.
The rest of this work is organized as follows. In the following section, using the above results, we provide an example of the case in which the asymptotic correlations are essentially thermal, which we show to be a consequence of the entanglement hamiltonians H A(B) to be proportional to the subsystem Hamiltonian, H A(B) . In section III, we provide a counter-example of the fact that bi-partite entanglement does not always lead to thermal correlations. This counter-example illustrates the observation that thermal correlations appear provided the bi-partite entanglement arises due from a gap in the spectrum of the Hamiltonian that determines the initial state. Finally, in section IV, we provide a discussion of our results and show that, even when the correlations look essentially thermal, there are certain observables like the energy fluctuations that still differ from their thermal values, a fact that signals a breakdown of the fluctuation-dissipation theorem in the asymptotic state. The appendix contains some technical details regarding a continuum version of the model discussed section II.

II. EXAMPLE
Let us first consider a model that has been numerically studied earlier by Rigol and coworkers [7]. The model describes a system of hard-core bosons in 1D that initially (i.e. for t ≤ 0) move in the presence of superlattice potential. The hard-core bosons in 1D can be treated using a Jordan Wigner transformation [29,30], which maps the hard-core bosons to non-interacting fermions and, in the case of a superlattice of strength ∆, leads to the following quadratic Hamiltonian: where f † j and f j are creation and annihilation operators for spinless fermions at site j (j = 1, . . . , L, for a lattice of L sites). Rigol et al. [7] numerically found that, starting from the ground state of H 0 , if the superlattice term ∝ ∆ is suddenly switched off at t = 0, and the system is allowed to evolve unitarity according to the long-time behavior of the momentum distribution can described by a standard Gibbs canonical ensemble, for which the effective temperature, T eff , was found to approach ∆/2 for ∆ 1. In what follows, we will analytically demonstrate that this numerical observation indeed follows from the existence of a strong bi-partite entanglement between two sets of eigenmodes of H.
We begin by Fourier transforming H 0 and H by using [32] with x j = j, the Hamiltonian (15) can be written as where ω k = −2 cos k and −π/2 < k ≤ π/2. The Hamiltonian describing the state the system at t < 0, namely H 0 , can be brought to diagonal form by means of the following canonical transformation: with tan 2θ k = ∆ ω k . Hence, where E k = ω 2 k + ∆ 2 . Note that the transformation in (21) and (22) implies the existence of strong bi-partite quantum correlations (i.e. entanglement) between the modes at k and k + π, which manifest in, e.g. f † k f k+π = − 1 2 sin 2θ k = − ∆ 2E k = 0. As discussed in Refs. [9,12], the asymptotic momentum distribution of the hardcore bosons for t → +∞ can be obtained from the Fourier transform of the oneparticle correlation function of the bosons, which in turn can be written as a Toeplitz determinant involving correlation two-point correlations of the Jordan-Wigner fermions: where The thermodynamic limit is implicitly understood in the above expressions; we have also introduced the notations: The above correlation function of the (local) operators A i and B j can be shown to be: Next, we use that (K = 0, π): where ρ 0 = |Φ 0 Φ 0 |, Φ 0 being the initial state, that is, the ground state of H 0 (cf. Eq.(15)) and the reduced density matrices: where S 0 = (−π/2, +π/2], and S π = (−π, +π] − S 0 . In other words, the asymptotic correlations can obtained from the reduced density matrix resulting from tracing one one the two sets of entangled modes with k belonging to either S 0 or S π . As explained in section I, the reduced matrices ρ K=0,π can be obtained analytically in terms of the occupation numbers of n(k + K) = f † k+K f k+K in the initial state. Using (21) and (22), we find: and, following the discussion in section I, the Lagrange multipliers determining the GGE density matrix read: For ∆ ≫ ω k , E k can be approximated by ∆, and therefore α(k) ≃ 2ω k /∆ ≃ −β(k). Thus, the GGE density matrix, Eq. (5) reduces to a standard Gibbs ensemble, Eq. (17) with which is in agreement with the numerical observations of Rigol and coworkers [7]. However, it is important to note that the above thermal ensemble and the result of Eq. (33) is only an approximation to the actual GGE ensemble determined by the Lagrange multipliers in Eq. (32). However, this approximation becomes better and better for larger values of ∆, which implies that numerically (and experimentally) the GGE and a standard thermal Gibbs ensemble will be essentially indistinguishable. It is worth noting that the above results are also relevant for a special limit an integrable field theory in one dimension, namely the sine-Gordon model: where a 0 is a short-distance cut-off, and the phase and density fields θ(x) and φ(x), are canonically conjugated in the sense that the obey: ; v is the speed of sound and K is a dimensionless parameter that determines the ground correlations of the system. Indeed, in equilibrium, the model exhibits two phases for g > 0, namely, a gapped phase for K < 2 and gapless phase for K ≥ 2 [29]. For K = 1, as described in the the Appendix, this model can be mapped onto a model of massive fermions in 1D: with ε 0 p = vp and ∆ ∝ g. The model in Eq. (36) can be regarded as the continuum limit of the model of Eq. (19).
For large values of the gap, ∆ ≫ va −1 0 , the multipliers α(p) and β(p) (see calculation in the Appendix) become proportional to the dispersion relation of the Hamiltonian that governs the time evolution, ±ε 0 p respectively, with an effective temperature T eff (it is the same for the two branches of fermions) given by We notice that the effective temperature depends on the temperature of the initial thermal state. This effective "final" temperature is always higher than the initial temperature because it contains the bi-partite entanglement between k modes. For high initial temperatures the effective temperature to which the system thermalizes is the same as the initial one. In the case of a zero-temperature initial state, the effective temperature results T eff = ∆/2 similarly to the case of XY model studied previously.
Generally speaking, the analysis described above is also related to the discussion of the conditions under which the entanglement Hamiltonian H A(B) and the Hamiltonian of the subsystem H A(B) are (approximately) proportional to each other. If this is the case, then, according to the discussion of section I, the GGE density matrix, Eq. 5, will be well approximated by the thermal density matrix of (17). Indeed, recently Peschel and Chung [27] addressed the problem of the proportionality between H A(B) and H A(B) . By considering a model of two two species of fermions with opposite dispersion ω k and coupling ∆, which leads to an energy spectrum for the coupled system with a gap of magnitude 2∆. Using perturbation theory, they showed for The thermal correlations obtained here can be thus regarded as direct consequence of this result when we apply it to a quantum quench where the coupling ∆ is switched off at t = 0 and we exploit the relationship between the GGE and the reduced density matrices ρ A(B) described in section I. Another interesting consequence of the above result is the possibility to use quantum quenches to prepare systems with exactly solvable dynamics in states whose correlations will become indistinguishable from thermal correlations after the quench. However, it is unfortunate that the requirement of a large gap (i.e. the condition that ∆ ≫ 1), implies that thermal states that can be thus obtained are characterized by extremely high effective temperatures (cf. Eq. 33).

III. COUNTER-EXAMPLE
The above result on the emergence of thermal behavior at long times stems from the existence of strong bi-partite entanglement in the initial state. However, as we show in this section, the existence of such entanglement is not a sufficient condition for the emergence of thermal correlations. Indeed, whenever the initial state is gapless, no thermal behavior can be expected even for large entanglement. The Luttinger model [29] is one such example, as we show below.
Let us consider a quantum quench in the Luttinger model (LM) [10]. The initial state is assume to be a mixed thermal state In the above expression, ρ R(L) is the density of the right (left) moving fermions in the Luttinger model, which propagate with Fermi velocity (v F ) [10]. The densities obey the commutation rule [ρ α (k), ρ β (k ′ )] = kL/2πδ k+k ′ δ α,β (α, β = R, L). Therefore we can define two pairs of bosonic operators: ρ L (k) = kL/2πa † k , ρ R (k) = kL/2πb † k and ρ L (−k) = ρ † L (k), ρ R (−k) = ρ † R (k). Thus, the Hamiltonian describing the system at t < 0 can be written as (40) Note that the LM is different from the example that has been discussed in section II because the term that couples the two subsystems is proportional to k, while in the previous example (cf. Eq. 15) it was a constant. This makes the Hamiltonian H LM 0 gapless, as we show in the following paragraph.
A standard way to diagonalize (40) is to introduce a bosonic canonical transformation: Choosing tanh (2φ k ) = −∆/2π. The initial Hamiltonian reads now where Ω k = v F k (1 − (∆/2πv F ) 2 ). Thus, the energy spectrum of H LM 0 is is gapless. According to the discussion of section I and using the methods of Ref. [12], after turning off the interaction described by H int ∝ ∆ at t = 0, the asymptotic behavior of the correlations can be described by a GGE matrix, which can be written as ρ GGE as in Eq. (5) with α(k) and β(k) given by the entanglement spectrum of the subsystems A and B, of modes a k , a † k and b k , b † k respectively. The occupation numbers are: n a (k) = n b (k) = sinh 2 (φ k ) = 1/2(v F k/Ω k − 1). Hence, α(k) and β(k) in Eq. (5) are: where ε = 2[ln (2πv F /∆ + (2πv F /∆) 2 − 1)] is a constant. Hence, the entanglement Hamiltonian take the form: . It then follows that he density matrix of GGE, ρ GGE = Z −1 GGE e −(HA+HB ) does no longer reduce to a thermal ensemble. Thus, we conclude that the existence of bipartite entanglement in the initial state is not a sufficient condition for the emergence of asymptotic thermal correlations. An additional condition, such us the existence of a gap, is appears to be required.

IV. DISCUSSION AND CONCLUSION
In this section, we would like to discuss a number of points concerning the above results. The first point is concerns the opposite situation to the one discussed in section II. We could ask what happens if initially the two subsystems A and B are not coupled and, at t = 0, they suddenly become coupled so that entanglement is created. Can we still expect the asymptotic correlations following such a quench to be described by an essentially thermal ensemble in a certain parameter regime? The answer is no, as we explain below.
To address the above question, let us imagine that, initially, the bosons are free to hop everywhere and there is no superlattice. However, at t = 0 a superlattice is imposed, say by the sudden application of an extra pair of counter-propagating laser beams. Mathematically, the system at t ≤ 0 is described by the Hamiltonian H (cf. Eq. 20) and its subsequent evolution at t > 0 is described by H 0 (cf. Eq. 19). The density matrix of GGE is determined by the occupation numbers n γ (k) = Tr ρ 0 γ † k γ k = sin 2 θ k and n δ (k) = Tr ρ 0 δ † k δ k = cos 2 θ k , for a half-filled lattice. Proceeding as above, the GGE density matrix describing the asymptotic correlations is: for ∆ ≫ 1. Although the above result may appear to be a thermal ensemble, we must recall the dispersion of the eigenmodes is not ω k = −2 cos k but E k = ω 2 k + ∆ 2 . Thus, the temperature becomes again mode dependent and equal to T (k) = ∆E k /ω k , which is consistent with the numerical results of Ref. [6], where lack of thermalization to the standard Gibbs ensemble but thermalization to the GGE was numerically found in this case.
Thus, it appears again that the emergence of thermal behavior in exactly solvable models requires, at least in the simplest case, the existence of an energy gap in the spectrum of the Hamiltonian that determines the initial state. However, as we have already briefly mentioned in section II, the thermal ensemble is just an approximation to the more general GGE, which applies in all circumstances. Indeed, the GGE can reproduce the behavior of the asymptotic correlation functions but it cannot reproduce the behavior of all observables [10]. This is because, in its simplest version of Eq. (5), the GGE does not capture all the correlations between the eigenmodes that exist in the initial state. For example, in the example of section II, I k I k+π GGE = I k GGE × I k+π GGE = Φ 0 |I k I k+π |Φ 0 , where |Φ 0 is the initial state (i.e. the ground state of H 0 , Eq. 19) and I k+K = f † k+K f k+K . This has important consequences, for example, when considering the energy fluctuations: On the other hand, if we compute the same quantity using the GGE, we find: As we have shown in section II, for ∆ ≫ 1, the GGE tends to a thermal Gibbs ensemble (TGE) with T eff = ∆/2 = β −1 eff (ρ GGE → ρ T GE ). And according to the fluctuation-dissipation theorem, for a thermal ensemble ρ T GE at a temperature T eff , where C V the heat capacity of the system. Yet, the actual energy fluctuations of the system following a quantum quench in which a superlattice of strength ∆ ≫ 1 is turned off at t = 0 are given by σ 2 = 2σ GGE ≃ 2σ T GE (cf. Eq. 47). Therefore, we conclude that the fluctuationdissipation theorem breaks down in the model of section II, in spite that the asymptotic correlations appear to be essentially thermal for ∆ ≫ 1. Nevertheless, although the effective temperature obtained in Eq. (33) has no thermodynamic meaning in the sense of the fluctuation-dissipation theorem, it still represents a quantity that is worth determining experimentally. The reason is that T eff is measure of the entanglement in the system. Indeed, determination of T eff from, say, a measurement of the boson momentum distribution, should allow for a determination of the entanglement spectrum of H A(B) , which, according to the discussion in sections I and II, is given by H A(B) ≃ H A(B) /T eff . Thus, an experimental determination of the von Neumann entropy, S A(B) = − Tr ρ A(B) log 2 ρ A(B) , (ρ A(B) = e −HAB /Z A(B) ) would be also possible. Similar remarks are applicable to the counter-example discussed in section III, provided we exchange the role of T eff by ǫ. In this case, however, we cannot expect thermal correlations.
In summary, we have presented a simple instance of a quantum quench in which a quantum quench in an exactly solvable system can produce essentially thermal correlations. The emergence of thermal correlations from the generalized Gibbs ensemble (GGE) has been related to the existence of bi-partite eigenmode entanglement and a gap in the spectrum of the Hamiltonian that describes the initial state. In this regard, we have also discussed a counter-example demonstrating that thermalization does not happen if the initial state is described by a gapless Hamiltonian. Our results allow to establish a link between the GGE and the entanglement spectrum in exactly solvable systems with bi-partite entanglement of the eigenmodes. Thus, it makes it possible an experimental measurement of the entanglement spectrum and other quantities derived from it (such like the von Neumann entropy), provided the asymptotic correlation functions of the system following a quantum quench can be measured experimentally. We have argued that this task becomes particularly simple when the GGE reduces to a thermal ensemble (or, when the entanglement spectrum has a relatively simple form, as in the counter-example discussed in section III). Finally, we have also shown that, even if correlations may become essentially thermal, other quantities, such as the energy fluctuations, are not. This is akin to a breakdown of the fluctuation-dissipation theorem.