{'fraction_non_alphanumeric': 0.04262358894511483, 'fraction_numerical': 0.016105488516932658, 'mean_word_length': 4.083436392118064, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reduced Density Matrix after a Quantum Quench

We consider the reduced density matrix (RDM) \rho_A(t) for a finite subsystem A after a global quantum quench in the infinite transverse-field Ising chain. It has been recently shown that the infinite time limit of \rho_A(t) is described by the RDM \rho_{GGE,A} of a generalized Gibbs ensemble. Here we present some details on how to construct this ensemble in terms of local integrals of motion, and show its equivalence to the expression in terms of mode occupation numbers widely used in the literature. We then address the question, how \rho_A(t) approaches \rho_{GGE,A} as a function of time. To that end we introduce a distance on the space of density matrices and show that it approaches zero as a universal power-law t^{-3/2} in time. As the RDM completely determines all local observables within A, this provides information on the relaxation of correlation functions of local operators. We then address the issue, of how well a truncated generalized Gibbs ensemble with a finite number of local higher conservation laws describes a given subsystem at late times. We find that taking into account only local conservation laws with a range at most comparable to the subsystem size provides a good description. However, excluding even a single one of the most local conservation laws in general completely spoils this agreement.


I. INTRODUCTION
Recent advances in systems of optically trapped ultracold atomic gases have made it possible to observe the nonequilibrium time evolution of isolated many particle systems over long time scales [1][2][3][4][5][6] . A key property of such cold atomic clouds is their weak coupling to the environment and resulting smallness of external dissipative processes. To a good approximation one is dealing with isolated quantum mechanical many particle systems, which are prepared in generally mixed states, and one is interested in the time dependence of observables, in particular at late times. The experimental results have stimulated theoretical efforts aimed at understanding the principles underlying the nonequilibrium dynamics of isolated many particle systems. Some of the most basic questions are whether observables generally relax to time-independent values, and if they do, whether their stationary values are described by a statistical ensemble. Other prominent issues concern the roles of dimensionality and conservations laws. Experiments on trapped 87 Rb atoms 2 established that three-dimensional condensates rapidly relax to a stationary state characterized by an effective temperature, whereas constraining the motion of atoms to one dimension leads to a much slower relaxation to a nonthermal distribution. It was argued that this observed difference has its origin in the presence of additional (approximate) conservation laws, related to quantum integrability, in the one dimensional case. Theoretical efforts aimed at understanding these and related questions  indicate that in translationally invariant models there are at least two basic types of behaviours at late times: subsystems either thermalize 46 , i.e. are characterized by a Gibbs distribution with an effective temperature, or they are described by a generalized Gibbs ensemble (GGE) 8 . There is evidence that the latter applies to integrable models, while the former constitutes the "generic" situation.
A popular protocol for analyzing nonequilibrium evolution is a so-called quantum quench: here the system is originally prepared in the ground state |Ψ 0 of some local, short ranged Hamiltonian H(h 0 ), where h 0 is a system parameter such as a magnetic field or an interaction strength. At time t = 0, h 0 is then suddenly "quenched" to h, and the subsequent time evolution under the new Hamiltonian H(h) is studied. Under this protocol the system remains in a pure state |Ψ t = exp(−iH(h)t)|Ψ 0 at all times, and as a whole can clearly never be described by a Gibbs or generalized Gibbs distribution. This can be seen by considering the hermitian operators O (n,m) = |n m| + |m n|, (1.1) where |n and |m are eigenstates of H(h) with energies E n and E m respectively. Then the expectation values Ψ t |O (n,m) |Ψ t = Ψ 0 |n m|Ψ 0 e i(En−Em)t + h.c. (1.2) are oscillating in time and never become stationary. A useful and intuitive point of view is to focus on local properties of a given system in the thermodynamic limit, i.e. ask questions only about observables contained in a finite subsystem A 13,31,32 . Here the (infinitely large) comple-mentĀ of A can act as an effective bath, and probability may freely dissipate from A toĀ. As a result A may be described by a mixed state. Arguably the most precise and convenient description of this situation is in terms of the reduced density matrix ρ A (t) of subsystem A. The latter is obtained from the density matrix ρ(t) = |Ψ t Ψ t | of the entire system as If (1.4) holds, then the system evolves towards a stationary state described by the distribution ρ stat . In particular (1.4) implies that the expectation values of any local operator O A acting only within subsystem A is given by (1.6) An efficient way of investigating whether a given RDM approaches a known stationary distribution at late times was introduced in Ref. [20] by considering the operator norm ρ A (t) − ρ stat,A op . If this approaches zero at late times, then the system relaxes locally to the stationary distribution ρ stat . Ref. [20] was concerned with the case where ρ stat describes a thermal ensemble with a given effective temperature, and considered very small subsystems. Here we are interested in a quench to a quantum integrable model. As we have alluded to before, the stationary state for such quenches is believed to be described locally by a generalized Gibbs ensemble (for the model we consider below this was established in Ref. [32]). More precisely, the density matrix of the entire system is expected to be of the form where Z is a normalization 1 , and I n are local conserved quantities, i.e. local operators such that We stress that the GGEs considered in the quench context are fundamentally different from thermal ensembles, because through the specific values of the Lagrange multipliers they retain an infinite amount of information about the initial state. Above we have stipulated that only local (in space) conservation laws I n are to be included in the definition of ρ GGE , but it is in fact a matter on ongoing debate whether locality is a necessary or even desirable requirement 47 . In this context a result obtained in Ref. [31] is rather illuminating: there it was demonstrated for a particular example, the transverse field Ising chain, that different statistical ensembles can have identical local properties. The two ensembles considered were a GGE of the form (1.7), and the so-called "pair ensemble" obtained by time averaging the quench density matrix ρ(t). Given that ρ stat is generally not unique, it is clearly desirable to identify the simplest description. To that end we introduce truncated generalized Gibbs ensembles of the form ρ (n0) tGGE = 1 Z e n<n 0 λnIn , (1.10) and investigate how well such ensembles describe the stationary state for quenches to integrable models. The outline of this paper is as follows. In section II we review some relevant results for the transverse field Ising chain. Local conservation laws are presented in section III and used in sections IV, V, VI to define several classes of generalized Gibbs ensembles. Properties of corresponding reduced density matrices are discussed in section VII. In section VIII we discuss general properties of distances on the space of reduced density matrices and introduce the distance used in the remainder of the paper. In sections IX, X, XI and XII we present results for the distance between quench and generalized Gibbs reduced density matrices. We summarize our results in section XIII. Various technical issues are discussed in four appendices.

II. SOME FACTS ABOUT THE TRANSVERSE-FIELD ISING CHAIN (TFIC)
Here we briefly review some relevant results on the TFIC. The latter is an important paradigm for quantum phase transitions in equilibrium 48 as well as nonequilibrium dynamics 14,24,29,31,32,36,49 . In the latter context experimental realizations range from cold atomic gases 50 to circuit QED 51 . The Hamiltonian of the model on a ring is where σ α L+1 = σ α 1 . The quantum spins can be mapped to (real) Majorana fermions by means of the Jordan-Wigner transformation where {a i , a j } = 2δ ij . In terms of the Majorana fermions (2.2), the Hamiltonian takes a block-diagonal form Here N is the number operator and by construction e iπN = j σ z j commutes with H R,NS . The two blocks H R and H NS correspond to periodic and antiperiodic boundary conditions on the fermions respectively. They can be diagonalized by Bogoliubov transformations where the Bogoliubov angle θ p is given by The diagonal form of the Hamiltonian is where the single-particle energy is given by The ground states of H R,NS (h) are the fermionic vacua Here the momenta are p = πn L , where n are even/odd integers for R and NS fermions respectively.

A. Quantum Quenches
Our quench protocol is as follows: we prepare the system in the ground state |Ψ 0 for an initial value h 0 of the transverse magnetic field. At time t = 0 we instantaneously change the field from h 0 to h. The state of the system at times t > 0 is obtained by evolving with respect to the new Hamiltonian H(h) (2.10) An important quantity is the difference ∆ k = θ k − θ 0 k of Bogolibov angles required to diagonalize H(h) and H(h 0 ) respectively As we are interested in obtaining results in the thermodynamic limit we have to distinguish between two cases. Here the initial state in a large, finite volume is simply the NS vacuum The time evolved state can then be written in the form 31 Here our initial state in a large, finite volume must reflect the spontaneous symmetry breaking of the Z 2 spinflip symmetry σ x,y → −σ x,y in the thermodynamic limit. The appropriate choice is 31 (2.14)

III. LOCAL CONSERVATION LAWS IN THE TFIC
We consider the one dimensional transverse field Ising chain in the thermodynamic limit Following Ref. [52] we can construct an infinite number of local conservation laws I ± n [I α n , I β m ] = 0 , n = 0, 1, . . . , α, β = ±, (3.2) where the Hamiltonian itself is H = I + 0 . Let us define operators and In terms of these operators the local conservation laws are They are local in the sense that the density of I α n involves only spins on n + 2 neighbouring sites. By virtue of their locality, the conservation laws can all be expressed in terms of Jordan-Wigner fermions (2.2) We now realize that all conservation laws (3.6) are in fact quadratic in Majorana fermions! It is then a simple matter to diagonalize them simultaneously by means of a Bogoliubov transformation (2.5), where on the infinite chain the Bogoliubov fermion operators have anticommutation relations The conservation laws take the simple form which furthermore shows that they are even/odd under spatial reflections. Interestingly the conservation laws I − n do not depend on the transverse field h and are therefore shared by the entire one-parameter family of Hamiltonians H(h). This seems to be a generic feature of models with a free fermion spectrum like the TFIC, see App. C and Ref. [58].

A. Local conservation laws for periodic boundary conditions
Above we focussed on the bulk contribution to the local conservation laws. For a finite system on a ring there are boundary contributions, which can be determined as follows. In terms of the Bogoliubov fermions, the conservation laws for periodic boundary conditions are where the momenta are taken to be either in the R or NS sectors. By inverting the Bogoliubov transformation and Fourier transforming back to position space, one obtains a representation of the conservation laws in terms of the Majorana fermions a j for periodic/antiperiodic boundary conditions respectively. Finally, inverting the Jordan-Wigner transformation gives the desired expression in terms of spins.

IV. GENERALIZED GIBBS ENSEMBLE
We now define the density matrix of a generalized Gibbs ensemble formally by the expression where Z is a normalization that ensures Trρ GGE = 1. In practice we need to regularize (4.1) in an asymptotically large, finite volume L. The Lagrange multipliers λ σ n are fixed through the requirements Using translational invariance we can alternatively work with the densities of the conservation laws to rewrite (4.2) as The solution to this system of equations is where l ≥ 0 and cos ∆ k is defined in (2.11). In Fig. 1 we show λ + l for a quench from h 0 = 0.1 to h = 0.7. The large l behaviour of Eq. (4.5) is determined by the regions k ∼ 0, π (where the integrand has a logarithmic singularity) and one can show that where the sign is + for quenches within the same phase and − otherwise. We see that the Lagrange multipliers λ + n decay rather slowly as a function of n.

A. GGE in terms of mode occupation numbers
In the literature the generalized Gibbs ensemble is often constructed from mode occupation numbers n k = α † k α k , see e.g. 8,12,18,53,54 . The latter are non-local (in space) as they involve a Fourier transform. We will now establish the relation between this and our definition (4.1). It follows from (3.8) that the density matrix can be rewritten in the form where This establishes the fact that the GGE density matrix can be constructed either from the local conservation laws (3.6), or from the mode occupation numbers n k . This relationship generalizes to interacting integrable models, where the appropriate GGE can be formulated either in terms of the local integrals of motion generated from the transfer matrix, or from the mode occupation numbers n a (k) = Z † a (k)Z a (k), where Z a (k) are Faddeev-Zamolodchikov operators 55,56 .

V. TRUNCATED GENERALIZED GIBBS ENSEMBLES
In order to assess the importance of the various conserved quantities, it is useful to define ensembles that interpolate between the Gibbs distribution and the GGE. We define particular such truncated GGEs as follows. Given that the densities of the conservation laws I ± n involve n + 2 neighbouring sites, it is natural to retain only the "most local" conservation laws, i.e. Here y is an integer and y = 1 (y = ∞) corresponds to the Gibbs ensemble (GGE). The Lagrange multipliers λ σ n,y are obtained from the requirements Tr (I σ n ) j,...,j+n+1 ρ (y) tGGE = Ψ 0 |(I σ n ) j,...,j+n+1 |Ψ 0 , (5.2) where 0 ≤ n < y. Eqns (5.2) are a consequence of [I σ n , H] = 0 and the assumption that the stationary state after the quench is described by RDMs based on (5.1). For transverse field quenches we have λ − n,y = 0, but the other Lagrange multipliers are different from their respective values in the full GGE, i.e. λ + n,y = λ + n .
We note that the correlation matrix of ρ (y) can be computed efficiently using FFT algorithms. This is in contrast to the case of theories with non-trivial scattering matrices, for which it is extremely difficult to reconstruct the Lagrange multipliers from the conservation laws 57 .

VI. DEFECTIVE GENERALIZED GIBBS ENSEMBLES
It is instructive to consider a second type of truncated GGE, where we retain an infinite, but incomplete set of integrals of motion. Such "defective" GGEs will allow us to ascertain the role of a particular local conservation law. We define the truncated defective GGE as the density matrix (q < y) in which the Lagrange multipliers λ + n,(+q) are fixed by the system (4.2) with n ≤ y, n = q, and we have used that the Lagrange multipliers λ − n,(+q),y must vanish as a consequence of reflection symmetry around the origin. We then define the defective GGE as the limit y → ∞ of truncated defective GGEs: In order to solve the system of equations (4.2) for the defective GGE it is useful to work in the mode occupation number representation (4.7), which reads where the Lagrange multipliers γ Guided by the fact that cos(nk) and sin((n + 1)k) form an orthonormal set of functions on [−π, π], we look for a solution of the form tanh γ (+q) k 2 = cos ∆ k − κ + q cos(qk) ε(k) , (6.5) where κ + q is a yet to be determined constant. We note that the value of κ + q affects the expectation values of local operators. For some cases κ + q can be easily determined as follows. Given that | tanh(x)| ≤ 1, (6.5) implies that Setting k = 0, π then gives Eqs (6.7) allow us to identify cases, in which κ + q = 0: 1. odd q and quenches within the same phase; 2. even q and quenches across the critical point.
the defective GGE is identical to the full GGE. This "GGE reconstruction" is a peculiarity of free-fermion models and can be traced back to the existence of conservation laws independent of the quench parameter, see also Ref. [58].
For general quenches and values of q, κ + q is determined by Eq. (6.2). We find that it takes the value corresponding to the maximal entanglement entropy (as shown in Fig. 9), although the entanglement entropy may be nonstationary under a variation of the excluded integral of motion (see Appendix D for further details).

VII. REDUCED DENSITY MATRICES IN THE TRANSVERSE FIELD ISING CHAIN
In this section we summarize some basic features of RDMs in the TFIC. We note that most of the following discussion generalizes straightforwardly to other spin chains with free fermionic spectra such as the quantum XY model. Our starting point is a density matrix ρ describing the entire system, which we take to be of size L with periodic boundary conditions. We are interested in the limit L → ∞, but it is convenient to start with a large, finite chain. The RDM of a subsystem consisting of ℓ spins 1/2 at sites x i , i = 1, . . . , ℓ can be expressed in the form Tr ρ σ α1 x1 · · · σ α ℓ x ℓ σ α1 x1 · · · σ α ℓ x ℓ , (7.1) where α i = 0, x, y, z and σ 0 ≡ I. The quantum spins are mapped to (real) Majorana fermions by the In the cases of interest to us, the correlation matrix is of block-Toeplitz form The TFIC Hamiltonian exhibits a Z 2 symmetry If the density ρ is invariant under the transformation (7.6) we have Tr [ρ a j ] = 0, and as Wick's theorem applies to the Jordan-Wigner fermions we can express (7.2) as a Gaussian 62 where Z ensures that Trρ ℓ = 1 and W is a skewsymmetric, 2ℓ-by-2ℓ Hermitian matrix related to Γ by We now turn to the three particular cases of interest, namely those where ρ in (7.1) is a thermal density matrix, a GGE density matrix, or the density matrix after a global quantum quench of the transverse field in the TFIC.

A. Thermal Density Matrix
On a very large ring, the Hamiltonian has a block diagonal structure, see section II. The thermal density matrix is a function of the Hamiltonian and therefore inherits the same block structure It follows from this that only even operators have nonvanishing expectation values, i.e.
In the thermodynamic limit the difference between expectation values of local operators with respect to the R and NS sectors tends to zero, so that we may work exclusively in e.g. the R sector. The resulting RDM of a contiguous block of spins is then Gaussian (7.7), (7.8) with It can be written in the form (7.4) with

B. GGE Density Matrix
It was shown in Ref. [32] (see also [24,31]) that the RDM of the generalized Gibbs ensemble (4.1),(4.7) is Gaussian and can be expressed in the form (7.7), (7.8). The correlation matrix is given by It can be written in the form (7.4) with Here the γ k 's are related to the λ σ m 's by (4.8) and the Bogoliubov angle θ k is given in (2.6).

C. Truncated GGE Density Matrix
The correlation matrix of the truncated GGE defined in Sec. V is given by 15) It can be written in the form (7.4) with f (k) = 0 , g(k) = −ie iθ k tanh P y−1 (cos(y))ε(k) .
(7.16) Here P y−1 (x) is a polynomial of order y − 1, which is computed numerically.

D. Defective GGE Density Matrix
In Sec. VI we defined the defective GGE ρ (q) dGGE as the ensemble that lacks in the conservation law I + q . Its correlation matrix is given by (7.17) It can be written in the form (7.4) with f (k) = 0 and (cf. Eq. (7.14)) where κ + q is computed numerically maximizing the entanglement entropy, which selects λ + q,(+q) = 0 whenever it is allowed. We note that the Fourier transform of Eq. (7.18), which is required to compute the correlation matrix (7.5), can be easily expressed in terms of the GGE correlators; for |ℓ| < q we have Since κ + q is a bounded function of q (cf. Eq. (6.7)), at fixed ℓ the fermionic correlators approach the GGE ones at least exponentially fast in q.

E. Quench Density Matrix
At zero temperature the ground state phase diagram of the TFIC exhibits ferromagnetic (h < 1) and paramagnetic (h > 1) phases, separated by a quantum critical point. In the ferromagnetic phase the Z 2 symmetry of the Hamiltonian is broken spontaneously. As we will see, this symmetry breaking has important effects on the time evolution of the density matrix.

Quenches originating in the paramagnetic phase
Here at t > 0 the full quench density matrix is where the state |Ψ t is given in (2.13). As a result of the squeezed-state form of |Ψ t , Wick's theorem applies to averages calculated with respect to ρ(t), and RDMs are Gaussians of the form (7.7), (7.8), with correlation matrix This is of the form (7.4) with where e iθ k is given by (2.6).

Quenches originating in the ferromagnetic phase
Given the initial state (2.14), the post-quench density matrix of the full system is Importantly, RDMs are no longer Gaussian in this case. We will discuss how to cope with this complication in section XI. It is known 32 that in the stationary state RDMs are Gaussian with a correlation matrix equal to the t → ∞ limit of (7.21).

VIII. DISTANCES ON THE SPACE OF RDMS
In the following we focus on RDMs for finite subsystems of lattice models with a finite dimensional Hilbert space at each site. In this case the RDMs are finite dimensional matrices and a simple way to define a distance between two density matrices is by means of a matrix norm Here the index a labels different matrix norms. As we are dealing with finite matrices, all norms are equivalent in the sense that where c ab and c ba are positive numbers that depend on the matrix dimension but are independent of ρ. One consequence of (8.2) is that if the distance between two matrices approaches zero when some external parameter p is tuned to a valuep, the dependence of the distance on p −p is almost independent of the norm. On the other hand, the dependence on matrix dimension is in general very different for different norms. This is important for our purposes, because the matrix dimension is related to the size of the subsystem under consideration, and it is principally desirable to be able to compare distances between different sizes. From a technical point of view, the distance induced by the Frobenius norm 2 is generally the easiest to calculate. On the other hand, it has the drawback that the physical interpretation of the distance is less transparent than for some other norms. For instance, given two density matrices ρ and ρ ′ , a very natural question is how different expectation values of local observables are in the two ensembles. We now discuss this question for the particular case of spin-1/2 quantum spin chains. Here the most important local observables are products of Pauli matrices. These are particular cases of involutionsÔ 2 = I, for which the following inequality holds Here is the trace norm. From Eq. (8.4) it is evident that the trace distance provides an upper bound for the difference between the expectation values of observables in the two states: if ρ − ρ ′ 1 < ǫ, then the expectation values of all (local) observables will agree in the two ensembles within accuracy ǫ. In terms of the Frobenius distance we have instead (here we use that the local Hilbert space is two-dimensional) On the other hand, we have where in the last step we have used that for involutionŝ O Combining (8.7) and (8.6) we see that as long as ρ − ρ ′ F 2 1−ℓ/2 , the Frobenius distance does not provide useful information about expectation values. It is shown in Appendix A that for sufficiently large ℓ this is always the case. A second problem with using the Frobenius norm as a distance is that the norms of RDMs at late times after a quantum quench, as well as the norms of RDMs describing Gibbs or generalized Gibbs ensembles, generally are exponentially small in the subsystem size. Given the upper bound derived in Appendix A 9) this implies that in these cases of interest ρ − ρ ′ F is exponentially small in subsystem size. This shows that the Frobenius norm itself is not a convenient measure for the distance between two density matrices. The same problem occurs for the operator norm ||A|| op = √ λ max , where λ max is the largest eigenvalue of A † A. This norm was used for example in Ref. 20 to analyze the relaxation properties of small subsystems after a quench into a nonintegrable model. Indeed we have and the maximal eigenvalues of RDMs for large subsystems are generally exponentially small in subsystem size.

A. Definition of the Distance
In order to circumvent the problem described above, we define our "distance" 3 on the space of RDMs as 1. An upper bound: Using the upper bound derived in (A1) we see that 2. A lower bound: A lower bound for D(ρ, ρ ′ ) can be established by means of the triangle inequality ||ρ− ρ ′ || F ≥ | ||ρ|| F − ||ρ ′ || F |. Using that the Frobenius norm of a RDM is related to the second Rényi entropy by we find that This provides the desired lower bound We note that the bound (8.15) is independent of subsystem size ℓ as long as the second Rényi entropies of the two ensembles differ at least by a constant (when viewed as functions of ℓ).

B. The Distance between two Thermal Ensembles
In order to establish a benchmark for (8.11), it is useful to consider the distance between the RDMs of two thermal ensembles at slightly different inverse temperatures β and β ′ (but the same Hamiltonian). Then For a sufficiently large subsystem (and a local Hamiltonian), the first factor on the right hand side can be expressed as For a large subsystem, this is proportional to the square of its size, and hence We conclude that the distance between two thermal RDMs on a subsystem of size ℓ, and β ≈ β ′ is As expected this is proportional to the difference in temperatures, but there is also a factor of ℓ. The latter is important if one is interested in comparing the distance between two ensembles for different subsystem sizes.

C. The Distance between two GGEs
The above discussion carries over to the case of two generailzed Gibbs ensembles (4.1), with slightly different values of Lagrange multipliers λ σ m . The leading contribution to the distance is given by A calculation similar to the thermal case shows that for large subsystem size

D. Information on observables contained in the distance
Let us consider the situation where the distance between two reduced density matrices ρ 1 and ρ 2 defined on an interval of length ℓ becomes small, and denote the corresponding averages of local operators on said interval by By expanding the density matrices in a complete basis of Hermitian involutions we can show that .
Defining an average we can express the distance (8.11) as is the relative difference between the ensembles described by ρ 1 and ρ 2 , respectively. This implies that D(ρ 1 , ρ 2 ) measures the mean relative difference of the expectation values of all local operators, averaged with respect to the probability distribution (8.24).

E. The Distance between two Gaussian Density Matrices
The distance (8.11) between two Gaussian RDMs ρ[Γ] and ρ[Γ ′ ] can be expressed in terms of their correlation matrices following Ref. [59]. Given the definition of the distance we require tractable expressions for the quantities This is achieved in two steps. First, we note that the product of two Gaussian RDMs (7.7) is itself Gaussian 59 This can be seen by expanding the left hand side of (8.29) by means of the Baker-Campbell-Hausdorff formula in terms of multiple commutators, and then observing that the commutator of quadratic operators is quadratic and gives rise to a commutator between the matrices W that characterize the 2-forms Given that the correlation matrices are only 2ℓ dimensional (with ℓ the subsystem size), (8.33) provides a very efficient way of computing distances for large subsystem sizes.

IX. SINGLE-SITE SUBSYSTEM
It is instructive to consider the time evolution of the RDM describing a single-site subsystem in some detail. In this case the RDM of site 1 can be expressed in the form where m(t) is the magnetization per site at time t after the quench, i.e. 2) The RDM of the generalized Gibbs ensemble describing the stationary state is Finally, the RDM of the thermal ensemble described by ρ β , whose inverse temperature β is fixed by the requirement is given by Here the transverse magnetization per site is A. Quenches originating in the paramagnetic phase Here the Z 2 symmetry enforces m x (t) = m y (t) = 0. For late times we may evaluate the integral by means of a stationary phase approximation, which gives where The distance between ρ 1 (t) and the generalized Gibbs RDM at late times then decays to zero like a power-law with exponent 3/2 On the other hand, the distance between ρ 1 (t) and the thermal RDM approaches a constant at late times . Here all three components of the magnetization per site are non-zero. The component along the transverse field direction is again given by (9.9), while the late-time asymptotics of m x (t) has been calculated in 31 (9.14) Here C x F F is a known amplitude and ε ′ k = dε h (k) dk . Finally, the Heisenberg equation of motion for σ x 1 (t) relates the y and x components m y (t) = 1 2Jh Importantly, m x,y (t) exhibit exponential decay in time.
In contrast, m z (t) again decays like a power law with exponent 3/2 and therefore will dominate the late time behaviour. Hence at sufficiently late times, the distances of ρ 1 (t) to GGE and thermal RDMs are again given by (9.13) and (9.12) respectively. So for a single site subsystem the spontaneous symmetry breaking only modifies the intermediate time behaviour of the distances. As we will see, this holds true also for larger subsystems.

X. LARGER SUBSYSTEMS FOR QUENCHES FROM THE PARAMAGNETIC PHASE
For quenches with h 0 > 1 and in the thermodynamic limit, we determine the distance between the quench RDM and that of an appropriate thermal or generalized Gibbs ensemble by means of relation (8.33). The correlation matrices for all cases are of the form (7.4), (7.5) with elements given in (7.12), (7.14) and (7.22) respectively. For a subsystem of size ℓ this requires the calculation of determinants of 2ℓ × 2ℓ matrices, which is done numerically. Results for a quench from h 0 = 1.2 to h = 3 and subsystem sizes ℓ = 10, 20, 30 . . . , 150 are shown in Figs 2 and 3 We see that the distance between quench and Gibbs RDMs tends to a ℓ-dependent constant at late times. This establishes that subsystems do not thermalize. On the other hand, as can be seen from Fig. 3, at sufficiently late times the distance between ρ ℓ (t) and ρ GGE,ℓ decays to zero in a universal power-law fashion large-ℓ asymptotics of the function k(ℓ) can be inferred as follows. On surfaces with constant, small D the time scales as t ∼ ℓ 4/3 as is shown in Fig. 5. This in turn implies that

A. Relaxation time
We may extract a relaxation time from the behaviour of the distance, by using the connection to averaged differences in the expectation values of local operators established in subsection VIII D. The distance can be written as where 4) and the bar denotes the average (8.24). Using that , (10.5) and then substituting the asymptotic behaviour (10.1), (10.2) into the right hand side, we obtain Bounding the right hand side by a (small) constant, we obtain a time scale t * rms associated with the relaxation of the average relative error with respect to the distribution (8.24) It is not simple to identify the observables that give significant contribution to the average, since it depends both on their "multiplicity" in the subsystem (produced by translational invariance and other symmetries) and on the expectation values. We note that the relaxation time t * rms is very different from the time scales identified in Ref. [32] in the time evolution of the two point functions of spin operators for quenches within the paramagnetic phase.

B. Distance from Truncated Generalized Gibbs Ensembles
Having established that the distance between quench and GGE reduced density matrices tends to zero as a universal power law at late times, a natural question is, how close the quench RDM is to the truncated GGEs (5.1), which retain only finite numbers of conservation laws. A representative example for a quench within the paramagnetic phase is shown in Fig. 6. We see that at sufficiently late times, the distances converge to constant values. However, increasing the range (and number) of conservation laws, the values of these plateaux decrease, signalling that retaining more conservation laws gives better descriptions. In an intermediate time window, the extent of which grows with y, the distance decays in a universal t −3/2 power-law fashion. In Fig. 7 we consider the distance 8) between the RDMs of the truncated and full generalized Gibbs ensembles as a function of the parameter y. For a given subsystem size ℓ, this corresponds to plotting the values of the plateaux seen in Fig. 6 against the corresponding values of y. The distance is seen to start decaying exponentially as a function of y as soon as y ℓ.
There are two main conclusions of the above analysis: 1. Including more local conservation laws improves the description of the stationary state.
tGGE,ℓ ) between the GGE and the truncated GGEs obtained by imposing local conservation laws with densities involving at most y + 1 consecutive sites. The quench is from h0 = 1.2 to h = 3 and the subsystem size ranges from ℓ = 5 to ℓ = 50. Colors and sizes change gradually as a function of the size ℓ. For y > ℓ, the distance starts decaying exponentially in y, with an ℓ-independent decay constant.
2. The description of the stationary state improves rapidly, once the range y + 1 of all conservation laws not included in the truncated GGE exceeds the subsystem size ℓ.

C. Distance from defective Generalized Gibbs Ensembles
We now turn to the role played by particular local conservation laws. We find that the distance between quench and defective GGE reduced density matrices for a given quench and subsystem size tends to a constant at late times, i.e. (10.9) The dependence of this asymptotic value on the subsystem size ℓ and the integer q is shown in Fig. 8 for a quench within the paramagnetic phase. We see that D hibits an exponential decay in q as soon as q ℓ. This is similar to the behaviour observed in the truncated GGE case. The decay length can be calculated from the large-q asymptotics of Eq. (7.19). By series expanding Eq. (8.33) to second order in Γ Numerically we find that κ + q ∼ 1/q 2 .

"GGE Reconstruction"
In section VI we discussed the issue that, for certain quenches and omitted conservation laws I + q , the corresponding defective GGE is identical to the full generalized Gibbs ensemble. We now return to this point. In Fig. 9 we consider the truncated, defective GGE for a quench across the critical point from h 0 = 2 to h = 0.5 for a subsystem of length ℓ = 5. We plot the distance between the reduced density matrices of the appropriate GGE and the truncated, defective GGE with y integrals of motion, where I + q (q < y) has been excluded, i.e. As discussed in section VI, for even q we expect this distance to approach zero, when the number y of conservation laws goes to infinity. This behaviour is clearly observed in Fig. 9. This implies that the corresponding conservation laws do not affect averages of local operators. As discussed before, this is a particular feature of free theories, where H(h 0 ) and H(h) generically share certain local conservation laws. (+q),y dtGGE,ℓ ) for a quench across the critical point between the GGE and the defective truncated GGE RDMs for a subsystem of 5 consecutive spins, as a function of the number y of retained conservation laws. Each symbol corresponds to a different excluded conservation law I + q (the legend indicates the value of q). The distance approaches zero for even q, whereas it remains finite for odd q, in agreement with the discussion of section VI. The lines are the distances from the corresponding defective generalized Gibbs ensemble ρ (+q) dGGE with maximal entanglement entropy.

More local conservation laws are more important
On the other hand, for odd q we find that D d(+q),y ∞ approaches constant values when y becomes large. This value agrees with the distance between the GGE and the defective GGE with maximal entanglement entropy (we stress that for the considered quench the defective GGE does not always correspond to a stationary point of the entanglement entropy under a variation of the excluded integral of motion, as shown in Fig. 20 of Appendix D).
The fact that D d(+q),y ∞ tends to a constant at large y shows that retaining an infinite number of local conservation laws while excluding one of them is insufficient for describing the stationary state. By comparing distances for different values of q we observe that for a given value of y, D d(+q),y ∞ decreases as a function of q. This implies the more local the conservation law, the more important it is for describing the stationary state.

XI. QUENCHES FROM THE FERROMAGNETIC PHASE: EFFECTS OF SPONTANEOUS SYMMETRY BREAKING.
We now turn to quenches originating in the ferromagnetic phase, i.e. h 0 < 1. In this case, the time evolved initial state is given by where |0; h 0 R/NS are the ground states of the Hamiltonian H(h 0 ) with periodic/antiperiodic boundary conditions. In order to analyze reduced density matrices after a quantum quench from the ferromagnetic phase we will make use of the following facts.
a) The fermion parity e iπN = j σ z j Eq. (2.4) is fully factorized in space.
b) The states |ψ t R and |ψ t NS are eigenstates of e iπN with eigenvalues 1 and −1 respectively.
c) The difference between the expectation values of local operators in the states |ψ t R and |ψ t NS tends to zero in the thermodynamic limit.
d) The RDMs TrĀ [|ψ t aa ψ t |], a = R, NS, where A is a single interval andĀ its complement, are Gaussian.
Property b) allows us to express the full density matrix in the form (cf. Eq. (7.1)) We will refer to O e/o as even and odd operators respectively. The main difference between even and odd operators, is that the latter are not local in terms of fermions: a Jordan-Wigner string is attached to them. We are interested in the RDM of a block A of ℓ contiguous spins, which is obtained by tracing out the degrees of freedom outside A ρ ℓ = TrĀ ρ . (11.4) A convenient representation for ρ ℓ is obtained by restricting the sums in Eq. (11.2) to involutions that act as the identity operator outside the interval A, i.e.

O → O (A) ⊗ I (Ā)
, (11.5) where the superscript (A) indicates that the operators act on the Hilbert space over all sites in A. As a result of property a), fermion parity has a simple restriction onto the interval A e iπNA ≡ l∈A σ z l , (11.6) and can be used to subdivide operators O (A) into even and odd ones This then implies that we can decompose the RDMs of (11.2) into even and odd parts as well In the thermodynamic limit we then may employ property c) to obtain the following expressions Importantly, the even part ρ ℓ,e (t) is Gaussian (7.7) by virtue of property d), and has the same structure as RDMs for quenches originating in the paramagnetic phase. On the other hand, the odd part ρ ℓ,o has its origin in the spontaneous breaking of the Z 2 symmetry. The commutation relations (11.7) imply that Tr ρ ℓ,o (t)ρ Ga ℓ = 0 for any Gaussian density matrix ρ Ga ℓ , because the latter is by construction even. As a result the odd part ρ ℓ,o of the RDM enters the distance from a Gaussian state only through its norm (11.10) We will be interested in the cases where ρ Ga ℓ describe Gibbs or (truncated) generalized Gibbs ensembles. The Frobenius norms ||ρ ℓ,e (t) − ρ Ga ℓ || F , ||ρ ℓ,e || F and ||ρ Ga ℓ || F can be efficiently evaluated by means of Eq. (8.33). What remains in order to determine the distance (11.10) is a method for calcuating the Frobenius norm ρ ℓ,o F . This is a somewhat involved technical problem, which is addressed in Sec. XI A and Appendix B. The basic idea is to utilize a cluster decomposition theorem at any finite time after the quench, see also Ref. [36]. The main difficulty in calculating the Frobenius norm of ρ ℓ,o is that the latter is not Gaussian. The idea is therefore to obtain ρ ℓ,o as a reduction of a Gaussian operator. To that end, we consider a composite system C = A ∪ n consisting of our subsystem A and a single site at position n, which is separated from A by a block B of length r, see Fig. 10.
The even part of the RDM ρ C (t) can be expanded in a complete basis of Hermitian involutions O e/o as (11.11) where · · · = Ψ t | · · · |Ψ t ≈ R ψ t | · · · |ψ t R , since both O e σ z n and O o σ α n in (11.11) are even operators with respect to fermion parity. In the limit of large separation r, we may use the cluster decomposition principle to simplify the expectation values This then leads to the following relation between RDMs in the limit or large separation lim r→∞ ρ C,e (t) = ρ ℓ,e (t)⊗ ρ 1,e (t)+ ρ ℓ,o (t)⊗ ρ 1,o (t), (11.13) where ρ 1 is the RDM of site n. The piece of interest to us is In the next step we move from spins to Majorana fermions by means of the Jordan-Wigner transformation (2.2) where p is a normalized, Gaussian operator (7.7) acting on the Hilbert space over sites A ∪ n p ≡ Tr A∪n e iπNB |ψ t RR ψ t | e iπNB . (11.17) In writing (11.15) we are assuming e iπNB = 0. The fact that p is Gaussian is a consequence of the particular form of |ψ t R (which is the analog of (2.13) in the R sector) and N B being quadratic in fermions. The odd part of the single-site RDM is of the form (11.18) and hence Here the late-time behaviour of m x,y (t) are given by (9.14) and (9.15) respectively, and following Ref. [31] they can be easily calculated numerically for all times. Combining (11.19) and (11.16) we obtain (11.20) Since both p and σ z n pσ z n are Gaussian, their moments can be written in terms of their respective correlation matrices We note that the correlation matrices are related byḠ = P n GP n , with P n the diagonal matrix that changes the sign of the last 2-by-2 block (I d is the d × d identity) While formally correct, (11.24) is not suitable for numerical computations, because at large distances e iπNB becomes very close to zero. A more convenient expression derived in Appendix B is (11.25) Here Γ A∪B∪n is the correlation matrix of the interval A ∪ B ∪ n and is given by (7.4), (7.5), (7.22). In order to utilize (11.25) we in principle have to consider infinite separations r and hence infinitely large matrices.
Crucially, in practice a finite separation r > 2v max t, where v max = max k ǫ h (k) is the maximal propagation velocity, is sufficient to recover the r → ∞ limit up to corrections that are exponentially small in r/ξ. Here ξ is the correlation length in the initial state. A representative example is shown in Fig. 11. In practice, using a finite r > 2v max t + ξδ with δ ≈ 20 provides an efficient way for calculating ||ρ ℓ,o (t)|| F and then by means of (11.10) distances D(ρ ℓ (t), ρ Ga ℓ ) for quenches originating in the ferromagnetic phase.

B. Results for quenches from the ferromagnetic phase
For quenches with h 0 < 1 and in the thermodynamic limit, we determine the distance between the quench RDM and that of an appropriate thermal or generalized Gibbs ensemble by means of relations (11.10) and (11.25). The correlation matrices for all cases are of the form (7.4), (7.5) with elements given in (7.12), (7.14) and (7.22) respectively. For a subsystem of size ℓ most terms require the calculation of determinants of 2ℓ × 2ℓ matrices, which is easily done numerically. The evaluation of ρ ℓ,0 F is significantly more costly, and in practice involves determinants of at most 2(ℓ + 2v max t + ξδ) × 2(ℓ + 2v max t + ξδ) matrices, as discussed above.
Results for a quench from h 0 = 1/3 to h = 2/3 and subsystem sizes ℓ = 10, 20, 30 . . . , 150 are shown in Figs 12 and 13. We see that the distance between quench and Distance D Gibbs = D(ρ ℓ (t), ρ β ℓ ) after a quench within the ferromagnetic phase for subsystem sizes ℓ = 10, 20, . . . , 150. As ℓ increases, the color fades from brown to green, the symbols become smaller and the curves narrower. At late times the distances tend to constants depending on subsystem size.
Gibbs RDMs tends to a ℓ-dependent constant at late times. On the other hand, as shown in Fig. 13, at sufficiently late times the distance between ρ ℓ (t) and ρ GGE,ℓ decays to zero in a universal power-law fashion The large-ℓ asymptotics of the function k(ℓ) can be inferred in the same way as for quenches within the paramagnetic phase. On surfaces with constant, small D, time scales as t ∼ ℓ 4/3 as is shown in Fig. 14, which implies that (11.27) We conclude that the late time behaviour of the distance between quench and generalized Gibbs RDMs is the same as for quenches within the paramagnetic phase. The mean relaxation time t * rms is therefore again given by (10.7). Interestingly this coincides with the result obtained in Ref. [32] for the relaxation of the order parameter two-point function after quenches within the ferromagnetic phase. The effects of the spontaneous symmetry breaking are important only at short and intermediate times. It is shown in the inset of Fig. 13 that there is a time window, in which the odd part of the RDM gives the dominant contribution to the distance, which decays exponentially.  The effects of the spontaneously broken Z 2 symmetry in the initial state make themselves felt through the Z 2 -odd part ρ ℓ,o of the density matrix. The relative importance of ρ ℓ,o for large ℓ can be estimated by considering the von Neumann entropy of subsystem A S vN [ρ ℓ ] = Tr [ρ ℓ ln (ρ ℓ )] . (11.28) We recall that the von Neumann entropy after a global quench grows linearly in time until the Fermi time t F = ℓ/(2v max ), and then saturates to a value proportional to the subsystem size ℓ 64,65 . Using the commutation relations (11.8) we see that the even part ρ ℓ,e can be expressed in terms of the full RDM ρ ℓ as follows ρ ℓ,e = ρ ℓ + e iπNA ρ ℓ e iπNA 2 . (11.29) Since for any set of density matrices ρ i the von Neumann entropy satisfies 66 30) the following bounds on the von Neumann entropy of subsystem A hold (11.31) This demonstrates that at any time after the quench the symmetry breaking contribution to the von Neumann entropy will be at most log 2. Given that for large subsystems the von Neumann entropy at late times is proportional to ℓ, we conclude that the relative contribution of the odd part of the RDM will be important only for small subsystem sizes.

A conjecture for ρ ℓ,o F in the limit of large ℓ and Jt
We now consider the space-time scaling limit 31 ℓ, Jt → ∞ , ℓ Jt fixed. (11.32) We observe that in this limit our numerical results for quenches within the ferromagnetic phase are in excellent agreement with the following relation Here we have highlighted the asymptotic nature of the relation and indicated by O(ℓ 0 , t 0 ), where the most im-portant corrections will arise. Since log ||ρ ℓ,e (t)|| F is pro-portional to the Rényi entropy S 2 (cf. Eq. (8.13)), we may use the known results 64 on the asymptotics of the latter log ρ ℓ,e (t) Combining (11.34) and (11.33) provides a conjecture for the asymptotic behaviour of ρ ℓ,o F . This conjecture is compared to numerical results in Fig. 15. The agreement is clearly quite good.

XII. QUENCHES ACROSS THE CRITICAL POINT
We now turn to quenches across the critical point. These are of particular interest 14,31,42 . In Fig. 16 we plot the distance between quench and GGE reduced density matrices for a quench from the ferromagnetic phase (h 0 = 1/2) to the paramagnetic phase (h = 3/2). The 15 data sets displayed correspond to subsystem sizes between ℓ = 10 and ℓ = 150. We find that the distance D GGE = D(ρ ℓ (t), ρ GGE ℓ ) again decays in a universal t −3/2 power law. In Fig. 16 we consider the same quench, but focus on very small subsystem sizes ℓ = 1, 2, 3, 4. We observe that the distance displays an oscillatory behaviour on top of a power-law decay in time. This is in agreement with the analytic results discussed in section IX B for the ℓ = 1 case. Increasing the subsystem size leads to a rapid suppression of the amplitude of the oscillations.
In Figs 18 and 19 we consider the reverse quenches, i.e. starting at h 0 = 3/2 in the paramagnetic phase, and quenching to h = 1/2 in the ferromagnetic phase. The behaviour of the distances is very similar to what we found for the quench from h 0 = 1/2 to h = 3/2: at late times the distance decays as a t −3/2 power law, and for small subsystem sizes we observe oscillatory behaviour on top of this decay.

XIII. SUMMARY AND CONCLUSIONS
In this work we have considered the evolution of reduced density matrices after a quantum quench in the transverse field Ising chain. The main result of our work is to demonstrate that where ρ ℓ (t) is the reduced density matrix of a subsystem consisting of ℓ adjacent spins after a quench of the transverse field, and ρ GGE,ℓ is the reduced density matrix of an appropriately defined generalized Gibbs ensemble. The derivation of (13.1) is based on defining an appropriate distance D(ρ, ρ ′ ) on the space of reduced density matrices, and then establishing that the distance between quench and GGE reduced density matrices approaches zero at late times. For our particular choice of distance we found that at late times this distance approaches zero as a universal power law in time We have presented a detailed construction of ρ GGE,ℓ in terms of the local (in space) integrals of motion I ± n of the TFIC. The densities of these conservation laws involve only spins on n + 2 consecutive sites. We proved that these local conservation laws are related in a linear fashion to the occupation numbers of the Bogoliubov fermions that diagonalize the Hamiltonian of the TFIC. This linear relation establishes the equivalence of our construction of the GGE to the one frequenctly used in the literature, which is based on mode occupation numbers.
We then have addressed the question, which of the conservation laws are most important for obtaining an accurate description of the stationary limit lim t→∞ ρ ℓ (t) of the quench RDM. To that end we introduced (defective) truncated generalized Gibbs ensembes, which are missing some of the local conservation laws. We found that the more local the conservation laws (i.e. the fewer consecutive spins their densities involve), the more important they are for describing the stationary state of a given subsystem. Loosely speaking we observed that in order to obtain a good description of the stationary state RDM of a subsystem of size ℓ, we need to retain all local conservation laws, whose densities involve at most ≈ ℓ + n 0 neighbouring spins, where n 0 is a constant depending on h 0 and h. Leaving out "highly local" conservation laws generally provides a very poor description of the stationary state. To the best of our knowledge this is the first such demonstration of a connection between locality of conservation laws and their importance in the GGE context.
Our work raises a number of issues. First and foremost is the dependence of the results obtained on the precise definition of the distance on the space of reduced density matrices. We have argued, that the "best" distance is the one based on the trace norm, because it provides the most direct and precise information on the time evolution of local observables. Unfortunately this distance is much harder to handle analytically. It would however be very interesting to implement it in purely numerical studies using iTEBD or related algorithms.
system. An upper bound is obtained as follows Here we have used that both ρ and ρ ′ are positive semidefinite and hence where 0 ≤ λ min ≤ λ j are the eigenvalues of ρ. To derive a lower bound we start by expressing the RDM of a block of ℓ spins in a spin-1 2 chain in the form where α i = 0, x, y, z with σ 0 ≡ I, and ρ is the density matrix of the full system; ρ ℓ is only function of the length because of translational invariance. By singling out the term with α ℓ = 0, we can express this in the form where ρ ℓ−1 is the RDM of the block consisting of sites 1, . . . , ℓ − 1. We also write the RDM of the ℓ th spin Tr δρ α ℓ ℓ−1 σ α ℓ ℓ (A5) and observe that Tr Ω α ℓ Here we have defined Ω α ℓ ℓ−1 = δρ α ℓ ℓ−1 − δρ ′α ℓ ℓ−1 . Using (A4) we have Tr Ω α ℓ ℓ−1 where we have used that for N × N matrices M we have N TrM 2 ≥ (TrM ) 2 in the second step, and (A6) in the last. Iterating Eq. (A7) ℓ − 1 times we obtain This implies that for sufficiently large subsystem size ℓ, the distance ρ ℓ − ρ ′ ℓ F will generally be larger than 2 1−ℓ/2 .
where G = P n GP n and P n is the diagonal involution defined in (11.22). We recall that {Q, Q 1 } denotes the product of the eigenvalues of (I + QQ 1 )/2 with halved degeneracy (the eigenvalues of QQ 1 are always double degenerate 67 for antisymmetric matrices Q and Q 1 ). The correlation matrix G defined in Eq. (11.21) turns out to be 59 the Schur complement of the block matrix Γ B of the matrix Γ A∪B∪n , i.e.
Here Γ R1, R2 denotes the matrix, whose rows and columns are associated with spatial regions R 1 and R 2 respectively, e.g.
Here ℓ is the size of subsystem A, while r is the distance between A and site n, see Fig. 10. The first step is to find determinant representations for {Q, Q} and {Q, P QP }, where P is a generic symmetric involution (P 2 = I and P t = P ).
We first consider {Q, Q}. Since Q is antisymmetric, its eigenvalues come in pairs ±q Both eigenvalues ±q give rise to the same eigenvalue 1+q 2 of I + Q 2 , and hence Here the product is over all positive eigenvalues of Q.
Using that Next we consider {Q, P QP }. (B11) where in the last step we used Eq. (B8). Since (P 1 2 ) 2 = P and P 2 = I, (B11) can be rewritten in the form Using (B8) and (B12), we can reexpress the quantities in (B2) as follows: Here we have use that the expectation value of the string operator in region B is related to the correlation matrix Γ B by e iπNB 2 = det |iΓ B |. A remaining problem is that lim r→∞ det |iΓ B | = 0, which precludes a numerical evaluation of (B1) on the basis of expressions (B13 (B16) The representations (B16) are suitable for numerical calculations even in the limit of large r. There is one further simplification: in the limit r → ∞ we have To see this, we expand the determinants in (B16) with respect to the last 2-by-2 block (from here on we omit the subscript in Γ A∪B∪n , i.e. Γ ≡ Γ A∪B∪n ) Using properties of the correlation matrix one could show that the determinants on the second line approach zero in the limit of large distance. For the sake of simplicity we propose a different proof, which is based on the assumption that the limit lim r→∞ det Γ A∪B + iI 2ℓ ⊕ 0 2r (B19) exists: we demonstrate that the limit cannot be infinite, so the expression in Eq. (B18) does tend to zero as r → ∞. To this end we consider the (2ℓ + 2r) × (2ℓ + 2r) correlation matrix G of a generic Gaussian density matrix, and show that the determinant det |G + iI 2ℓ ⊗ 0 2r | has an upper bound independent of r. Hence it cannot diverge in the limit r → ∞. Our proof is based on the following facts: (a.) G op ≤ 1, and hence G 2 op ≤ 1 and G + iI 2ℓ ⊗ 0 2r op ≤ G op +1 ≤ 2; (b.) G + iI 2ℓ ⊗ 0 2r cannot have more than 2ℓ eigenvalues with absolute values exceeding 1.
Property (a.) is a consequence of G being the correlation matrix of a positive semidefinite Gaussian. Property (b.) can be proved as follows: Let w a normalized vector with w i = 0 for any i ≤ 2ℓ. Then where the inequality follows from property (a.). If there were more than 2ℓ eigenvalues λ of G + iI 2ℓ ⊗ 0 2r with modulus larger than 1, we could find a linear combination W = i c i v i of the corresponding normalized eigenvectors v i with the property W i = 0 for any i ≤ 2ℓ; this leads to a contradiction with (B20) since This completes the proof of property (b.). Properties (a.) and (b.) imply that which establishes that the determinants in (B18) remain finite in the limit r → ∞. Concomitantly the expression in Eq. (B18) approaches zero as r → ∞. This establishes (B17). Putting everything together we see that (B16) can be written as which is Eq. (11.25).
We stress that our assumption reagrding the limit (B19) is equivalent to the existence of the limit in (B23). From a numerical point of view, this can be inferred from the scaling analysis of det I 2ℓ ⊕ 0 2r ⊕ I 2 + iΓ A∪B∪n which is still required to check the cluster decomposition hypothesis (see Fig. 11). The magnetization |m ⊥ (t)| can be computed writing a self-consistent equation for Eq. (B23) in the case ℓ = 1: From Eq. (11.19) we have which together with Eq. (B23) gives 4m 2 ⊥ (t) = lim r→∞ det I 2 ⊕ 0 2r ⊕ I 2 + iΓ 1 ∪B ∪n . (B26) Appendix C: Conservation laws in spin models with free fermion spectra In this appendix we present a simple construction of the bulk contribution to local conservation laws of the TFIC on the infinite line. Our method readily generalizes to other models with free fermionic spectrum such as the XY chain. Ignoring boundary conditions, we can use the Jordan-Wigner transformation to express the Hamiltonian as a quadratic form in Majorana fermions H = 1 2 l,n a l H ln a n .
Here H is a skewsymmetric block-circulant matrix where Y n = −Y T L−n are 2-by-2 matrices. In Fourier space we have where (Y k ) jn = −(Y −k ) nj . One can show that a complete set of local conservation laws is obtained by taking I r = 1 2 l,n a l I r;ln a n .
where ω We define the 'range' of a local conservation as the maximal number of neighbouring spins involved in its density minus one. By construction, the range is equal to the maximal |n| such thatȲ (r) n is nonzero (cf. Eqs (C1), (C2)). For the TFIC one finds that Y n = 0 for |n| > 1, and concomitantly the range of the Hamiltonian is r H = 1. It is straightforward to identify the conservation laws with ranges ≤ r + 1: from Eq. (C8) they are such that ω k = r+1 n=1 c − n sin(nk) , q k = r+1−rH n=0 c + n cos(nk) .
(C9) They can be divided in two classes: one with q k = 0, which we denote by I − , and one with ω k = 0, which we denote by I + . Finally, a complete set of conservation laws is given by These are exactly the conservation laws reported in Eq. (2.13). We note that the conservation laws I − r are independent of the system details, and can be found in any noninteracting model with a block circulant structure (see also Ref. [58]). Indeed they are originated from the trivial solution of Eq. (C6), namely the identity. In this appendix we discuss some properties of the defective generalized Gibbs ensembles defined in Section VI. We start by recalling the standard variational approach for deriving statistical ensembles in quantum mechanics. One generally seeks the density matrix that maximizes the entropy under a given set of constraints on independent, additive conservation laws I j δTr −ρ log ρ − λρ − j λ j I j ρ = 0 . (D1) The solution of (D1) is of the form ρ ∝ exp j λ j I j , which shows that the ensemble is a function only of the conservation laws appearing in Eq. (D1).
We now consider the density matrix after a quench. All the ensembles defined in the main text are compatible with the principle of maximal entanglement entropy, and the GGE, the truncated GGE, and the truncated defective GGE can be obtained (a posteriori) by means of the variational approach (D1).
Some complications arise when we consider defective GGEs, in which we exclude a single integral of motion. From Eq. (6.5) we find that the entanglement entropy density σ dGGE(+q) vN of the defective GGE ρ dGGE can be obtained from Eq. (D1). In the absence of peculiar constraints, one would expect the maximum to be generally a stationary point of the entanglement entropy. However, quenches in translationally invariant noninteracting models are very special since the initial state is a simultaneous eigenstate of an infinite number of local conservation laws. This substantially reduces the degrees of freedom, and can result in an exceptionally small domain for k + q (which may not include a stationary point). In Fig. 20 we show this paradoxical behaviour for the same set of parameters used in Fig. 9. Besides the pathological cases of even q, in which the curves collapse to the point κ + q = 0, the effect of the reduction of degrees of freedom is reflected in the "truncated" shape of the curves for q = 1, which turn out to as a function of the parameter κ + q for the same quench shown in Fig. 9 (the legend indicates the value of q). The points have the maximal entropy and correspond to the lines plotted in Fig. 9. Only for q = 1 the entanglement entropy is maximal at a stationary point.
be strictly decreasing functions of κ + q . The limiting procedure (6.2) selects the value of κ + q corresponding to the maximal entanglement entropy (the circles in Fig. 20).