Quench dynamics and relaxation in isolated integrable quantum spin chains

The previous argument shows that an isolated quantum system prepared in a pure state cannot relax in its entirety. However, it can and does relax locally in space [26, 28, 126, 155–157]. To explain this concept let us consider the example of a one-dimensional spin system with Hamiltonian H with only short-range interactions. An example would be the TFIC (6). We partition the entire system into an arbitrary but finite subsystem B and its complement A. Eventually we will take the thermodynamic limit while keeping B fixed. Let us prepare our spin system at time t  =  0 in some initial density matrix $\rho (0)$ that is not an eigenstate of H. At later times the entire system is characterized by the density matrix

$$\begin{eqnarray}&&\rho (t)={{\text{e}}^{-\text{i}Ht}}\rho (0){{\text{e}}^{\text{i}Ht}}.\end{eqnarray} \tag{ 18 }$$

The reduced density matrix of the subsystem B is obtained by tracing out the degrees of freedom in A

$$\begin{eqnarray}&&{{\rho}_{B}}(t)=\text{T}{{\text{r}}_{A}}\rho (t).\end{eqnarray} \tag{ 19 }$$

It is instructive to give an explicit representation of the reduced density matrix for our spin chain example. An orthonormal basis of states is given by

$$\begin{eqnarray}&&|{{\sigma}_{1}},{{\sigma}_{2}},\ldots,{{\sigma}_{L}}\rangle =\underset{j=1}{\overset{L}{{\otimes}}}|{{\sigma}_{j}}{{\rangle}_{j}},\quad {{\sigma}_{j}}=\pm 1,\end{eqnarray} \tag{ 20 }$$

where $|{{\sigma}_{j}}{{\rangle}_{j}}$ denote the two eigenstates of $S_{j}^{z}$. Let us take our subsystem B to consist of sites $1,2,\ldots \ell$ for simplicity. The reduced density matrix is then

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\rho}_{B}}(t) & = & \underset{{{\sigma}_{\ell +1}},\ldots,{{\sigma}_{L}}}{\sum}\left(\underset{k=\ell +1}{\overset{L}{{\otimes}}}{{}_{k}}\langle {{\sigma}_{k}}|\right)\rho (t)\left({{\underset{n=\ell +1}{\overset{L}{{\otimes}}}|{{\sigma}_{n}}\rangle}_{n}}\right) \\ {} & = & {{2}^{-\ell}}\underset{{{\alpha}_{1}},\ldots,{{\alpha}_{\ell}}=0,x,y,z}{\sum}\text{Tr}[\rho (t)\sigma _{1}^{{{\alpha}_{1}}}\sigma _{2}^{{{\alpha}_{2}}}\ldots \sigma _{\ell}^{{{\alpha}_{\ell}}}]~\sigma _{1}^{{{\alpha}_{1}}}\sigma _{2}^{{{\alpha}_{2}}}\ldots \sigma _{\ell}^{{{\alpha}_{\ell}}}, \end{array}\end{eqnarray} \tag{ 21 }$$

where we have defined $\sigma _{j}^{0}={{{\mathbb 1}}_{j}}$. Equation (21) shows that the matrix elements of ${{\rho}_{B}}(t)$ are equal to particular correlation functions of local operators acting non-trivially only on B, and that ${{\rho}_{B}}(t)$ in fact encodes all such correlators.

Definition 3. Local Relaxation.

We say that our system relaxes locally if the limit

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}\underset{L\to \infty}{{\lim}}{{\rho}_{B}}(t)={{\rho}_{B}}\left(\infty \right)\end{eqnarray} \tag{ 22 }$$

exists for any finite subsystem B (see figure 1).

Zoom In Zoom Out Reset image size

Definition 4. Stationary State.

Consider a system that relaxes locally in the sense just defined. Then its stationary state is defined as a time-independent density matrix ${{\rho}^{\text{SS}}}$ for the full system such that for any finite subsystem B

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\text{T}{{\text{r}}_{A}}\left({{\rho}^{\text{SS}}}\right)={{\rho}_{B}}\left(\infty \right),\end{eqnarray} \tag{ 23 }$$

where A is the complement of B. We stress that this equivalence applies only at the level of finite subsystems in the thermodynamic limit and not for the density matrices of the full system.

Definition 5. Local equivalence of ensembles.

Let ${{\rho}_{1}}$ and ${{\rho}_{2}}$ be two density matrices. We call the corresponding two ensembles locally equivalent, if in the thermodynamic limit the reduced density matrices for any finite subsystem B coincide, i.e.

$$\begin{eqnarray}&&\underset{|A|\to \infty}{{\lim}}\text{T}{{\text{r}}_{A}}\left({{\rho}_{1}}\right)=\underset{|A|\to \infty}{{\lim}}\text{T}{{\text{r}}_{A}}\left({{\rho}_{2}}\right).\end{eqnarray} \tag{ 24 }$$

Here A is the complement of B and $|A|$ denotes its volume. We denote local equivalence by

$$\begin{eqnarray}&&{{\rho}_{1}}{{=}_{\text{loc}}}{{\rho}_{2}}.\end{eqnarray} \tag{ 25 }$$

A key problem in the field of non-equilibrium dynamics is the determination of the stationary state density matrix. A crucial role is played by symmetries of the Hamiltonian. In translationally invariant systems there are essentially two paradigms for local relaxation, and we turn to their respective descriptions next. Before doing so an important comment is in order. By design our focus is on local operators, which in turn leads us to define local relaxation in the sense formulated above. In practice one might be interested in observables that are formally not captured by our framework, and yet might relax to stationary values. An example would be expectation values of the momentum distribution function, which is the Fourier transform of the single particle Green's function [24, 158]. The relaxational properties of such quantities are at present not understood in general, and further investigation is called for.

As we are dealing with an isolated system, energy is always a conserved quantity

$$\begin{eqnarray}&&E=\text{Tr}\left(\rho (t)H\right)=\text{Tr}\left({{\text{e}}^{-\text{i}Ht}}\rho (0){{\text{e}}^{\text{i}Ht}}H\right)=\text{Tr}\left(\rho (0)H\right).\end{eqnarray} \tag{ 26 }$$

In absence of other conserved quantities isolated systems are believed to locally relax to thermal equilibrium. This is known as thermalization [1, 159, 160]. In the framework we have set up, this means that our stationary state is described by a Gibbs ensemble

$$\begin{eqnarray}&&{{\rho}^{\text{SS}}}{{=}_{\text{loc}}}{{\rho}^{\text{Gibbs}}}=\frac{{{\text{e}}^{-{{\beta}_{\text{eff}}}H}}}{\text{Tr}\left({{\text{e}}^{-{{\beta}_{\text{eff}}}H}}\right)}.\end{eqnarray} \tag{ 27 }$$

Here the inverse effective temperature ${{\beta}_{\text{eff}}}$ is fixed by the initial value of the energy density

$$\begin{eqnarray}&&e\equiv \underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}\left(\rho (0)H\right)=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}\left({{\rho}^{\text{Gibbs}}}H\right).\end{eqnarray} \tag{ 28 }$$

We once again stress that (27) implies only that in the thermodynamic limit the stationary state is locally indistinguishable from a Gibbs ensemble in the sense that expectation values of all operators that act non-trivially only in a finite subsystem are identical to the corresponding thermal expectation values. The physical picture underlying thermalization is that the infinitely large complement of our subsystem acts like a heat bath with an effective inverse temperature ${{\beta}_{\text{eff}}}$.

We note that ${{\beta}_{\text{eff}}}$ defined in this way can be either positive or negative. The meaning of a negative temperature in this context is as follows. Let us consider a Hamiltonian with bounded spectrum (e.g. a quantum spin model). Then the average energy associated with H is

$$\begin{eqnarray}&&\bar{e}=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}(H).\end{eqnarray} \tag{ 29 }$$

Depending on our initial density matrix we now have two possibilities

$$\begin{eqnarray}\begin{array}{*{35}{l}} e & < & \bar{e}\Rightarrow {{\beta}_{\text{eff}}}>0, \\ e & > & \bar{e}\Rightarrow {{\beta}_{\text{eff}}}<0. \end{array}\end{eqnarray} \tag{ 30 }$$

We note that negative temperature ensembles have been observed experimentally in [161].

We now turn to the case of systems with additional local conservation laws I⁽ⁿ⁾. These are operators that commute with the Hamiltonian

$$\begin{eqnarray}&&[H,{{I}^{(n)}}]=0,\end{eqnarray} \tag{ 31 }$$

and have the property that their densities are local operators. In most of the cases of interest these conservation laws also commute with one another

$$\begin{eqnarray}&&[{{I}^{(n)}},{{I}^{(m)}}]=0,\end{eqnarray} \tag{ 32 }$$

and we will assume this to be the case in the following. We note that exceptions to (32) are known [51, 92, 93] (see also section 3.3.1). What we mean by local conservation laws is best explained by an example. The TFIC (6) has infinitely many local conservation laws in the thermodynamic limit [162]

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{I}^{(1,+)}} & =H(h), \\ {{I}^{(n,+)}} & =-J\underset{j}{\sum}(S_{j,j+n}^{xx}+S_{j,j+n-2}^{yy})+h(S_{j,j+n-1}^{xx}+S_{j,j+n-1}^{yy})\equiv \underset{j}{\sum}\mathcal{I}_{j}^{(n,+)}, \\ {{I}^{(n,-)}} & =-J\underset{j}{\sum}(S_{j,j+n}^{xy}-S_{j,j+n}^{yx})\equiv \underset{j}{\sum}\mathcal{I}_{j}^{(n,-)}, \end{array}\end{eqnarray} \tag{ 33 }$$

where $n\geqslant 1$ and we have defined

$$\begin{eqnarray}&&S_{j,j+\ell}^{\alpha \beta}=\sigma _{j}^{\alpha}\left[\underset{k=1}{\overset{\ell -1}{\prod}}\sigma _{j+k}^{z}\right]\sigma _{j+\ell}^{\beta},\quad S_{j,j}^{yy}=-\sigma _{j}^{z}.\end{eqnarray} \tag{ 34 }$$

The conservation laws themselves are extensive, but their densities $\mathcal{I}_{j}^{\left(n,\pm \right)}$ involve only finite numbers of neighbouring sites (n  +  1 sites for ${{I}^{\left(n,\pm \right)}}$). As a consequence of (31), expectation values of the conservation laws as well as their densities are time independent

$$\begin{eqnarray}&&\frac{1}{L}\text{Tr}(\rho (t){{I}^{(n,\pm )}})=\text{Tr}(\rho (t)\mathcal{I}_{j}^{(n,\pm )})=\frac{1}{L}\text{Tr}(\rho (0){{I}^{(n,\pm )}})\equiv \frac{{{E}^{(n,\pm )}}}{L},\end{eqnarray} \tag{ 35 }$$

where in the first step we used translational invariance. An obvious question at this point is why we insist on conservation laws to be local. The answer lies in the simple fact that any Hamiltonian has at least as many nonlocal conservation laws as there are basis states in the Hilbert space. Let us consider one-dimensional projectors on energy eigenstates

$$\begin{eqnarray}&&H|n\rangle ={{E}_{n}}|n\rangle,\quad {{P}_{n}}=|n\rangle \langle n|.\end{eqnarray} \tag{ 36 }$$

Then we clearly have

$$\begin{eqnarray}&&[{{P}_{n}},{{P}_{m}}]=0=[{{P}_{n}},H],\end{eqnarray} \tag{ 37 }$$

which imply that the expectation values of all P_n are time independent. Importantly, these conservation laws are not local, and as they exist for any Hamiltonian they are not expected to have any profound effect. In contrast, local conservation laws are rather special and, as is clear from (35), have important consequences for local relaxation after quantum quenches.

An immediate consequence of (35) is that systems with local conservation laws cannot thermalize after a quantum quench, because the system retains memory of the initial expectation values of all conserved quantities at all times. The maximum entropy principle [163] then suggests that the stationary state density matrix should be given by a GGE [164]

$$\begin{eqnarray}&&{{\rho}^{\text{GGE}}}=\frac{{{\text{e}}^{-\sum\nolimits_{n}{{\lambda}_{n}}{{I}^{(n)}}}}}{\text{Tr}\left({{\text{e}}^{-\sum\nolimits_{n}{{\lambda}_{n}}{{I}^{(n)}}}}\right)}.\end{eqnarray} \tag{ 38 }$$

Here ${{\lambda}_{n}}$ are Lagrange multipliers that are fixed by the initial conditions (35), as we require that

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\frac{{{E}^{(n,\pm )}}}{L}=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}({{\rho}^{\text{GGE}}}{{I}^{(n)}}).\end{eqnarray} \tag{ 39 }$$

3.3.1. Local conservation laws versus mode occupation numbers

In free theories GGEs are often formulated using conserved mode occupation numbers [164]. This is usually equivalent to our formulation in terms of local conservation laws in these cases [28] as we now demonstrate for a simple example. Let us consider a fermionic tight binding model

$$\begin{eqnarray}&&H=-J\underset{j}{\sum}c_{j}^{\dagger}{{c}_{j+1}}+c_{j+1}^{\dagger}{{c}_{j}}-\mu \underset{j}{\sum}c_{j}^{\dagger}{{c}_{j}},\quad \{{{c}_{j}},c_{n}^{\dagger}\}={{\delta}_{j,n}}.\end{eqnarray} \tag{ 40 }$$

The Hamiltonian is easily diagonalized by going to momentum space

$$\begin{eqnarray}&&H=\underset{k}{\sum}[-2J\cos (k)-\mu ]{{c}^{\dagger}}(k)c(k).\end{eqnarray} \tag{ 41 }$$

The mode occupation operators $n(k)={{c}^{\dagger}}(k)c(k)$ commute with the Hamiltonian and with one another. However, the n(k) are neither local nor extensive. An equivalent set of local, extensive conserved quantities is easily found

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{I}^{(n,+)}} & = & 2J\underset{k}{\sum}\cos (nk)~{{c}^{\dagger}}(k)c(k)=J\underset{j}{\sum}c_{j}^{\dagger}{{c}_{j+n}}+c_{j+n}^{\dagger}{{c}_{j}}, \\ {{I}^{(n,-)}} & = & 2J\underset{k}{\sum}\sin (nk)~{{c}^{\dagger}}(k)c(k)=\text{i}J\underset{j}{\sum}c_{j}^{\dagger}{{c}_{j+n}}-c_{j+n}^{\dagger}{{c}_{j}}. \end{array}\end{eqnarray} \tag{ 42 }$$

The crucial point is that these conservation laws are linearly related to the mode occupation numbers. This implies that the GGE describing the stationary state after a quantum quench to the Hamiltonian H can be constructed either from the local conservation laws, or from the mode occupation numbers

$$\begin{eqnarray}&&{{\rho}^{\text{GGE}}}=\frac{{{\text{e}}^{-\sum\nolimits_{n}{{\lambda}_{n,+}}{{I}^{\left(n,+\right)}}+{{\lambda}_{n,-}}{{I}^{\left(n,-\right)}}}}}{\text{Tr}\left({{\text{e}}^{-\sum\nolimits_{n}{{\lambda}_{n,+}}{{I}^{\left(n,+\right)}}+{{\lambda}_{n,-}}{{I}^{\left(n,-\right)}}}}\right)}=\frac{{{\text{e}}^{-\sum\nolimits_{k}{{\mu}_{k}}~n(k)}}}{\text{Tr}({{\text{e}}^{-\sum\nolimits_{k}{{\mu}_{k}}~n(k)}})}.\end{eqnarray} \tag{ 43 }$$

At this point a word of caution is in order. There are cases in which local conservation laws exist that cannot be expressed in terms of mode occupation operators. In non-interacting theories this occurs if the dispersion relation in the finite volume has additional symmetries. An example is provided by setting $\mu =0$ in our tight-binding model (40), in which case the dispersion has the symmetry $\epsilon \left(\,p\right)=-\epsilon \left(\pi -p\right)$. In this case the following operator commutes with the Hamiltonian

$$\begin{eqnarray}&&\underset{j}{\sum}{{(-1)}^{j}}[{{c}_{j}}{{c}_{j+1}}-c_{j}^{\dagger}c_{j+1}^{\dagger}]=\underset{k}{\sum}{{\text{e}}^{-\text{i}k}}c(\pi -k)c(k)+\text{h}.\text{c}.\qquad (\mu =0),\end{eqnarray} \tag{ 44 }$$

but is clearly neither expressible in terms of mode occupation operators, nor does it commute with the latter. In cases like this the stationary state is not always locally equivalent to a GGE of the form (43) since the constraints arising from the additional local charges must be imposed as well [51, 93, 165].

In equilibrium Statistical Physics it is well known that in the thermodynamic limit the Gibbs ensemble becomes equivalent to the microcanonical ensemble. In the quench context analogous equivalences of ensembles hold, and we now turn to their discussion.

3.4.1. Generic systems

We have stated above that generic systems thermalize when their stationary state density matrix is equal to a Gibbs ensemble in the sense of (23) and (27), where the effective temperature is fixed by the energy density e given by (28). There is strong numerical evidence [160, 166–175] that an equivalent microcanonical ensemble can be constructed in the form

$$\begin{eqnarray}&&{{\rho}^{\text{MC}}}=\,|n\rangle \langle n|,\end{eqnarray} \tag{ 45 }$$

where $|n\rangle$ is any energy eigenstate (of the post-quench Hamiltonian) with energy density e. Importantly, no averaging over a microcanonical energy shell is required. This implies that in generic systems all energy eigenstates at a given energy density are locally indistinguishable and thermal. The microcanonical description for the stationary state of a thermalizing system is then

$$\begin{eqnarray}&&{{\rho}^{\text{SS}}}{{=}_{\text{loc}}}{{\rho}^{\text{MC}}}.\end{eqnarray} \tag{ 46 }$$

3.4.2. Systems with local conservation laws

If a system has local conservation laws other than energy the above construction needs to be modified. This was first discussed in the non-interacting case by Cassidy et al [22] and subsequently generalized to interacting integrable models by Caux and Essler [176]. Let us recall that the 'initial data' for a (post quench) system with local conservation laws I ⁽ⁿ⁾ and Hamiltonian H  =  I ⁽¹⁾ is

$$\begin{eqnarray}&&{{\text{e}}^{(n)}}=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}(\rho (0){{I}^{(n)}}).\end{eqnarray} \tag{ 47 }$$

The stationary state is then locally equivalent to the density matrix

$$\begin{eqnarray}&&{{\rho}^{\text{SS}}}{{=}_{\text{loc}}}{{\rho}^{\text{GMC}}}=\,|\Phi\rangle \langle \Phi|,\end{eqnarray} \tag{ 48 }$$

where $|\Phi\rangle$ is a simultaneous eigenstate of all local conservation laws such that

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\left[\frac{1}{L}{{I}^{(n)}}-{{\text{e}}^{(n)}}\right]|\Phi\rangle =0.\end{eqnarray} \tag{ 49 }$$

Conditions (49) ensure that all I ⁽ⁿ⁾ have the correct expectation values (47) in the stationary state. An obvious question is how to construct $|\Phi\rangle$ in practice. For non-interacting lattice models this is often straightforward as we now discuss. In free theories the Hamiltonian is diagonalized in terms of mode occupation operators $n(k)={{c}^{\dagger}}(k)c(k)$

$$\begin{eqnarray}&&H=\underset{k}{\sum}\epsilon (k)n(k),\quad \left[n(k),n(q)\right]=0=\left[n(k),H\right].\end{eqnarray} \tag{ 50 }$$

As discussed above the n(k) are in one-to-one correspondence with local conservation laws and can therefore be used to construct GGEs. The initial conditions (47) therefore fix the densities $\rho (k)$ of particles with momentum k

$$\begin{eqnarray}&&\text{Tr}\left(\rho (0)n(k)\right)=\rho (k).\end{eqnarray} \tag{ 51 }$$

In a large, finite volume L, the state $|\Phi\rangle$ can then be taken as a Fock state

$$\begin{eqnarray}&&|\Phi\rangle =\underset{j=1}{\overset{N}{\prod}}{{c}^{\dagger}}\left({{k}_{j}}\right)|0\rangle,\quad c(q)|0\rangle =0.\end{eqnarray} \tag{ 52 }$$

The k_j need to be chosen to reproduce the correct particle density $\rho (k)$ in the thermodynamic limit $N,L\to \infty$, n  =  N/L fixed. Noting that $\Delta n=\rho (k)\Delta k$ we see that a possible choice is

$$\begin{eqnarray}&&{{k}_{j+1}}={{k}_{j}}+\frac{1}{L\rho \left({{k}_{j}}\right)}.\end{eqnarray} \tag{ 53 }$$

In practice we can determine $\rho (k)$ from (51), construct the state $|\Phi\rangle$ from (52) and (53), and then use it to form ${{\rho}^{\text{GMC}}}$. In interacting theories the construction is analogous but considerably more involved and will be discussed in section 7.

In contrast to the non-integrable case, energy eigenstates at a given energy density e are not all locally indistinguishable and thermal. While the typical state at a given e is thermal in this case as well, there exist other, non-thermal, macro-states at a given e. However, their entropies are smaller than the one of the thermal equilibrium state. A more detailed discussion is given in appendix B.

In the literature a different definition of relaxation after quantum quenches in finite systems is widely used [118, 121]. Let us for simplicity consider a one-dimensional system on a ring of length L, which is initially prepared in a state with density matrix $\rho (0)$. The time average of an operator $\mathcal{O}$ is then defined as

$$\begin{eqnarray}&&{{\bar{\mathcal{O}}}_{L}}=\underset{T\to \infty}{{\lim}}\frac{1}{T}{\int}_{0}^{T}\text{d}t~\text{Tr}\left(\rho (t)~\mathcal{O}\right),\end{eqnarray} \tag{ 54 }$$

where we have introduced the subscript L to indicate that the system size is kept fixed at L, and where we assumed the limit $T\to \infty$ to exist. The operator $\mathcal{O}$ is then said to relax if its expectation value is very close to its time average most of the time, i.e. if

$$\begin{eqnarray}&&\frac{1}{T}{\int}_{0}^{T}\text{d}t{{[\text{Tr}(\rho (t)~\mathcal{O})-{{\bar{\mathcal{O}}}_{L}}]}^{2}}\ll \bar{\mathcal{O}}_{L}^{2}.\end{eqnarray} \tag{ 55 }$$

Physically this way of thinking about relaxation is very different from ours. This is most easily understood for a non-interacting system, which by construction features stable particle and hole excitations. Let us denote their maximal group velocity by ${{v}_{\text{max}}}$. In time averaged relaxation one considers times t such that ${{v}_{\text{max}}}t\gg L$, i.e. the stable excitations traverse the entire system many times. In contrast, our way of defining relaxation is based on taking the thermodynamic limit first, and then considering late times.

A natural question to ask is what statistical ensemble describes the averages $\bar{\mathcal{O}}$. To see this let us consider time evolution with a Hamiltonian H that has non-degenerate eigenvalues E_n. Expanding the density matrix in the basis of energy eigenstates gives

$$\begin{eqnarray}&&\rho (t)=\underset{n,m}{\sum}\langle n|\rho (0)|m\rangle ~{{\text{e}}^{-\text{i}({{E}_{n}}-{{E}_{m}})t}}|n\rangle \langle m|.\end{eqnarray} \tag{ 56 }$$

Substituting this back into (54) and using that the energy eigenvalues are non-degenerate we arrive at

$$\begin{eqnarray}&&{{\bar{\mathcal{O}}}_{L}}=\underset{n}{\sum}\langle n|\rho (0)|n\rangle ~\langle n|\mathcal{O}|n\rangle.\end{eqnarray} \tag{ 57 }$$

This shows that time averages involve only the diagonal elements (in the energy eigenbasis) of the density matrix. The description (57) is known as diagonal ensemble. Is there a relation between the subsystem view of relaxation and the time-averaged one? It is believed that for local operators $\mathcal{O}$ the two viewpoints are in fact equivalent in the sense that

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}{{\bar{\mathcal{O}}}_{L}}=\text{Tr}[{{\rho}^{SS}}\mathcal{O}]\,.\end{eqnarray} \tag{ 58 }$$

In the case of the TFIC this has been shown in appendix E of [26] and elaborated on in [30].

An interesting question in the quantum quench context concerns problems where the Hamiltonian is invariant under a symmetry operation, but the initial state after the quench is not. If we denote the symmetry operation by U, then we have

$$\begin{eqnarray}&&[H,U]=0,\qquad U|\Psi(0)\rangle \ne |\Psi(0)\rangle.\end{eqnarray} \tag{ 59 }$$

The question is then, whether or not the symmetry associated with U is restored in the stationary state

$$\begin{eqnarray}&&[{{\rho}^{\text{SS}}},U]\overset{?}{{=}}0.\end{eqnarray} \tag{ 60 }$$

The answer to this is quite straightforward for the stationary states we have considered here, which are described by (generalized) Gibbs ensembles

$$\begin{eqnarray}&&{{\rho}^{\text{SS}}}{{=}_{\text{loc}}}{{\rho}_{\text{GGE}}}=\frac{1}{{{Z}_{\text{GGE}}}}{{\text{e}}^{-\underset{n}{\sum}{{\lambda}_{n}}{{I}^{(n)}}}}.\end{eqnarray} \tag{ 61 }$$

Clearly, if U is a symmetry of all local conservation laws

$$\begin{eqnarray}&&[{{I}^{(n)}},U]=0,\quad \forall n,\end{eqnarray} \tag{ 62 }$$

the symmetry associated with U will be restored in the stationary state. Conversely, if there is at least one conservation law for which $[{{I}^{(s)}},U]\ne 0$, then the symmetry will remain broken in the stationary state. It is useful to consider some explicit examples.

3.6.1. Spin-flip symmetry

The TFIC Hamiltonian (6) is invariant under rotations in spin space around the z-axis by 180 degrees

$$\begin{eqnarray}&&U\sigma _{\ell}^{\alpha}{{U}^{\dagger}}=-\sigma _{\ell}^{\alpha},\quad \alpha =x,y,\qquad U\sigma _{\ell}^{z}{{U}^{\dagger}}=\sigma _{\ell}^{z}.\end{eqnarray} \tag{ 63 }$$

A quantum quench of the transverse field starting in the ordered phase h  <  1 leads to an initial state that breaks this symmetry. However, the local conservation laws (33) that characterize the stationary state are invariant with respect to spin reflection. As a result the spin flip symmetry is restored in the stationary state [28], as will be shown in section 6.

3.6.2. Parity symmetry

The Hamiltonian of the tight-binding model (40) is invariant under the parity transformation (reflection across a link)

$$\begin{eqnarray}&&Uc_{\ell}^{\dagger}{{U}^{\dagger}}=c_{1-\ell}^{\dagger}.\end{eqnarray} \tag{ 64 }$$

On the other hand, the set of charges I ^(n,−) in (42) is odd under U, and hence a quantum quench from a state that is not parity even does not generally result in a parity-symmetric stationary state. As an example we may consider the time evolution of the ground state $|\Psi(0)\rangle$ of the following Hamiltonian

$$\begin{eqnarray}&&{{H}_{0}}=\underset{\ell}{\sum}\frac{\text{i}}{4}(3c_{\ell}^{\dagger}{{c}_{\ell +1}}-3c_{\ell +1}^{\dagger}{{c}_{\ell}}-c_{\ell}^{\dagger}c_{\ell +1}^{\dagger}+{{c}_{\ell +1}}{{c}_{\ell}})+{{\mu}_{0}}c_{\ell}^{\dagger}{{c}_{\ell}}.\end{eqnarray} \tag{ 65 }$$

For $|{{\mu}_{0}}|<\sqrt{2}$, $|\Psi(0)\rangle$ is not parity-symmetric and the expectation value of I^(n,−) is nonzero

$$\begin{eqnarray}&&\langle \Psi(0)|{{I}^{\left(n,-\right)}}|\Psi(0)\rangle =J\frac{1-{{(-1)}^{n}}}{2}\frac{\cos \left(n\arcsin \left(\frac{|{{\mu}_{0}}|}{\sqrt{2}}\right)\right)}{\pi n}\,\ne 0.\end{eqnarray} \tag{ 66 }$$

As I ^(n,−) are conserved we must have

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}({{\rho}^{\text{SS}}}~{{I}^{(n,-)}})=\underset{L\to \infty}{{\lim}}\frac{1}{L}\langle \Psi(0)|{{I}^{(n,-)}}|\Psi(0)\rangle,\end{eqnarray} \tag{ 67 }$$

and consequently the stationary state cannot be not symmetric under (64). We note that the above symmetry argument does not guarantee a non-zero value for the expectation value in (66), as there could be other symmetries that force it to vanish.

3.6.3. Translational symmetry

If we allow the initial state to partially break translational symmetry the situation can become significantly more complicated. An interesting example is provided by the quantum XY model [177]

$$\begin{eqnarray}&&{{H}_{\text{XY}}}=J\underset{\ell =1}{\overset{L}{\sum}}\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y},\end{eqnarray} \tag{ 68 }$$

where we take L to be even. As long as we consider quenches from translationally invariant initial states, the stationary behaviour is described by a GGE constructed from local conservation laws ${{Q}^{\left(n,\pm \right)}}$ along the lines described above. The situation changes for initial states $|{{\Psi}_{2}}(0)\rangle$ that are invariant only under translations by two sites. One might naively expect that translational symmetry gets restored in the stationary state, but this is in fact not the case [51]. This can be understood by noting that H has local conservation laws that are not translationally invariant, an example being

$$\begin{eqnarray}&&J_{1}^{+}=\underset{\ell =1}{\overset{L}{\sum}}{{(-1)}^{\ell}}\left[\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}-\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}\right],\quad \left[J_{1}^{+},{{H}_{\text{XY}}}\right]=0.\end{eqnarray} \tag{ 69 }$$

Generically one will have $\langle {{\Psi}_{2}}(0)|J_{1}^{+}|{{\Psi}_{2}}(0)\rangle \ne 0$, and as $J_{1}^{+}$ is conserved, translational invariance must remain broken at all times. The construction of a GGE in this case is complicated by the fact that $J_{1}^{+}$ does not commute with all ${{Q}^{\left(n,\pm \right)}}$, and involves an initial state dependent selection of mutually commuting local conservation laws.

As we have seen above, GGEs in integrable models involve infinite numbers of conservation laws in the thermodynamic limit. A natural question to ask is whether all of these are equally important for the description of the stationary state, or whether certain classes are more important than others. A general framework for investigating this question was developed in [28] and then applied to the specific example of quenches in the TFIC. The main findings are conjectured to apply much more generally to quenches in integrable models. A key measure of importance of a given conservation law is its degree of locality [28] D_loc. This is straightforward to quantify for conservation laws like (33), which have densities that act non-trivially only on a fixed number of n  +  1 neighbouring sites,

$$\begin{eqnarray}&&{{I}^{(n)}}=\underset{j}{\sum}\mathcal{I}_{j}^{(n)},\quad \text{range}(\mathcal{I}_{j}^{(n)})=n+1\qquad \Rightarrow {{D}_{\text{loc}}}({{I}^{(n)}})=n+1.\end{eqnarray} \tag{ 70 }$$

The idea of [28] was to select various finite subsets $\{{{J}_{m}}|m=1,\ldots,y\}\subset \{{{I}^{(n)}}\}$, and ask how well the density matrices

$$\begin{eqnarray}&&{{\rho}^{\text{pGGE},y}}=\frac{{{\text{e}}^{-\sum\nolimits_{n=1}^{y}{\lambda}_{n}^{\left(\,y\right)}{{J}_{n}}}}}{\text{Tr}\left({{\text{e}}^{-\sum\nolimits_{n=1}^{y}{\lambda}_{n}^{\left(\,y\right)}{{J}_{n}}}}\right)}\end{eqnarray} \tag{ 71 }$$

approximate the full GGE density matrix. Here the Lagrange multipliers ${\lambda}_{n}^{\left(\,y\right)}$ are fixed by imposing the appropriate initial conditions

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\frac{\langle \Psi(0)|{{J}_{n}}|\Psi(0)\rangle}{L}=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}\left({{\rho}^{\text{pGGE},y}}{{J}_{n}}\right),\quad n=1,\ldots,y.\end{eqnarray} \tag{ 72 }$$

A useful way of comparing the full ${{\rho}^{\text{GGE}}}$ and partial ${{\rho}^{\text{pGGE},y}}$ GGEs is through a distance $\mathcal{D}\left({{\rho}_{\ell}},{\rho}_{\ell}^{\prime}\right)$ on the space of reduced density matrices on an interval of length $\ell$. In order for ${{\rho}^{\text{pGGE},y}}$ to be a good approximation to ${{\rho}^{\text{GGE}}}$ it should be possible to make the distance between their respective reduced density matrices arbitrarily small

$$\begin{eqnarray}&&\mathcal{D}({\rho}_{\ell}^{\text{GGE}},{\rho}_{\ell}^{\text{pGGE},y})\overset{!}{{<}}\epsilon,\quad \forall \ell.\end{eqnarray} \tag{ 73 }$$

The findings of [28] show that if even a single conservation law with small degree of locality is excluded, the approximation (71) immediately becomes very poor, and the only way to achieve (73) is by retaining all conservation laws with the smallest degrees of locality (70), i.e. by considering 'truncated GGEs' of the form

$$\begin{eqnarray}&&{{\rho}^{\text{tGGE},\text{y}}}=\frac{{{\text{e}}^{-\sum\nolimits_{n=1}^{y}{\lambda}_{n}^{\left(\,y\right)}{{I}^{(n)}}}}}{\text{Tr}\left({{\text{e}}^{-\sum\nolimits_{n=1}^{y}{\lambda}_{n}^{\left(\,y\right)}{{I}^{(n)}}}}\right)}.\end{eqnarray} \tag{ 74 }$$

The value of y is then determined by the required degree of accuracy in (73). A consequence of these findings is that the full GGE can in fact be defined as the limit of truncated GGEs

$$\begin{eqnarray}&&{{\rho}^{\text{GGE}}}\equiv \underset{y\to \infty}{{\lim}}{{\rho}^{\text{tGGE},\text{y}}}.\end{eqnarray} \tag{ 75 }$$

The limiting procedure (75) was originally introduced in a non-interacting model (the transverse-field Ising) [28], and has since proved very useful for the construction of GGEs and the calculation of the stationary state values of local observables in interacting integrable models [17, 58, 59, 64, 178].

To exhibit the above ideas more explicitly it is useful to consider the example of a transverse field quench in the disordered phase of the TFIC [28]. Figure 2 shows results for an appropriately defined distance between the reduced density matrices of the full and truncated GGEs as a function of the number y of conservation laws retained in (74). Ten different subsystem sized $\ell =5,10,\ldots,50$ are shown.

Zoom In Zoom Out Reset image size

We see that for a given subsystem size $\ell$ the distance $\mathcal{D}({\rho}_{\ell}^{\text{GGE}},{\rho}_{\ell}^{\text{pGGE},\text{y}})$ starts decaying exponentially in y above an $\ell$-dependent value. This means that conservation laws I ⁽ⁿ⁾ with $n\gg \ell$ play a negligible role in describing stationary state properties in a subsystem of size $\ell$. The main message is then as follows:

The smaller the degree of locality of a conservation law is, the more important it is for describing the stationary properties of local observables.

It is clearly of interest to go beyond the equal time correlators we have discussed so far and consider dynamical (non equal time) correlations. In particular one can envisage using them to characterize the stationary state in the same way they are used in thermal equilibrium. An example would be the measurement of dynamical response functions in the stationary state by e.g. photoemission spectroscopy [179]. For such purposes the objects of interest are of the form

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}\langle \Psi(t)|{{\mathcal{O}}_{1}}\left({{t}_{1}}\right)\cdots {{\mathcal{O}}_{n}}\left({{t}_{n}}\right)|\Psi(t)\rangle,\qquad {{\mathcal{O}}_{j}}(t)={{\text{e}}^{\text{i}Ht}}{{\mathcal{O}}_{j}}{{\text{e}}^{-\text{i}Ht}},\end{eqnarray} \tag{ 76 }$$

where ${{\mathcal{O}}_{j}}$ are local observables. The problem one is faced with is that the descriptions of the stationary state by statistical ensembles (GGE, GMC) a priori hold only at the level of finite subsystems in the thermodynamic limit. On the other hand, time dependent operators act by construction non-trivially on the entire system, which moves them beyond the remit of applicability of the framework set out above. This issue was addressed in [29], which established that the stationary state density matrix ${{\rho}^{\text{SS}}}$ that describes the local relaxation in fact provides a correct description of dynamical correlations as well, i.e.

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}\langle \Psi(t)|{{\mathcal{O}}_{1}}\left({{t}_{1}}\right)\cdots {{\mathcal{O}}_{n}}\left({{t}_{n}}\right)|\Psi(t)\rangle =\text{Tr}\left[{{\rho}^{\text{SS}}}{{\mathcal{O}}_{1}}\left({{t}_{1}}\right)\cdots {{\mathcal{O}}_{n}}\left({{t}_{n}}\right)\right].\end{eqnarray} \tag{ 77 }$$

The proof of this statement is based on the Lieb–Robinson bound (see section 5.1) and more specifically on a theorem by Bravyi, Hastings, and Verstraete [180], who showed that time-evolving operators (in the Heisenberg picture) are well approximated by local operators with a range that increases linearly in time. We stress that (77) does not depend on which statistical ensemble describes the local relaxation in the stationary state, i.e. for generic systems ${{\rho}^{\text{SS}}}$ would be a Gibbs or a microcanonical density matrix, while for integrable systems it would be that of a GGE or GMC.

The issue of whether or not a fluctuation dissipation relation (FDR) holds in the stationary state was investigated in [29, 38, 41, 181]. In thermal equilibrium, the FDR connects the imaginary part of the linear response function ${{\chi}_{AB}}\left(\omega,\vec{q}|\rho \right)$ of two observables A, B to the corresponding spectral function ${{S}_{AB}}\left(\omega,\vec{q}|\rho \right)$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} -\frac{1}{\pi}\text{Im}~{{\chi}_{AB}}\left(\omega,\vec{q}~|{{\rho}^{\text{Gibbs}}}\left(\beta \right)\right)=\left(1-{{\text{e}}^{-\beta \omega}}\right){{S}_{AB}}\left(\omega,\vec{q}~|{{\rho}^{\text{Gibbs}}}\left(\beta \right)\right). \end{array}\end{eqnarray} \tag{ 78 }$$

Here we have defined

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\chi}_{AB}}(\omega,\vec{q}~|\rho ) & = & -\frac{\text{i}}{L}\underset{j,\ell}{\sum}{\int}_{0}^{\infty}\text{d}t~{{\text{e}}^{\text{i}\omega t-\text{i}\vec{q}\cdot ({{{\vec{r}}}_{j}}-{{{\vec{r}}}_{\ell}})}}\,\,\,\,\text{Tr}(\rho ~[{{A}_{j}}(t),{{B}_{\ell}}]), \\ {{S}_{AB}}(\omega,\vec{q}~|\rho ) & = & \frac{1}{L}\underset{j,\ell}{\sum}{\int}_{-\infty}^{\infty}\frac{\text{d}t}{2\pi}~{{\text{e}}^{\text{i}\omega t-\text{i}\vec{q}\cdot ({{{\vec{r}}}_{j}}-{{{\vec{r}}}_{\ell}})}}\,\,\,\,\,\text{Tr}(\rho ~{{A}_{j}}(t){{B}_{\ell}}). \end{array}\end{eqnarray} \tag{ 79 }$$

In the derivation of (78) one uses that in thermal equilibrium the spectral function for negative and positive frequencies are related in a simple way

$$\begin{eqnarray}&&{{S}_{BA}}(-\omega,-\vec{q}~|{{\rho}^{\text{Gibbs}}}(\beta ))={{\text{e}}^{-\beta \omega}}{{S}_{AB}}(\omega,\vec{q}~|{{\rho}^{\text{Gibbs}}}(\beta )).\end{eqnarray} \tag{ 80 }$$

Let us now turn to FDRs in the steady state after quantum quenches. Clearly, if the system thermalizes (27), the equilibrium FDR (78) with inverse temperature ${{\beta}_{\text{eff}}}$ applies. If on the other hand the system locally relaxes to a GGE, the imaginary part of the linear response function ceases to be proportional to the spectral function [29, 38, 41]. This can be traced back to the absence of a simple relationship between the spectral functions at positive and negative frequencies. However, the basic form of the FDR still holds [29]

$$\begin{eqnarray}&&-\frac{1}{\pi}\text{Im}~{{\chi}_{AB}}(\omega,\vec{q}~|{{\rho}^{\text{GGE}}})={{S}_{AB}}(\omega,\vec{q}~|{{\rho}^{\text{GGE}}})-{{S}_{BA}}(-\omega,-\vec{q}~|{{\rho}^{\text{GGE}}}).\end{eqnarray} \tag{ 81 }$$

According to our previous discussion the steady state should be described by the GGE (43), where the Lagrange multipliers are fixed by

$$\begin{eqnarray}&&\text{Tr}[{{\rho}^{\text{GGE}}}{{c}^{\dagger}}(k)c(k)]=\frac{1}{1+{{\text{e}}^{{{\mu}_{k}}}}}=\langle \Psi(0)|{{c}^{\dagger}}(k)c(k)|\Psi(0)\rangle ={{\sin}^{2}}({{\Theta}_{k}}/2).\end{eqnarray} \tag{ 94 }$$

As the GGE density matrix is Gaussian a Wick's theorem holds, and as a consequence of (94) the two-point functions in the GGE coincide with those of our steady state. As we have a Wick's theorem in the steady state as well, the GGE correctly reproduces all multi-point correlation functions. This proves that the steady state in our example is locally equivalent to the GGE (43) and (94).

As shown by Lieb and Robinson [193, 194], the velocity of information transfer in quantum spin chains is effectively bounded. More precisely, there exists a causal structure in commutators of local operators at different times

$$\begin{eqnarray}&&\Vert \left[{{\mathcal{O}}_{A}}(t),{{\mathcal{O}}_{B}}(0)\right]\Vert \leqslant c~\text{min}\left(|A|,|B|\right)~\Vert {{\mathcal{O}}_{A}}\Vert ~\Vert {{\mathcal{O}}_{B}}\Vert ~{{\text{e}}^{-\frac{L-vt}{\xi}}}.\end{eqnarray} \tag{ 101 }$$

Here ${{\mathcal{O}}_{A}}$ and ${{\mathcal{O}}_{B}}$ are local operators acting non-trivially only in two subsystems A and B that are spatially separated by a distance L, $\Vert.\Vert$ denotes the operator norm and $|A|$ the number of sites in subsystem A. Finally, c, v and ξ are constants. More recently, the Lieb–Robinson bounds have been refined [195, 196] and extended to mixed state dynamics in open quantum systems [196, 197].

The Lieb–Robinson bound has important consequences for quantum quenches starting in initial states with finite correlation lengths, and time evolving under a short-ranged Hamiltonian. It was shown in [180] that (101) implies a bound on the connected two-point correlation functions after such quenches

$$\begin{eqnarray}&&{{\langle \Psi(t)|{{\mathcal{O}}_{A}}~{{\mathcal{O}}_{B}}|\Psi(t)\rangle}_{\text{conn}}}<\bar{c}(|A|\,+\,|B|)~{{\text{e}}^{-\frac{L-2vt}{\chi}}}.\end{eqnarray} \tag{ 102 }$$

Here $\bar{c}$, v and χ are constants. The bound (102) shows that connected two-point functions of local operators after quantum quenches in spin chains are exponentially small up to times t  =  L/2v (see figure 4). This tallies very nicely with the light cone effects discussed above. We note that the bound does not provide values for the velocity v or the length χ.

Zoom In Zoom Out Reset image size

Throughout our discussion we have stressed that we are ultimately interested in taking the thermodynamic limit. In a large but finite system local observables can never truly relax. There always will be recurrences [198] such that the return amplitude $\mathcal{F}(t)=\,|\langle \Psi(0)|\Psi(t)\rangle |$ is arbitrarily close to 1

$$\begin{eqnarray}&&|1-\mathcal{F}(t)|<\epsilon.\end{eqnarray} \tag{ 103 }$$

However, in many-particle systems these typically occur only at astronomically late times. The exception to this rule are cases in which the spectrum of the post-quench Hamiltonian in a (large) finite volume for some reason has a highly commensurate structure. A very simple example are Hamiltonians with equidistant energy levels ${{E}_{n}}={{E}_{0}}+n\delta$, for which we have $\mathcal{F}\left(2\pi j/\delta \right)=1$.

A more relevant effect in practice are revivals, which refers to situations where $\mathcal{F}(t)=\mathcal{O}(1)$. This again requires the finite size energy spectrum to have a certain regularity, as is the case for example in conformal field theories [88, 91].

A different finite-size effect affects all observables and is related to the light cone repeatedly traversing the system [37, 50, 104, 106]. We refer to this effect as a traversal. As an example we consider the connected two point function (95) on a finite ring. A density plot is shown for a system of size L  =  100 in figure 5. We see that the light cone traverses the system and induces a signal in ${{S}_{L}}\left(\ell,t\right)$ at a time $\left(L\,-\,|\ell |\right)/2{{v}_{\text{max}}}$ after the light cone first reaches. Clearly at times $t>\left(L\,-\,|\ell |\right)/2{{v}_{\text{max}}}$ correlation functions of the finite system look very different from the ones in the thermodynamic limit. As is shown in figure 6, the traversal in not associated with any revival, because the return amplitude remains exponentially small in the system size.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Spin chains that can be mapped to free fermions differ in two important aspects from free fermion models like the one we considered in section 4:

• two-point functions of spin operators map onto n-point correlation functions of fermions, where n is related to the distance between the two spins and can be arbitrarily large;
• The ground state of the TFIC in the ordered phase is not a Fock state (as a result of spontaneous symmetry breaking).

The Hamiltonian (104) can be mapped to a fermionic theory by a Jordan–Wigner transformation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \sigma _{\ell}^{z}=i{{a}_{2\ell}}{{a}_{2\ell -1}},\quad \sigma _{\ell}^{x}=\left(\underset{j=1}{\overset{\ell -1}{\prod}}\left(i{{a}_{2j}}{{a}_{2j-1}}\right)\right){{a}_{2\ell -1}},\quad \sigma _{\ell}^{y}=\left(\underset{j=1}{\overset{\ell -1}{\prod}}\left(i{{a}_{2j}}{{a}_{2j-1}}\right)\right){{a}_{2\ell}}, \end{array}\end{eqnarray} \tag{ 106 }$$

Here ${{a}_{\ell}}$ are Majorana fermions satisfying the anti-commutation relations $\left\{{{a}_{\ell}},{{a}_{n}}\right\}=2{{\delta}_{\ell n}}$. The usual spinless fermions are obtained by taking linear combinations $c_{\ell}^{\dagger}=\left({{a}_{2\ell -1}}+\text{i}{{a}_{2\ell}}\right)/2$. It is now straightforward to see that spin–spin correlation functions map onto expectation values of strings of fermions, e.g.

$$\begin{eqnarray}&&\langle \sigma _{\ell}^{x}\rangle ={{(-i)}^{\ell -1}}\left\langle \underset{j=1}{\overset{2\ell -1}{\prod}}{{a}_{j}}\right\rangle \qquad {{\langle \sigma _{\ell}^{x}{{\sigma}_{\ell +n}}\rangle}^{x}}={{(-i)}^{n}}\left\langle \underset{j=2\ell}{\overset{2\ell +2n-1}{\prod}}{{a}_{j}}\right\rangle.\end{eqnarray} \tag{ 107 }$$

Application of the Jordan–Wigner transformation to the TFIC Hamiltonian (104) results in a fermion Hamiltonian of the form

$$\begin{eqnarray}\begin{array}{r l l} H & = & \frac{\text{I}-{{\text{e}}^{\text{i}\pi \mathcal{N}}}}{2}{{H}_{\text{R}}}+\frac{\text{I}+{{\text{e}}^{\text{i}\pi \mathcal{N}}}}{2}{{H}_{\text{NS}}}, \\ {{H}_{\text{NS}/\text{R}}} & = & \text{i}J\underset{j=1}{\overset{L-1}{\sum}}{{a}_{2j}}\left[{{a}_{2j+1}}-h{{a}_{2j-1}}\right]-\text{i}J{{a}_{2L}}\left[h{{a}_{2L-1}}\mp {{a}_{1}}\right]. \end{array}\end{eqnarray} \tag{ 108 }$$

Here ${{\text{e}}^{\text{i}\pi \mathcal{N}}}$ is the fermion parity operator with eigenvalues $\pm 1$

$$\begin{eqnarray}&&{{\text{e}}^{\text{i}\pi \mathcal{N}}}=\underset{\ell =1}{\overset{L}{\prod}}\sigma _{\ell}^{z}={{(-i)}^{L}}\underset{j=1}{\overset{2L}{\prod}}{{a}_{j}},\qquad {{\text{e}}^{\text{i}\pi \mathcal{N}}}{{a}_{j}}=-{{a}_{j}}{{\text{e}}^{\text{i}\pi \mathcal{N}}}.\end{eqnarray} \tag{ 109 }$$

${{H}_{\text{R},\text{NS}}}$ commute with the fermion parity operator, and the full Hamiltonian (108) is therefore block-diagonal: H_R (H_NS) describes the action on states with an odd (even) number of fermions. We note that the free fermion Hamiltonians ${{H}_{\text{NS}/\text{R}}}$ are closely related to the pairing model (83) considered earlier

$$\begin{eqnarray}&&{{H}_{\text{R}}}=-H\left(-J,-2hJ\right)+JhL.\end{eqnarray} \tag{ 110 }$$

6.1.1. Ground states

The Hamiltonians ${{H}_{\text{NS}/\text{R}}}$ can be diagonalized by Bogoliubov transformations to canonical momentum space fermion operators b_p (details can be found in e.g. appendix A of [25])

$$\begin{eqnarray}&&{{H}_{\text{a}}}(h)=\underset{p\in \text{a}}{\sum}{{\varepsilon}_{h}}\left(\,p\right)\left(b_{p}^{\dagger}{{b}_{p}}-\frac{1}{2}\right),\quad \text{a}=\text{R},~\text{NS},\end{eqnarray} \tag{ 111 }$$

where the single-particle energy is given by

$$\begin{eqnarray}&&{{\varepsilon}_{h}}(k)=2J\sqrt{1+{{h}^{2}}-2h\cos k}.\end{eqnarray} \tag{ 112 }$$

The difference between R and NS sectors enters via the allowed values of the momenta, which are $p=\frac{\pi n}{L}$, where n are even/odd integers for R and NS fermions respectively. The ground states of ${{H}_{\text{R},\text{NS}}}(h)$ are the fermionic vacua

$$\begin{eqnarray}&&{{b}_{p}}|GS{{\rangle}_{\text{a}}}=0~\forall p\in \text{a},\quad \text{a}=\text{R},\text{NS}.\end{eqnarray} \tag{ 113 }$$

These vacuum states are also eigenstates of the fermion parity operator

$$\begin{eqnarray}&&{{\text{e}}^{\text{i}\pi \mathcal{N}}}|GS{{\rangle}_{\text{NS}}}=\,|GS{{\rangle}_{\text{NS}}},\qquad {{\text{e}}^{\text{i}\pi \mathcal{N}}}|GS{{\rangle}_{\text{R}}}=\text{sgn}(h-1)|GS{{\rangle}_{\text{R}}}.\end{eqnarray} \tag{ 114 }$$

From (108) it follows that in the ferromagnetic phase h  <  1 both fermion vacua are eigenstates of the full Hamiltonian H. Their respective energies are exponentially (in system size) close, and they become degenerate in the thermodynamic limit. Spin-flip symmetry then gets spontaneously broken, and the ground state is either the symmetric or the antisymmetric combination of the two vacuuum states. In the paramagnetic phase h  >  1 the ground state of H is given by the NS vacuum state. In summary, we have

$$\begin{eqnarray}|GS\rangle =\left\{\begin{array}{*{35}{l}} \frac{|GS{{\rangle}_{\text{NS}}}\pm |GS{{\rangle}_{\text{R}}}}{\sqrt{2}} & h<1 \\ |GS{{\rangle}_{\text{NS}}} & h>1. \end{array}\right. \end{eqnarray} \tag{ 115 }$$

This shows that for h  <  1 the ground state of H is not a Fock state.

We now consider the following quench protocol. We prepare the system in the ground state of H(h₀), and at time t  =  0 quench the transverse field to a new value h. At times t  >  0 we time evolve with the new Hamiltonian H(h). All local operators in spin basis can be classified according to their fermion parity, and this turns out to be very useful. We have

$$\begin{eqnarray}&&{{\text{e}}^{\text{i}\pi \mathcal{N}}}{{\mathcal{O}}_{e/o}}{{\text{e}}^{-\text{i}\pi \mathcal{N}}}=\pm {{\mathcal{O}}_{e/o}}.\end{eqnarray} \tag{ 116 }$$

6.2.1. Quenches originating in the paramagnetic phase

For a quench starting in the paramagnetic phase it follows from (115), (108) and (114) that the state of the system at times t  >  0 is given by

$$\begin{eqnarray}&&|\Psi(t)\rangle ={{\text{e}}^{-\text{i}Ht}}|GS\rangle \equiv {{\text{e}}^{-\text{i}{{H}_{\text{NS}}}t}}|GS{{\rangle}_{\text{NS}}}.\end{eqnarray} \tag{ 117 }$$

This is even under fermion parity. Hence expectation values of odd operators must vanish

$$\begin{eqnarray}&&\langle \Psi(t)|{{O}_{o}}|\Psi(t)\rangle =0.\end{eqnarray} \tag{ 118 }$$

By virtue of the simple form of both the time evolution operator H_NS and the initial state $|\text{GS}{{\rangle}_{\text{NS}}}$ in (117), expectation values of even operators can be calculated by applying Wick's theorem. This allows the reduction of expectation values of strings of fermion operators to Pfaffians involving only the two-point functions $\langle {{a}_{i}}{{a}_{j}}\rangle$ [46, 47, 177]. For example we have

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle \sigma _{\ell}^{x}\sigma _{\ell +2}^{x}\rangle & = & -\langle {{a}_{2\ell}}{{a}_{2\ell +1}}{{a}_{2\ell +2}}{{a}_{2\ell +3}}\rangle \\ {} & = & -\langle {{a}_{2\ell}}{{a}_{2\ell +1}}\rangle \langle {{a}_{2\ell +2}}{{a}_{2\ell +3}}\rangle +\langle {{a}_{2\ell}}{{a}_{2\ell +2}}\rangle \langle {{a}_{2\ell +1}}{{a}_{2\ell +3}}\rangle -\langle {{a}_{2\ell}}{{a}_{2\ell +3}}\rangle \langle {{a}_{2\ell +1}}{{a}_{2\ell +2}}\rangle. \end{array}\end{eqnarray} \tag{ 119 }$$

6.2.2. Quenches originating in the ferromagnetic phase

As a result of spontaneous symmetry breaking in the ground state, the situation for quenches originating in the ferromagnetic phase h₀  <  1 is more complicated. The time evolved initial state is

$$\begin{eqnarray}&&|\Psi(t)\rangle =\frac{{{\text{e}}^{-\text{i}{{H}_{\text{NS}}}t}}|GS{{\rangle}_{\text{NS}}}\pm {{\text{e}}^{-\text{i}{{H}_{\text{R}}}t}}|GS{{\rangle}_{\text{R}}}}{\sqrt{2}}.\end{eqnarray} \tag{ 120 }$$

The expectation values of even operators can be expressed in the form

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle \Psi(t)|{{\mathcal{O}}_{e}}|\Psi(t)\rangle & = & \frac{1}{2}\underset{\text{a}\in \left\{\text{R},\text{NS}\right\}}{\sum}{{}_{\text{a}}}{{\langle GS|{{\text{e}}^{\text{i}{{H}_{\text{NS}}}t}}{{\mathcal{O}}_{e}}{{\text{e}}^{-\text{i}{{H}_{\text{NS}}}t}}|GS\rangle}_{\text{a}}} \\ {} & \underset{L\to \infty}{{\longrightarrow}} & {{}_{\text{NS}}}{{\langle GS|{{\text{e}}^{\text{i}{{H}_{\text{NS}}}t}}{{\mathcal{O}}_{e}}{{\text{e}}^{-\text{i}{{H}_{\text{NS}}}t}}|GS\rangle}_{\text{NS}}}. \end{array}\end{eqnarray} \tag{ 121 }$$

In the last line we have used that the expectation values in the R and NS sectors become equal in the thermodynamic limit. The last line in (121) can again be evaluated by application of Wick's theorem.

In contrast, expectation of odd operators cannot be simplified in this way. Instead we have

$$\begin{eqnarray}&&\langle GS|{{\text{e}}^{\text{i}Ht}}{{\mathcal{O}}_{o}}{{\text{e}}^{-\text{i}Ht}}|GS\rangle =\pm \text{Re}\left[{{}_{\text{NS}}}{{\langle GS|{{\text{e}}^{\text{i}{{H}_{\text{NS}}}t}}{{\mathcal{O}}_{o}}{{\text{e}}^{-\text{i}{{H}_{\text{R}}}t}}|GS\rangle}_{\text{R}}}\right].\end{eqnarray} \tag{ 122 }$$

Wick's theorem does not apply here, and in order to proceed one commonly resorts to one of the following methods [202]:

• 1.  
  Using the cluster decomposition property we can obtain (122) from the expectation value of an even operator by considering the limit

  $$\begin{eqnarray}&&|\langle GS|{{\text{e}}^{\text{i}Ht}}{{\mathcal{O}}_{o}}{{\text{e}}^{-\text{i}Ht}}|GS\rangle |\,=\underset{d\to \infty}{{\lim}}\sqrt{{{}_{\text{NS}}}{{\langle GS|{{\text{e}}^{\text{i}{{H}_{\text{NS}}}t}}{{\mathcal{O}}_{o}}(r){{\mathcal{O}}_{o}}(r+d){{\text{e}}^{-\text{i}{{H}_{\text{NS}}}t}}|GS\rangle}_{\text{NS}}}}.\end{eqnarray} \tag{ 123 }$$
• 2.  
  One imposes open boundary conditions on the spins (in that case the Hamiltonian is mapped into a purely quadratic form of fermions) and considers the expectation value of ${{\mathcal{O}}_{o}}$ asymptotically far away from the boundaries.

An alternative method that applies more generally to integrable models was developed in [25–27]. It is based on the observation that the initial state after the quench can be expressed in a squeezed state form, e.g.

$$\begin{eqnarray}&&\langle\Psi(0)=\exp \left[\text{i}\underset{0 < p\in \text{NS}}{\sum}K(\,p)b_{p}^{\dagger}b_{-p}^{\dagger}\right]G{{S}\rangle_{\text{NS}}}.\end{eqnarray} \tag{ 124 }$$

Here K(p) is a function that depends on the quench (see (132)), $b_{p}^{\dagger}$ are the aforementioned momentum space Bogoliubov fermion operators. Two point correlation functions can then be written in a Lehmann representation based on energy eigenstates. The matrix elements ('form factors') in the Lehmann representation are known [203], and it is possible to obtain explicit results in the framework of a low-density expansion. We refer the reader to [26, 27] for further details and the explicit calculations. Here we only review and discuss the main results.

Stationary state properties were analyzed in detail in [27]. An important simplification at late times is that expectation values of odd operators go to zero. As discussed above, for quenches originating in the paramagnetic phase they vanish identically because the ${\mathbb{Z}_{2}}$ symmetry remains unbroken. On the other hand, for quenches originating in ferromagnetic phase (h₀  <  1) their expectation value is generally nonzero (see (122)). However, as shown in [28], expectation values of all odd local operators decay exponentially in time to zero. This leaves us with expectation values for even operators, which can be analyzed by standard free fermion methods. The basic object is the fermion two-point function, which in the thermodynamic limit can be written in the form (see (91))

$$\begin{eqnarray}&&\left. \text{i}\langle {{a}_{\ell}}{{a}_{n}}\rangle \right|_{\ell \ne n}={\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}[{{A}_{\ell -n}}(k)+{{B}_{\ell -n}}(k){{\text{e}}^{2\text{i}{{\varepsilon}_{h}}(k)t}}+B_{\ell -n}^{\ast}(k){{\text{e}}^{-2\text{i}{{\varepsilon}_{h}}(k)t}}].\end{eqnarray} \tag{ 125 }$$

Here A(k) and B(k) are smooth functions (we are considering noncritical quenches) that depend on the quench details and ${{\varepsilon}_{h}}(k)$ is the dispersion relation (112). By the Riemann-Lebesgue lemma, in the infinite time limit the fermion two-point functions approach stationary values

$$\begin{eqnarray}&&\left. \underset{t\to \infty}{{\lim}}\text{i}\langle {{a}_{\ell}}{{a}_{n}}\rangle \right|_{\ell \ne n}={\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}{{A}_{\ell -n}}(k).\end{eqnarray} \tag{ 126 }$$

Since the expectation value of any even local operator can always be written as a finite sum of finite products of fermion two-point functions (see (119)), the infinite time limit exists and is obtained by replacing (125) with (126). Courtesy of Wick's theorem, the stationary properties of all other even operators can be expressed in terms of (126).

6.3.1. Description of the steady state by a GGE

Since expectation values of odd operators vanish at infinite times, the steady state can be constructed in complete analogy to our fermionic example considered in section 4. The appropriate GGE density matrix is of the form

$$\begin{eqnarray}&&{{\rho}_{\text{GGE}}}=\frac{{{\text{e}}^{-\sum\nolimits_{j=1}^{\infty}{\lambda}_{n}^{+}{{I}^{\left(n,+\right)}}+{\lambda}_{n}^{-}{{I}^{\left(n,-\right)}}}}}{{{Z}_{\text{GGE}}}},\end{eqnarray} \tag{ 127 }$$

where the conservation laws ${{I}^{\left(n,\pm \right)}}$ have been reported earlier in (33). The Lagrange multipliers for transverse field quenches were determined in [28]

$$\begin{eqnarray}&&{\lambda}_{j}^{+}=\left\{\frac{\text{sgn}\left[(h-1)\left({{h}_{0}}-1\right)\right]}{{{\varepsilon}_{h}}(0)}+\frac{{{(-1)}^{j}}}{{{\varepsilon}_{h}}\left(\pi \right)}\right\}\frac{2}{j},\qquad {\lambda}_{j}^{-}=0.\end{eqnarray} \tag{ 128 }$$

The local equivalence of the steady state density matrix to this GGE was demonstrated in [25, 28]

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}|\Psi(t)\rangle \langle \Psi(t)|{{=}_{\text{loc}}}{{\rho}_{\text{GGE}}}.\end{eqnarray} \tag{ 129 }$$

The GGE can be interpreted as a Gibbs ensemble at inverse temperature $\beta ={{J}^{-1}}$ for an effective 'GGE Hamiltonian', defined as

$$\begin{eqnarray}&&{{H}_{\text{GGE}}}\equiv J\underset{j=1}{\overset{\infty}{\sum}}{\lambda}_{j}^{+}I_{j}^{+}.\end{eqnarray} \tag{ 130 }$$

By virtue of (128) this 'Hamiltonian' is long-ranged. It can be diagonalized by combined Jordan–Wigner and Bogoliubov transformations [28], which for a quench from h₀ to h results in

$$\begin{eqnarray}&&{{H}_{\text{GGE}}}={\int}_{-\pi}^{\pi}\frac{\text{d}p}{2\pi}J\log ({{K}^{2}}(\,p))\left[{{b}^{\dagger}}(\,p)b(\,p)-\frac{1}{2}\right],\end{eqnarray} \tag{ 131 }$$

where

$$\begin{eqnarray}&&{{K}^{2}}\left(\,p\right)=\frac{{{\sin}^{2}}\left(\,p\right)~{{\left(h-{{h}_{0}}\right)}^{2}}}{{{\left({{\varepsilon}_{{{h}_{0}}}}\left(\,p\right){{\varepsilon}_{h}}\left(\,p\right){{(2J)}^{-2}}+1+h{{h}_{0}}-\left(h+{{h}_{0}}\right)\cos \left(\,p\right)\right)}^{2}}}.\end{eqnarray} \tag{ 132 }$$

The 'dispersion relation' $J\log {{K}^{2}}\left(\,p\right)$ diverges logarithmically at momenta zero and π. This is related to the fact that the mode occupation numbers at these momenta are independent of h, and ultimately produces the algebraic decay (128) of the Lagrange multipliers. These logarithmic singularities do not compromise the cluster decomposition properties of the steady state.

6.3.2. Connected spin–spin correlation functions

In the stationary state the connected two-point correlators decay exponentially with distance [26]

$$\begin{eqnarray}&&{\rho}_{c}^{\alpha \alpha}\left(\ell \right)=\langle \sigma _{j}^{\alpha}\sigma _{j+\ell}^{\alpha}\rangle -\langle \sigma _{j}^{\alpha}\rangle \langle \sigma _{j+\ell}^{\alpha}\rangle \simeq {{C}^{\alpha}}\left(\ell \right){{\text{e}}^{-\ell /{{\xi}_{\alpha}}}}.\end{eqnarray} \tag{ 133 }$$

The correlation lengths are given by [26]

$$\begin{eqnarray}\begin{array}{*{35}{l}} \xi _{z}^{-1} & = & |\ln {{h}_{0}}|\,+\,\min \left(|\ln {{h}_{0}}|,|\ln h|\right), \end{array}\end{eqnarray} \tag{ 134 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \xi _{x}^{-1} & = & {{\theta}_{\text{H}}}(h-1){{\theta}_{\text{H}}}\left({{h}_{0}}-1\right)\ln \left[\min \left({{h}_{0}},{{h}_{1}}\right)\right]-\ln \left[{{x}_{+}}+{{x}_{-}}+{{\theta}_{\text{H}}}\left((h-1)\left({{h}_{0}}-1\right)\right)\sqrt{4{{x}_{+}}{{x}_{-}}}\right], \end{array}\end{eqnarray} \tag{ 135 }$$

where ${{\theta}_{\text{H}}}(x)$ is the Heaviside step function and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{x}_{\pm}}=\frac{\left[\text{min}\left(h,{{h}^{-1}}\right)\pm 1\right]\left[\text{min}\left({{h}_{0}},h_{0}^{-1}\right)\pm 1\right]}{4},\quad {{h}_{1}}=\frac{1+h{{h}_{0}}+\sqrt{\left({{h}^{2}}-1\right)\left(h_{0}^{2}-1\right)}}{h+{{h}_{0}}}. \end{array}\end{eqnarray} \tag{ 136 }$$

The prefactors ${{C}^{\alpha}}\left(\ell \right)$ depend on the details of the quench and are listed in appendix C.

Having established the stationary behaviour of spin correlations, we now turn to their dynamics at late times. The first question of interest is how they relax towards their stationary values.

For even operators we can use Wick's theorem to express spin–spin correlators in terms of the fermion two-point functions (125). At late times the resulting expression can be evaluated by a stationary phase approximation, which gives

$$\begin{eqnarray}&&\langle {{\mathcal{O}}_{e}}\rangle \simeq {{\langle {{\mathcal{O}}_{e}}\rangle}_{\text{GGE}}}+O\left({{t}^{-n\left[{{\mathcal{O}}_{e}}\right]/2}}\right).\end{eqnarray} \tag{ 137 }$$

Here the exponent $n\left[{{\mathcal{O}}_{e}}\right]$ is an integer that depends on the particular even operator under consideration.

Expectation values of odd operators were argued in [28] to decay exponentially in time

$$\begin{eqnarray}&&\langle {{\mathcal{O}}_{o}}\rangle \simeq O\left({{\text{e}}^{-t/\tau \left[{{\mathcal{O}}_{o}}\right]}}\right),\end{eqnarray} \tag{ 138 }$$

where $\tau \left[{{\mathcal{O}}_{o}}\right]$ denotes a relaxation time. Having established the gross structure of the late time dynamics, we now turn to a more quantitative description.

6.4.1. One-point functions

Longitudinal spin operator

The longitudinal spin operator $\sigma _{j}^{x}={{(-i)}^{\ell -1}}\prod\nolimits_{j=1}^{2\ell -1}{{a}_{j}}$ is the simplest and most important example of an odd operator. Its expectation value is the order parameter in the ferromagnetic phase. As we are dealing with an odd operator, its expectation value is identically zero for quenches originating in the paramagnetic phase. For quenches from the ferromagnetic phase (h₀  <  1) it was shown in [25, 26] that

$$\begin{eqnarray}&&\langle \Psi(t)|\sigma _{j}^{x}|\Psi(t)\rangle \simeq {{C}^{x}}(t){{\text{e}}^{-t/{{\tau}_{x}}}},\qquad {{h}_{0}}<1,\end{eqnarray} \tag{ 139 }$$

where the inverse decay time is given by

$$\begin{eqnarray}&&\tau _{x}^{-1}={\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\varepsilon _{h}^{\prime}(k)\ln \left|\frac{1-{{K}^{2}}(k)}{1+{{K}^{2}}(k)}\right|.\end{eqnarray} \tag{ 140 }$$

Here the function K ²(k) has been previously defined in (132). The prefactor C ^x(t) was calculated in [26]

$$\begin{eqnarray}{{C}^{x}}(t)=\left\{\begin{array}{*{35}{l}} \sqrt{\frac{1-h{{h}_{0}}+\sqrt{\left(1-{{h}^{2}}\right)\left(1-h_{0}^{2}\right)}}{2\sqrt{1-h{{h}_{0}}}{{\left(1-h_{0}^{2}\right)}^{\frac{1}{4}}}}} & h<1 \\ {{\left[\frac{h\sqrt{1-h_{0}^{2}}}{h+{{h}_{0}}}\right]}^{\frac{1}{4}}}{{\left[1+\cos \left(2{{\varepsilon}_{h}}\left({{k}_{0}}\right)t+\alpha \right)+\ldots \right]}^{\frac{1}{2}}} & h>1. \end{array}\right. \end{eqnarray} \tag{ 141 }$$

The oscillatory behaviour for h  >  1 can be related to the presence of a gapless mode with momentum k₀ in the GGE Hamiltonian (130) (${{K}^{2}}\left({{k}_{0}}\right)=1$). It was noted in [42] that the period of the oscillations in (141) coincides with cusps in the time evolution of the logarithm of the return probability per unit length $f(t)=-\lim\nolimits_{L\to \infty}\frac{1}{L}\log |\langle \Psi(0)|{{\text{e}}^{-\text{i}Ht}}|\Psi(0)\rangle {{|}^{2}}$ after quenches between the phases. The result (139) has been derived in the late time limit $Jt\gg 1$. However, it gives an excellent account of the full answer except at very short times. This can be shown by comparing (139) to a numerical solution based on free fermion methods [26]. The latter works directly in the thermodynamic limit and does not suffer from finite-size effects. Figure 7 shows such comparisons for two different quenches within the ferromagnetic phase, and for two quenches from the ferromagnetic to the paramagnetic phase. The agreement is visibly excellent even at moderate times.

Zoom In Zoom Out Reset image size

At first sight the exponential decay (139) of the order parameter for quenches within the ferromagnetic phase may look surprising. Even a very small quench will lead to the eventual disappearance of the order parameter. A simple way of understanding this is to note that the ferromagnetic order persists only at zero temperature T  =  0, and melts for any T  >  0. By means of our quantum quench we deposit a finite energy density into the system, which is very similar to imposing a finite temperature. This consideration provides an intuitive explanation for why even small quenches wipe out the long range order present in the initial state. We note that this behaviour is specific to one dimensional systems, where discrete symmetries can be spontaneously broken only at T  =  0. In higher dimensional systems we expect order parameters to be generally stable when subjected to sufficiently small quantum quenches.

Finally we note that the expectation value $\langle \sigma _{j}^{y}(t)\rangle$ can be obtained in a simple way by considering the Heisenberg equations of motion for $\sigma _{j}^{x}$, and is given by

$$\begin{eqnarray}&&\langle \Psi(t)|\sigma _{j}^{y}|\Psi(t)\rangle =\frac{\partial}{\partial t}\frac{\langle \Psi(t)|\sigma _{j}^{x}|\Psi(t)\rangle}{2Jh}.\end{eqnarray} \tag{ 142 }$$

Transverse spin operator

The spin operator $\sigma _{j}^{z}=\text{i}{{a}_{2j}}{{a}_{2j-1}}$ is even and its expectation value can be straightforwardly calculated with free-fermion techniques [45, 47]. It decays like a t^−3/2 power law in time towards its stationary value ${{\langle \cdot \rangle}_{\text{GGE}}}$

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\langle \Psi(t)|\sigma _{l}^{z}|\Psi(t)\rangle {{|}_{Jt\gg 1}}={{\langle \sigma _{l}^{z}\rangle}_{\text{GGE}}}\\ \,\,\,\,\,\,+\frac{{{h}_{0}}-h}{4\sqrt{\pi}{{(2hJt)}^{\frac{3}{2}}}}\left[\frac{\sin (4J|1-h|t+\pi /4)}{\sqrt{|1-h|}|1-{{h}_{0}}|}-\frac{\sin (4J(1+h)t-\pi /4)}{\sqrt{1+h}(1+{{h}_{0}})}\right]+O({{(Jt)}^{-\frac{5}{2}}}).\end{array}\end{eqnarray} \tag{ 143 }$$

The relaxation to the stationary value is only algebraic, in agreement with (137).

6.4.2. Spin–spin correlators in the 'space-time scaling limit'

A particularly useful way of describing the time dependence of two-point functions after a quantum quench is by considering an asymptotic expansion around the so-called space-time scaling limit [25, 26]. The latter refers to the behaviour along a particular ray in space-time

$$\begin{eqnarray}&&t,\ell \to \infty,\quad \frac{{{v}_{\text{max}}}t}{\ell}=\kappa =\text{fixed}.\end{eqnarray} \tag{ 144 }$$

Here ${{v}_{\text{max}}}=\text{ma}{{\text{x}}_{k}}\frac{\text{d}{{\varepsilon}_{h}}(k)}{\text{d}k}$ is the maximal group velocity of elementary excitations of the post-quench Hamiltonian.

Transverse spin–spin correlator

In the space-time scaling limit, the asymptotic behavior of ${\rho}_{c}^{zz}\left(\ell,t\right)$ can be evaluated by means of Wick's theorem, followed by a stationary phase approximation. The leading behaviour is a t⁻¹ power-law decay, and the subleading corrections are power laws as well [26]

$$\begin{eqnarray}&&{\rho}_{c}^{zz}\left(\ell =\frac{{{v}_{\text{max}}}t}{\kappa},t\right)\sim \frac{{{D}^{z}}(t)}{{{\kappa}^{2}}t}+o\left({{t}^{-1}}\right).\end{eqnarray} \tag{ 145 }$$

Here D^z(t) is the sum of a constant contribution and oscillatory terms with constant amplitudes.

Longitudinal spin–spin correlator

In the space-time scaling limit the order parameter two-point function ${{\rho}^{xx}}\left(\ell,t\right)$ takes the form

$$\begin{eqnarray}&&{{\rho}^{xx}}(\ell,t)\simeq {{\mathcal{C}}^{x}}(\ell,t)\exp \left[{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\ln \left|\frac{1-{{K}^{2}}(k)}{1+{{K}^{2}}(k)}\right|\min (2\varepsilon _{h}^{\prime}(k)t,\ell )\right].\end{eqnarray} \tag{ 146 }$$

The function ${{\mathcal{C}}^{x}}\left(\ell,t\right)$ has been determined in [26].

• 1.  
  For quenches within the ferromagnetic phase, h₀, h  <  1, ${{\mathcal{C}}^{x}}\left(\ell,t\right)$ equals the constant denoted by $\mathcal{C}_{\text{FF}}^{x}$ in (C.3). For times smaller than the Fermi time

  $$\begin{eqnarray}&&{{t}_{\text{F}}}=\frac{\ell}{2{{v}_{\text{max}}}},\end{eqnarray} \tag{ 147 }$$

  equation (146) equals the square of the one-point function (139). Thus, in the space-time scaling limit, connected correlations vanish identically for times $t<{{t}_{\text{F}}}$ and begin to form only after the Fermi time. We stress that this does not imply that the connected correlations are exactly zero for $t<{{t}_{\text{F}}}$: in any model, both on the lattice or in the continuum there are exponentially suppressed terms (in $\ell$), which however vanish in the scaling limit.
• 2.  
  For quenches from the ferromagnetic phase to the paramagnetic phase the prefactor is given by

  $$\begin{eqnarray}&&{{\mathcal{C}}^{x}}\left(\ell,t\right)=\mathcal{C}_{\text{FP}}^{x}\left[1+{{\theta}_{\text{H}}}\left({{t}_{\text{F}}}-t\right)\left(\cos \left(2{{\varepsilon}_{h}}\left({{k}_{0}}\right)t+\alpha \right)+\ldots \right)\right],\end{eqnarray} \tag{ 148 }$$

  where $\mathcal{C}_{\text{FP}}^{x}$ is the constant defined in (C.4), while k₀ and α are the constants appearing in the one-point function (141). For $t<{{t}_{\text{F}}}$, (146) is simply the square of the corresponding one-point function, which ensures that connected correlations vanish for $t<{{t}_{\text{F}}}$ in the space-time scaling regime. We note that the expression for $t<{{t}_{\text{F}}}$ is a conjecture [26].
• 3.  
  For quenches within the paramagnetic phase one has

  $$\begin{eqnarray}&&{{\mathcal{C}}^{x}}\left(\ell,t\right)\simeq \mathcal{C}_{\text{PP}}^{x}\left(\ell \right)+{{\left({{h}^{2}}-1\right)}^{\frac{1}{4}}}\sqrt{4{{J}^{2}}h}{\int}_{-\pi}^{\pi}\frac{\text{d}k}{\pi}\frac{K(k)}{{{\varepsilon}_{k}}}\sin \left(2t{{\varepsilon}_{k}}-k\ell \right)+\ldots,\end{eqnarray} \tag{ 149 }$$

  where $\mathcal{C}_{\text{PP}}^{x}\left(\ell \,\right)$ is the function defined in (C.6). Equation (149) constitutes the leading order in a low-density expansion computed within the form-factor formalism. The exact expression for a generic (not small) quench is not known.
• 4.  
  For quenches from the paramagnetic to the ferromagnetic phase, for $t>{{t}_{\text{F}}}$, ${{\mathcal{C}}^{x}}\left(\ell,t\right)$ is independent of time and is given by $\mathcal{C}_{\text{PF}}^{x}\left(\ell \right)$ of (C.5). For $t<{{t}_{\text{F}}}$ the correlator is exponentially small and, to the best of our knowledge, there are no analytic predictions for its behaviour.

6.4.3. Long time asymptotics of connected spin–spin correlators at a fixed separation ℓ

As we have seen above, the quench dynamics of one and two point functions of quantum spins is rich and interesting. However, these are nonetheless very special observables. Ideally one would like to have access to the full reduced density matrix for a given subsystem size, as its matrix elements encode the evolution of all correlation functions, see (21). In practice this is only possible in very simple non-interacting examples, or for very small subsystem sizes. An example are quenches in the disordered phase of the TFIC [28, 176]. Here the reduced density matrix on the interval $[1,\ell ]$ is given by

$$\begin{eqnarray}&&{{\rho}_{\ell}}(t)=\frac{1}{Z}\exp \left(\frac{1}{4}\underset{\ell,n}{\sum}{{a}_{l}}{{W}_{lm}}{{a}_{m}}\right),\end{eqnarray} \tag{ 154 }$$

where a_2n and a_2n−1 are the Majorana fermion operators (106) and the factor Z ensures $\text{Tr}\left({{\rho}_{\ell}}(t)\right)=1$. The matrix W is related to the two-point function of Majorana fermions by [123, 124, 204]

$$\begin{eqnarray}&&\tanh \frac{W}{2}=\Gamma,\qquad {{\Gamma}_{jk}}=\langle \Psi(t)|{{a}_{k}}{{a}_{j}}|\Psi(t)\rangle -{{\delta}_{j,k}}=-{{\Gamma}_{kj}},\end{eqnarray} \tag{ 155 }$$

where the time evolved initial state $|\Psi(t)\rangle$ is given by (117). The matrix elements of the correlation matrix are simple (single) integrals and can be found in [26].

As we have argued above, at late times after global quantum quenches isolated quantum systems relax locally towards some steady states ${{\rho}^{\text{SS}}}$. How quickly this relaxation occurs can be efficiently measured by considering the distance of the time evolving reduced density matrix ${{\rho}_{B}}(t)$ from its steady state value ${\rho}_{B}^{\text{SS}}$, where B is a subsystem of a given size. This diagnostic can be implemented quite generally numerically as long as the subsystem size is small [205]. For models that can be mapped to non-interacting theories, it is possible to go considerably further. An example is the TFIC, which was considered in [28]. The first step is to introduce a measure of distance on the space of RDMs on a given subsystem. A convenient choice was introduced in [28]

$$\begin{eqnarray}&&\mathcal{D}({{\rho}_{1}},{{\rho}_{2}})=\sqrt{\frac{\text{Tr}[{{({{\rho}_{1}}-{{\rho}_{2}})}^{2}}]}{\text{Tr}[{{({{\rho}_{1}})}^{2}}]+\text{Tr}[{{({{\rho}_{1}})}^{2}}]}}.\end{eqnarray} \tag{ 156 }$$

Figure 8 shows results for the distance between the RMDs for the time-evolving and stationary states for quantum quenches within the disordered and ordered phases in the TFIC. Subsystems consisting of $\ell$ neighbouring sites are considered, where $\ell$ ranges from 10 to 150. In both cases the distance is seen to eventually decay as a power law in time

$$\begin{eqnarray}&&\mathcal{D}\left({{\rho}_{\ell}}(t),{\rho}_{\ell}^{\text{GGE}}\right)\propto {{\ell}^{2}}{{(Jt)}^{-\frac{3}{2}}}.\end{eqnarray} \tag{ 157 }$$

Zoom In Zoom Out Reset image size

The power-law decay of the distance in time can be interpreted by relating it to expectation values of local operators in the subsystem. Using (21) one can show that

$$\begin{eqnarray}&&\mathcal{D}({{\rho}_{1}},{{\rho}_{2}})={{(\overline{R{{(\mathcal{O})}^{2}}})}^{\frac{1}{2}}},\qquad R(\mathcal{O})=\frac{|\text{Tr}[({{\rho}_{1}}-{{\rho}_{2}})\mathcal{O}]|}{\sqrt{\text{Tr}{{[{{\rho}_{1}}\mathcal{O}]}^{2}}+\text{Tr}{{[{{\rho}_{2}}\mathcal{O}]}^{2}}}},\end{eqnarray} \tag{ 158 }$$

where the bar denotes an average on the space of operators acting on the spins in the subsystem, taken with respect to the probability distribution

$$\begin{eqnarray}&&P(\mathcal{O})=\frac{\text{Tr}{{[{{\rho}_{1}}\mathcal{O}]}^{2}}+\text{Tr}{{[{{\rho}_{2}}\mathcal{O}]}^{2}}}{\sum\nolimits_{{{\mathcal{O}}^{\prime}}}\text{Tr}{{[{{\rho}_{1}}{{\mathcal{O}}^{\prime}}]}^{2}}+\text{Tr}{{[{{\rho}_{2}}{{\mathcal{O}}^{\prime}}]}^{2}}}.\end{eqnarray} \tag{ 159 }$$

Here the sum is over all operators $\sigma _{1}^{{{\alpha}_{1}}}\sigma _{2}^{{{\alpha}_{2}}}\ldots \sigma _{\ell}^{{{\alpha}_{\ell}}}$, where ${{\alpha}_{j}}=0,x,y,z$ and where we have assumed for simplicity that the subsystem is the interval $[1,\ell ]$. The contribution from a given operator $\mathcal{O}$ to the distance is weighted by the square of its expectation value. This shows that (158) measures a mean relative difference between the expectation values of local operators in the two states.

As $\overline{R(\mathcal{O})}\leqslant {{(\overline{R{{(\mathcal{O})}^{2}}})}^{1/2}}=\mathcal{D}({{\rho}_{\ell}}(t),{\rho}_{\ell}^{\text{GGE}})$ we may use (158) to identify a time scale $t_{\text{rms}}^{\ast}$ associated with the relaxation of the 'typical' operator (with respect to the probability distribution (159))

$$\begin{eqnarray}&&Jt_{\text{rms}}^{\ast}\sim \ell _{B}^{\frac{4}{3}}.\end{eqnarray} \tag{ 160 }$$

The time scale $t_{\text{rms}}^{\ast}$ is very different from the ones governing the time evolution of the two-point functions of spin operators.

The von Neumann entropy (also known as entanglement entropy) of a density matrix $\rho$ is defined as

$$\begin{eqnarray}&&{{S}_{\text{vN}}}[\rho ]=-\text{Tr}[\rho \log \rho ].\end{eqnarray} \tag{ 161 }$$

If ρ is a reduced density matrix in a system that is in a pure state, S_vN measures the entanglement between the subsystem and its complement. Entanglement entropies have become a standard diagnostic for detecting and identifying quantum phase transitions. In the context of quantum quenches the time evolution of the von Neumann entropy and other entanglement measures provides very useful information about the spreading of correlations [16, 31, 48, 89, 109, 127, 182, 183, 206–209]. A key result obtained in [48] is that after quenches to conformal field theories the von Neumann entropy of a subsystem of length $\ell$ increases linearly in time until it eventually saturates (see the review by Calabrese and Cardy [88] in this volume)

$$\begin{eqnarray}{{S}_{\text{vN}}}\left[\rho \right]\sim \left\{\begin{array}{*{35}{l}} \frac{\pi cvt}{6\epsilon} & vt<\frac{\ell}{2} \\ \frac{\pi c\ell}{12\epsilon} & vt\gtrsim \frac{\ell}{2}. \end{array}\right. \end{eqnarray} \tag{ 162 }$$

Here c is the central charge of the CFT, is a constant with dimensions of length that depends on the initial state and v is the speed of light. The two behaviours in (162) connect smoothly over a region $|vt-\frac{\ell}{2}|\sim \epsilon$. A physical interpretation of the result (162) is provided by the Calabrese–Cardy quasi-particle picture [89] we already encountered in section 5. Its application to the time evolution of the entanglement entropy in integrable models initialized in squeezed states with finite correlation lengths proceeds as follows. The idea is that in a squeezed state correlations spread via the propagation of pairs of quasi-particles with equal but opposite momenta. At time t  =  0 the quantum quench generates such quasi-particles pairs throughout the system. Correlations between quasi-particles produced at a distance larger than the correlation length in the initial state can be neglected. The entanglement between a given region B and its complement is generated by quasi-particle pairs. The entanglement entropy is interpreted as a measure of the number of correlated pairs such that, at a given time, one quasi-particle is inside B and one outside, see figure 9. Entanglement is initially generated at the boundaries of B, and the entangled region spreads outwards in the form of two light cones. This picture suggests the following semiclassical expression for the von Neumann entropy

$$\begin{eqnarray}&&{{S}_{\text{vN}}}[{{\rho}_{B}}]\overset{\text{sc}}{{=}}{\int}^{}\text{d}k~f(k)\min ({{\ell}_{B}},2|v(k)|t).\end{eqnarray} \tag{ 163 }$$

Zoom In Zoom Out Reset image size

Here v (k) is the semiclassical group velocity $v(k)=\frac{\text{d}\varepsilon (k)}{\text{d}k}$, $\varepsilon (k)$ is the dispersion relation of the quasiparticles and f (k) is an unknown function that contains information on the initial state.

This behaviour has been observed in a variety of lattice and continuum models [31, 48, 89, 109, 127, 182, 183, 206, 209], including non-integrable cases [183], in which the quasi-particle picture does not apply by virtue of the finite quasi-particle life time. The semiclassical interpretation has been generalized to prethermalized regimes in models with weak integrability breaking [127].

The linear entanglement growth after quantum quenches has important ramifications. It is a crucial limiting factor for applying matrix-product state methods such as t-DRMG [210] and i iTEBD [211] algorithms to the computation of the dynamics after global quantum quenches. We note that translational invariance is essential for the quasi-particle picture to hold, and the time evolution of the entanglement entropy in e.g. disordered models [206] is very different.

Exact results for the evolution of the entanglement entropy after quenches in the transverse field Ising chain [48] are in accordance with the structure (163) suggested by the quasi-particle picture. In the limit $1\ll \ell,Jt$, the entanglement entropy of a block of $\ell$ neighbouring spins is

$$\begin{eqnarray}&&{{S}_{\text{vN}}}\left[{{\rho}_{B}}\right]\simeq {\int}_{\pi}^{\pi}\frac{\text{d}k}{2\pi}~w\left(\langle \Psi(0)|n(k)|\Psi(0)\rangle \right)~\min \left(\ell,2|{{\varepsilon}^{\prime}}(k)|t\right)+o\left(\ell \right),\end{eqnarray} \tag{ 164 }$$

where $w(x)=-x\log x-(1-x)\log (1-x)$ is the entropy per site (see appendix B), and $\langle \Psi(0)|n(k)|\Psi(0)\rangle$ is the (conserved) density of elementary excitations of the post-quench Hamiltonian H with momentum k at times t  >  0. It is given by

$$\begin{eqnarray}&&\langle \Psi(0)|n(k)|\Psi(0)\rangle =\frac{{{K}^{2}}(k)}{1+{{K}^{2}}(k)},\end{eqnarray} \tag{ 165 }$$

where K ²(k) was defined previously in (132). It follows from (164) that the entanglement entropy increases linearly until the Fermi time (147), and then slowly approaches its stationary value set by the GGE. The latter equals the entropy per site of the GGE for the entire system [109]

$$\begin{eqnarray}&&{\int}_{\pi}^{\pi}\frac{\text{d}k}{2\pi}~w(\langle \Psi(0)|n(k)|\Psi(0)\rangle )=\underset{L\to \infty}{{\lim}}\frac{1}{L}{{S}_{vN}}[{{\rho}^{\text{GGE}}}].\end{eqnarray} \tag{ 166 }$$

These observations persist for quantum quenches in the TFIC starting in excited states [31].

Interestingly, the stationary value of the entropy density in the diagonal ensemble differs from that in the GGE [30, 212]. This is not a problem, because the entropy per site is a global property of the system, while the equivalence between ensembles only holds for (finite) subsystems in the thermodynamic limit. A detailed explanation of the origin of this difference for the case of the TFIC was provided in [30].

We now turn to dynamical correlation functions

$$\begin{eqnarray}&&{{\rho}^{\alpha \alpha}}(\ell,t+{{t}_{1}},t+{{t}_{2}})=\langle \Psi(t)|\sigma _{j+\ell}^{\alpha}({{t}_{1}})\sigma _{j}^{\alpha}({{t}_{2}})|\Psi(t)\rangle.\end{eqnarray} \tag{ 167 }$$

The transverse correlator ($\alpha =z$ in (167)) can be calculated by elementary means [45], as $\sigma _{j}^{z}$ is quadratic in Jordan–Wigner fermions. The order parameter correlator ($\alpha =x$ in (167)) is much more difficult to evaluate. In [29] it was determined by means of a generalization of the form factor methods developed in [26], and by exploiting exact results in particular limits [27]. These methods are so far restricted to quenches within either the paramagnetic or the ferromagnetic phase, and lead to answers of the form [29]

$$\begin{eqnarray}&&{{\rho}^{xx}}\left(\ell,t+\tau,t+{{t}_{2}}\right)\simeq {{C}^{x}}\left(\ell,\tau,t\right)\exp \left[{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\log \left|\frac{1-{{K}^{2}}(k)}{1+{{K}^{2}}(k)}\right|\min \left\{\max \left[\varepsilon _{h}^{\prime}(k)\tau,\ell \right],\varepsilon _{h}^{\prime}(k)\left(2t+\tau \right)\right\}\right].\end{eqnarray} \tag{ 168 }$$

Here the functions ${{\epsilon}_{h}}(k)$ and K ²(k) are given in (112) and (132) respectively, and h₀ and h are the initial and final values of the transverse field. The function ${{C}^{(x)}}\left(\ell,\tau,t\right)$ depends on the phase in which the quench is performed:

• 1.  
  In the ferromagnetic phase, h₀, h  <  1, ${{C}^{x}}\left(\ell,\tau,t\right)$ equals the constant $\mathcal{C}_{\text{FF}}^{x}$ in (C.3).
• 2.  
  In the paramagnetic phase (h₀, h  >  1) ${{C}^{x}}\left(\ell,\tau,t\right)$ is given (to leading order in the form factor expansion) by [29]

  $$\begin{eqnarray}\begin{array}{*{35}{l}} {{C}^{x}}(\ell,\tau,t) & = & {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{J{{\text{e}}^{\text{i}\ell k}}}{{{\varepsilon}_{h}}(k)}[{{\text{e}}^{-\text{i}{{\varepsilon}_{h}}(k)t}}+2iK(k)\cos (2{{\varepsilon}_{h}}(k)(2t+\tau ))\text{sgn}(\ell -\varepsilon _{h}^{\prime}(k)\tau )] \\ {} & \times & \frac{h{{h}_{0}}-1+\sqrt{({{h}^{2}}-1)(h_{0}^{2}-1)}}{\sqrt{{{h}_{0}}-1/h}\sqrt[4]{h_{0}^{2}-1}},\qquad \ell <{{v}_{\text{max}}}(2t+\tau ). \end{array}\end{eqnarray} \tag{ 169 }$$

  In the complementary regime ${{v}_{\text{max}}}\left(2t+\tau \right)<\ell$ the correlator is exponentially small and this expression no longer applies.

The form factor result gives an excellent approximation to the exact answer (which can be computed numerically using free fermion techniques) for 'small' quenches [29]. These are defined as being characterized by having low densities of excitations in the initial state.

It is possible to obtain some of these results in an alternative way by generalizing the semiclassical approach of [39] to the non-equal time case, and then elevating it using exact limiting results derived in [25, 26]. This method fails to reproduce the result for quenches in the disordered phase outside the light cone ${{v}_{\text{max}}}\tau <\ell$, but is significantly simpler.

Integrability allows the construction of a basis of simultaneous eigenstates of the Hamiltonian and all its conservation laws. Unlike in the non-interacting case these states have a very complicated structure described by the Bethe Ansatz [20, 213, 214]. When we consider a quench between two integrable Hamiltonians $H\left({{h}_{0}}\right)\to H(h)$, we are thus faced with the problem of how to translate between the eigenbases of the two integrable theories. This is a difficult undertaking [215, 216], and no general formalism for achieving it is currently known. Progress has however been made in cases where the initial state has a simpler structure, in particular for (matrix) product states in either position [58, 63, 65, 66, 68, 69] or momentum/rapidity space [95, 97, 105, 107, 111, 113, 114]. There are two main methods for encoding the relevant information contained in the initial state.

• 1.  
  Let us denote the eigenstates of the post-quench Hamiltonian H(h) by $|n\rangle$. One way to implement the initial conditions is via the overlaps $\langle n|\Psi(0)\rangle$ [67]. If these are known, the initial state can be translated into the eigenbasis of the time evolution operator. This method is used in the quench action approach [176] (see the review by Caux [217] in this volume).
• 2.  
  Let us denote the local conservation laws of H(h) by {I ⁽ⁿ⁾}. If the set {I ⁽ⁿ⁾} is in some sense complete [112, 218, 219], then the initial conditions can be encoded in the constraints

  $$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\frac{\langle \Psi(0)|{{I}^{(n)}}|\Psi(0)\rangle}{L}=\underset{L\to \infty}{{\lim}}\frac{\text{Tr}\left({{\rho}^{\text{SS}}}{{I}^{(n)}}\right)}{L},\end{eqnarray} \tag{ 170 }$$

  where ${{\rho}^{\text{SS}}}$ is one of the density matrices (GGE, GMC, diagonal ensemble) that describes the local properties of the stationary state.

In the following we will discuss implementations of the second approach.

In free theories a convenient way for constructing the GGE is by exploiting the linear relation between the local conservation laws and the mode occupation operators, see (42). In interacting integrable models the situation is different. Like free theories they feature stable excitations. In the thermodynamic limit these can be described by creation and annihilation operators $Z_{a}^{\dagger}\left(\lambda \right)$, ${{Z}_{a}}\left(\lambda \right)$ (the index a labels different particle species) fulfilling the Faddeev–Zamolodchikov algebra [220, 221]

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{Z}_{a}}\left({{\lambda}_{1}}\right){{Z}_{b}}\left({{\lambda}_{2}}\right) & = & S_{ab}^{cd}\left({{\lambda}_{1}},{{\lambda}_{2}}\right){{Z}_{d}}\left({{\lambda}_{2}}\right){{Z}_{c}}\left({{\lambda}_{1}}\right), \\ {{Z}_{a}}\left({{\lambda}_{1}}\right)Z_{b}^{\dagger}\left({{\lambda}_{2}}\right) & = & 2\pi \delta \left({{\lambda}_{1}}-{{\lambda}_{2}}\right){{\delta}_{a,b}}+S_{bc}^{da}\left({{\lambda}_{2}},{{\lambda}_{1}}\right)Z_{d}^{\dagger}\left({{\lambda}_{2}}\right){{Z}_{c}}\left({{\lambda}_{1}}\right). \end{array}\end{eqnarray} \tag{ 171 }$$

Here λ parametrizes the momenta ${{p}_{a}}\left(\lambda \right)$ and $S_{ab}^{cd}\left({{\lambda}_{1}},{{\lambda}_{2}}\right)$ is the purely elastic two-particle S-matrix. The generalized mode occupation operators ${{N}_{a}}\left(\lambda \right)=Z_{a}^{\dagger}\left(\lambda \right){{Z}_{a}}\left(\lambda \right)$ then indeed provide a set of mutually commuting conserved charges. The problem is that, due to the interacting nature of the stable excitations, there is no simple way of defining such operators in the finite volume [95], which is the standard way of making the theory well defined (working directly in the thermodynamic limit requires the regularization of very complicated singularities [40, 95, 222], which appears impractical in general). The problem lies in the nature of the quantization conditions in the finite volume, which on a ring of length L read

$$\begin{eqnarray}&&{{\text{e}}^{\text{i}L{{p}_{a}}({\lambda}_{j}^{(a)})}}=-\underset{b,k}{\prod}{{S}_{ab}}({\lambda}_{j}^{(a)},{\lambda}_{k}^{(b)}).\end{eqnarray} \tag{ 172 }$$

The solutions to this complicated system of coupled equations are such that the possible values of ${\lambda}_{j}^{(a)}$, and hence ${{p}_{a}}({\lambda}_{j}^{(a)})$, depend in a very sensitive way on all the other particles present in a given excitation. This is fundamentally different from the non-interacting case, where the momenta are simply given by

$$\begin{eqnarray}&&{{p}_{a}}({\lambda}_{j}^{(a)})=\frac{2\pi}{L}\times ~\text{integer},\end{eqnarray} \tag{ 173 }$$

and are independent of the particle content of a given excitation. This makes it clear that defining finite volume analogues of ${{N}_{a}}(\lambda )$ is difficult.

Our paradigm for an interacting integrable model will be the spin-1/2 Heisenberg XXZ chain. Its Hamiltonian on a ring with L sites is

$$\begin{eqnarray}&&{{H}_{\text{XXZ}}}=\frac{J}{4}\underset{j=1}{\overset{L}{\sum}}\sigma _{j}^{x}\sigma _{j+1}^{x}+\sigma _{j}^{y}\sigma _{j+1}^{y}+\Delta [\sigma _{j}^{z}\sigma _{j+1}^{z}-1],\end{eqnarray} \tag{ 174 }$$

where we will assume for definiteness that

$$\begin{eqnarray}&&\Delta =\cosh (\eta )\geqslant 1.\end{eqnarray} \tag{ 175 }$$

From now on we set J  =  1.

7.3.1. Generalized microcanonical ensemble

For interacting integrable models the GMC is easier to work with than the GGE. It is based on working with macro-states obtained by taking the thermodynamic limit of eigenstates constructed from the Bethe Ansatz. This procedure is an essential ingredient of the Thermodyamic Bethe Ansatz and is reviewed in several monographs [213, 225]. A very brief summary is given in appendix B.2. The upshot is that macro-states in integrable models are characterized by an (infinite) set of densities $\{{{\rho}_{n,p}}(\lambda )|n=1,2,\ldots \}$, where n labels all distinct stable species of excitations in the model. A given macro-state corresponds to a set of micro-states $|\Phi\rangle$, called representative states in [176]. These are by construction simultaneous eigenstates of all local conservation laws. For the macro-state describing the stationary state after our quench, they satisfy the initial conditions

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\frac{\langle {{\Phi}^{\text{SS}}}|{{I}^{(n)}}|{{\Phi}^{\text{SS}}}\rangle}{L}=\underset{L\to \infty}{{\lim}}\frac{\langle \Psi(0)|{{I}^{(n)}}|\Psi(0)\rangle}{L}.\end{eqnarray} \tag{ 176 }$$

The GMC density matrix is then defined [176] in terms of a single such 'representative' micro-state $|{{\Phi}^{\text{SS}}}\rangle$

$$\begin{eqnarray}&&{{\rho}^{\text{GMC}}}=\,|{{\Phi}^{\text{SS}}}\rangle \langle {{\Phi}^{\text{SS}}}|.\end{eqnarray} \tag{ 177 }$$

Here we have assumed that the stationary state is given in terms of a single macro-state constructed from the Bethe Ansatz. In principle it is possible that the steady state has a more complicated structure and requires a description in terms of a sum of several density matrices of the form (177).

7.3.2. Transfer matrix and 'ultra-local' conservation laws

According to our general discussion, local observables should relax to an appropriate GGE after quenches to the XXZ chain. In order to construct this GGE, we need to know the required set of local conservation laws of (174). One family of conservation laws has been known for a long time and is most conveniently constructed by exploiting the relation of the Heisenberg Hamiltonian to the transfer matrix of the six-vertex model [20, 223, 224]. The fundamental building block of the six-vertex model is the L-operator

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{L}_{n}}\left(\lambda \right)=\frac{1}{\sinh \left(\eta +\text{i}\lambda \right)}\left[\sinh \left(\frac{\eta}{2}+\text{i}\lambda \right)\cosh \left(\frac{\eta}{2}\right)+\cosh \left(\frac{\eta}{2}+\text{i}\lambda \right)\sinh \left(\frac{\eta}{2}\right){{\tau}^{z}}\sigma _{n}^{z}\right. \\ \left. +\sinh \left(\eta \right)\left({{\tau}^{-}}\sigma _{n}^{+}+{{\tau}^{+}}\sigma _{n}^{-}\right)\right],\end{array}\end{eqnarray} \tag{ 178 }$$

which acts on the tensor product ${\mathbb{C}^{2}}\otimes {\mathbb{C}^{2}}$ of 'auxiliary' and 'quantum' spaces through the Pauli matrices ${{\tau}^{\alpha}}$ and ${{\sigma}^{\alpha}}$ respectively. Matrix elements in the auxiliary/quantum spaces are denotes by Roman/Greek letters respectively, e.g.

$$\begin{eqnarray}&&\left({{\tau}^{-}}\sigma _{n}^{+}\right)_{\alpha \beta}^{ab}=\tau _{ab}^{-}{{\left(\sigma _{n}^{+}\right)}_{\alpha \beta}}.\end{eqnarray} \tag{ 179 }$$

The vertex weights of the six-vertex model are obtained by taking matrix elements $L\left(\lambda \right)_{\alpha \beta}^{ab}$ and have a graphical representation as shown in figure 10(a). The row-to-row transfer matrix is obtained as shown in figure 10(b)

$$\begin{eqnarray}&&\left({{\tau}_{\frac{1}{2}}}(\lambda )\right)_{{{\alpha}_{1}}\ldots {{\alpha}_{L}}}^{{{\beta}_{1}}\ldots {{\beta}_{L}}}=({{L}_{1}}(\lambda ))_{{{\alpha}_{1}}{{\beta}_{1}}}^{{{a}_{1}}{{a}_{2}}}({{L}_{2}}(\lambda ))_{{{\alpha}_{2}}{{\beta}_{2}}}^{{{a}_{2}}{{a}_{3}}}\ldots ({{L}_{L}}(\lambda ))_{{{\alpha}_{L}}{{\beta}_{L}}}^{{{a}_{L}}{{a}_{1}}}.\end{eqnarray} \tag{ 180 }$$

Zoom In Zoom Out Reset image size

The partition function of the 6-vertex model on an $L\times M$ rectangular lattice with periodic boundary conditions is then

$$\begin{eqnarray}&&{{Z}_{6-\text{vertex}}}=\text{Tr}\left[{{\left({{\tau}_{\frac{1}{2}}}\left(\lambda \right)\right)}^{M}}\right],\end{eqnarray} \tag{ 181 }$$

where the trace is over the quantum space. As a consequence of the Yang–Baxter relation for the L-operators [20], the transfer matrices form a commuting family

$$\begin{eqnarray}&&[{{\tau}_{\frac{1}{2}}}\left(\lambda \right),{{\tau}_{\frac{1}{2}}}\left(\mu \right)]=0.\end{eqnarray} \tag{ 182 }$$

The Heisenberg Hamiltonian is related to the transfer matrix by taking a logarithmic derivative

$$\begin{eqnarray}&&{{H}_{\text{XXZ}}}=-\text{i}\frac{\sinh \eta}{2}\left. \frac{\partial}{\partial \lambda}\right|_{\lambda =0}\ln \left[{{\tau}_{\frac{1}{2}}}\left(\lambda \right)\right].\end{eqnarray} \tag{ 183 }$$

By virtue of the commutation relations (182) it is clear that a set of mutually commuting operators can be obtained by taking higher derivatives, i.e.

$$\begin{eqnarray}&&{{H}^{(\frac{1}{2},k)}}=\left. \text{i}{{\left(-\frac{\sinh \eta}{2}\frac{\partial}{\partial \lambda}\right)}^{k}}\right|_{\lambda =0}\ln \left[{{\tau}_{\frac{1}{2}}}\left(\lambda \right)\right].\end{eqnarray} \tag{ 184 }$$

Crucially, these conservation laws have the form

$$\begin{eqnarray}&&{{H}^{(\frac{1}{2},k)}}=\underset{j}{\sum}H_{j,j+1,\ldots,j+k}^{(\frac{1}{2},k)},\end{eqnarray} \tag{ 185 }$$

where the densities $H_{j,j+1,\ldots,j+k}^{(\frac{1}{2},k)}$ act non-trivially only on the k  +  1 consecutive sites $j,j+1,\ldots,j+k$. These conservation laws are sometimes referred to as ultra-local. They have been studied extensively in the literature [20, 226–228]. We note that the above construction is not restricted to Heisenberg models, but works much more generally [20].

7.3.3. 'Ultra-local' GGE

According to our general discussion, the GGE describing the steady state after a quench to the Heisenberg model should contain all of the conservation laws (184). An important question is whether these conservation laws are also sufficient. This was investigated in [58, 59, 63]. The basic idea is as follows. One considers time evolution induced by the Hamiltonian (174) starting from an initial state $|\Psi(0)\rangle$. The quantities of interest are the matrix elements of the reduced density matrix on a short interval in the steady state, see (21)

$$\begin{eqnarray}&&{{g}_{{{\alpha}_{1}},\ldots,{{\alpha}_{n}}}}=\underset{t\to \infty}{{\lim}}\underset{L\to \infty}{{\lim}}\langle \Psi(t)|\sigma _{1}^{{{\alpha}_{1}}}\sigma _{2}^{{{\alpha}_{2}}}\ldots \sigma _{n}^{{{\alpha}_{n}}}|\Psi(t)\rangle.\end{eqnarray} \tag{ 186 }$$

The question is whether these expectation values can be obtained, to a given accuracy, from a GGE density matrix of the form

$$\begin{eqnarray}&&{\rho}_{\text{ulGGE}}^{\left(\,y\right)}=\frac{1}{{{Z}_{\text{ulGGE}}}}\exp \left(-\underset{k=1}{\overset{y}{\sum}}{\lambda}_{k}^{\left(\,y\right)}{{H}^{(k,\frac{1}{2})}}\right),\end{eqnarray} \tag{ 187 }$$

which takes into account the first y conservation laws in the series (184). The Lagrange multipliers ${\lambda}_{k}^{\left(\,y\right)}$ are in principle fixed by the requirements

$$\begin{eqnarray}&&\langle \Psi(0)|{{H}^{\left(k,\frac{1}{2}\right)}}|\Psi(0)\rangle =\text{Tr}\left({\rho}_{\text{ulGGE}}^{\left(\,y\right)}\,{{H}^{\left(k,\frac{1}{2}\right)}}\right),\quad k=1,\ldots,y.\end{eqnarray} \tag{ 188 }$$

In practice it is very difficult to determine the Lagrange multipliers from these conditions, even for very simple initial states $|\Psi(0)\rangle$. A method to circumvent this problem was developed in [58]. The idea is to define a generating function for the initial values (188)

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Omega _{\Psi(0)}^{\left(\frac{1}{2}\right)}\left(\lambda \right) & = & \underset{L\to \infty}{{\lim}}\frac{\text{i}}{L}\langle \Psi(0)|\tau _{\frac{1}{2}}^{\prime}\left(\lambda \right)\tau _{\frac{1}{2}}^{-1}\left(\lambda \right)|\Psi(0)\rangle \\ {} & = & \underset{L\to \infty}{{\lim}}\frac{1}{L}\underset{k=1}{\sum}{{\left(-\frac{2}{\sinh \eta}\right)}^{k}}\frac{{{\lambda}^{k-1}}}{(k-1)\text{!}}\frac{\langle \Psi(0)|{{H}^{\left(k,\frac{1}{2}\right)}}|\Psi(0)\rangle}{L}. \end{array}\end{eqnarray} \tag{ 189 }$$

For large L (and real λ) the inverse of the transfer matrix becomes

$$\begin{eqnarray}&&\tau _{\frac{1}{2}}^{-1}\left(\lambda \right)\simeq {{\left({{\tau}_{\frac{1}{2}}}\left(\lambda \right)\right)}^{\dagger}}={{\left[\frac{\sinh \left(-\text{i}\lambda \right)}{\sinh \left(\eta -\text{i}\lambda \right)}\right]}^{L}}{{\tau}_{\frac{1}{2}}}\left(\lambda +\text{i}\eta \right).\end{eqnarray} \tag{ 190 }$$

Using the expression (180) of the transfer matrix as a product of L-operators, the generating function can thus be expressed in the form

$$\begin{eqnarray}&&\Omega _{\Psi(0)}^{(\frac{1}{2})}\left(\lambda \right)=\left. \underset{L\to \infty}{{\lim}}\frac{\text{i}}{L}\frac{\partial}{\partial x}\right|_{x=\lambda}\text{Sp}\langle \Psi(0)|{{V}_{L}}\left(x,\lambda \right)\ldots {{V}_{1}}\left(x,\lambda \right)|\Psi(0)\rangle,\end{eqnarray} \tag{ 191 }$$

where ${{V}_{n}}\left(x,\lambda \right)$ are $4\times 4$ matrices with entries $\left({{V}_{n}}\left(x,\lambda \right)\right)_{cd}^{ab}$ that are operators acting on the 2-dimensional quantum space at site n, and $\text{Sp}$ denotes the usual trace for $4\times 4$ matrices. The explicit expression is

$$\begin{eqnarray}&&{{[({{V}_{n}}(x,\lambda ))_{cd}^{ab}]}_{{{\alpha}_{n}}{{\beta}_{n}}}}=\frac{\sinh (-\text{i}\lambda )}{\sinh (\eta -\text{i}\lambda )}\underset{{{\gamma}_{n}}}{\sum}({{L}_{n}}(x))_{{{\alpha}_{n}}{{\gamma}_{n}}}^{ab}({{L}_{n}}(\lambda +\text{i}\eta ))_{{{\gamma}_{n}}{{\beta}_{n}}}^{cd}.\end{eqnarray} \tag{ 192 }$$

The advantage of representation (191) is that it can be efficiently evaluated for initial states $|\Psi(0)\rangle$ of matrix product form [58, 63]. To understand the principle behind this let us consider a translationally invariant product state

$$\begin{eqnarray}&&|\Psi(0)\rangle =\otimes _{j=1}^{L}|\psi {{\rangle}_{j}}.\end{eqnarray} \tag{ 193 }$$

In this case the generating function is obtained from the eigenvalues of the $4\times 4$ matrix $U(x,\lambda )=\,{{}_{1}}{{\langle \psi |{{V}_{1}}(x,\lambda )|\psi \rangle}_{1}}$ as

$$\begin{eqnarray}&&\Omega _{\Psi(0)}^{(\frac{1}{2})}(\lambda )=\left. \underset{L\to \infty}{{\lim}}\frac{\text{i}}{L}\frac{\partial}{\partial x}\right|_{x=\lambda}\text{Sp}(U{{(x,\lambda )}^{L}}).\end{eqnarray} \tag{ 194 }$$

An efficient algorithm for calculating $\Omega _{\Psi(0)}^{(\frac{1}{2})}(\lambda )$ for matrix-product states was given in [63]. Having encoded the 'initial data' of our quantum quench in the generating function $\Omega _{\Psi(0)}^{(\frac{1}{2})}(\lambda )$, we now move on to the calculation of expectation values of local operators in the state given by (187). A very useful observation is that ${\rho}_{\text{ulGGE}}^{(\,y)}$ can be viewed as a Gibbs ensemble for the 'Hamiltonian'

$$\begin{eqnarray}&&{{H}_{\text{eff}}}=\underset{k=1}{\overset{y}{\sum}}{\lambda}_{k}^{(\,y)}{{H}^{(k,\frac{1}{2})}}~\end{eqnarray} \tag{ 195 }$$

at an effective inverse temperature $\beta =1$. As all ${{H}^{(k,\frac{1}{2})}}$ commute and are obtained by taking logarithmic derivatives of the transfer matrix, finite temperature properties of (195) can be studied by standard methods [226, 229]. The Quantum Transfer Matrix approach (QTM) [230] is particularly useful in this regard, as it provides an efficient way to obtain explicit results for thermal averages of local operators [231, 232] (see also [233]). [58, 59, 63] employed the QTM approach to the calculation of steady state properties of the density matrix (187) for quenches from simple initial states. By employing the generating function (189), it is possible to take into account all ultra-local conservation laws, and arrive at explicit results for local observables without having to determine the Lagrange multipliers ${\lambda}_{k}^{(\infty )}$ [58]. In the QTM approach the stationary state is described in terms of the solution of a system of coupled, nonlinear integral equations. Remarkably, the information on the initial state enters this system only via the function $\Omega _{\Psi(0)}^{(\frac{1}{2})}(\lambda )$. Results for spin correlators obtained in this way were compared to t-DMRG computations for quenches from a variety of initial states in [63], and found to be compatible within the limitations of the numerical analysis.

The subsequent application of the quench action approach [176] (reviewed by Caux [217] in this volume) to the same problem revealed that the ultra-local GGE in fact does not correctly describe the steady state for quenches to the spin-1/2 Heisenberg chain [65, 66, 68–70], although it does provide a very good approximation for e.g. quenches from the Néel state. This suggested the existence of hitherto unknown conservation laws in the Heisenberg chain, which need to be taken into account in the construction of the GGE.

7.3.4. 'Quasi-local' GGE

The 'missing' conservation laws for the spin-1/2 Heisenberg XXZ chain were discovered in [64, 234, 237] (see the review by Ilievski, Medenjak, Prosen and Zadnik [235] in this volume). Their structure is quite different from that of the ultra-local conservation laws discussed above: their densities are not local in the sense that they act non-trivially only on a finite number of neighbouring sites, but quasi-local. Similar conservation laws had been identified earlier in relation to transport properties of the Heisenberg chain [236]. In order to define the concept of quasi-locality one introduces an inner product on the space of operators by

$$\begin{eqnarray}&&(A,B)={{\langle {{A}^{\dagger}}B\rangle}_{\infty}},\qquad {{\langle A\rangle}_{\infty}}=\frac{1}{{{2}^{L}}}\text{Tr}(A).\end{eqnarray} \tag{ 196 }$$

Definition 6. Quasi-local operators [234].

Let us consider an operator Q and expand it in terms of mutually orthogonal local operators q_j,r of range r

$$\begin{eqnarray}&&Q=\underset{j}{\sum}\underset{r}{\sum}{{q}_{j,r}}.\end{eqnarray} \tag{ 197 }$$

Q is called quasi-local if it fulfils the following three conditions:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \left(\text{QL}1\right): & \qquad & \underset{L\to \infty}{{\lim}}\frac{1}{L}\left(Q-{{\langle Q\rangle}_{\infty}},Q-{{\langle Q\rangle}_{\infty}}\right)=\text{const}; \\ \left(\text{QL}2\right): & \qquad & \underset{L\to \infty}{{\lim}}\left(Q,{{B}_{k}}\right)~~\text{exists}, \\ \left(\text{QL}3\right): & \qquad & \left({{q}_{j,r}},{{q}_{j,r}}\right)<C{{\text{e}}^{-r/\xi}}, \end{array}\end{eqnarray} \tag{ 198 }$$

where B_k is any operator that acts non-trivially only on a fixed number of k sites, and ξ and C are positive constants.

In the anisotropic Heisenberg chain quasi-local charges can be constructed as follows [64, 234]. It is well known that the six-vertex model transfer matrix (180) is part of a much larger family, built from the L-operators

$$\begin{eqnarray} \begin{array}{*{35}{l}} L_{n}^{(S)}\left(\lambda \right)=\frac{1}{\left. \sinh \left(\frac{\eta}{2}(1+s)+\text{i}\lambda \right)\right)} & \left[ \sinh \left(\frac{\eta}{2}+\text{i}\lambda \right)\cosh \left(\eta {{\mathcal{S}}^{z}}\right)+\cosh \left(\frac{\eta}{2}+\text{i}\lambda \right)\sinh \left(\eta {{\mathcal{S}}^{z}}\right)\sigma _{n}^{z}\right.\\ {} & \quad \left. +\sinh \left(\eta \right)\left({{\mathcal{S}}^{-}}\sigma _{n}^{+}+{{\mathcal{S}}^{+}}\sigma _{n}^{-}\right)\right), \end{array}\end{eqnarray} \tag{ 199 }$$

which again acts on the tensor product of auxiliary and quantum spaces, but now the auxiliary space is 2S  +  1 dimensional. Here S is an arbitrary half integer. The operators ${{\mathcal{S}}^{\alpha}}$ obey a q-deformed SU(2) algebra

$$\begin{eqnarray}&&\left[{{\mathcal{S}}^{+}},{{\mathcal{S}}^{-}}\right]={{\left[2{{\mathcal{S}}^{z}}\right]}_{q}},\qquad \left[{{\mathcal{S}}^{z}},{{\mathcal{S}}^{\pm}}\right]=\pm {{\mathcal{S}}^{\pm}},\end{eqnarray} \tag{ 200 }$$

where ${{\left[x\right]}_{q}}=\sinh \left(\eta x\right)/\sinh \left(\eta \right)$, and act on a q-deformed spin-S representation as

$$\begin{eqnarray}&&{{\mathcal{S}}^{z}}|k\rangle =k|k\rangle,\qquad {{\mathcal{S}}^{\pm}}|k\rangle =\sqrt{{{\left[S+1\pm k\right]}_{q}}{{\left[S\mp k\right]}_{q}}}|k\pm 1\rangle,\quad k=-S,\ldots,S.\end{eqnarray} \tag{ 201 }$$

A family of row-to-row transfer matrices is then obtained as

$$\begin{eqnarray}&&({{\tau}_{S}}(\lambda ))_{{{\alpha}_{1}}\ldots {{\alpha}_{L}}}^{{{\beta}_{1}}\ldots {{\beta}_{L}}}=(L_{1}^{(S)}(\lambda ))_{{{\alpha}_{1}}{{\beta}_{1}}}^{{{a}_{1}}{{a}_{2}}}(L_{2}^{(S)}(\lambda ))_{{{\alpha}_{2}}{{\beta}_{2}}}^{{{a}_{2}}{{a}_{3}}}\ldots (L_{L}^{(S)}(\lambda ))_{{{\alpha}_{L}}{{\beta}_{L}}}^{{{a}_{L}}{{a}_{1}}}.\end{eqnarray} \tag{ 202 }$$

All ${{\tau}_{S}}(\lambda )$ are operators on the same quantum space (a tensor product of L spin-1/2's), and as a consequence of the Yang–Baxter relation form a commuting family

$$\begin{eqnarray}&&[{{\tau}_{S}}(\lambda ),{{\tau}_{{{S}^{\prime}}}}(\mu )]=0.\end{eqnarray} \tag{ 203 }$$

By virtue of the commutation relations (182) it is clear that a set of mutually commuting operators can be obtained by taking higher derivatives, i.e.

$$\begin{eqnarray}&&{{H}^{(S,k)}}=\left. \text{i}{{\left({{C}_{S}}\frac{\partial}{\partial \lambda}\right)}^{k}}\right|_{\lambda =0}\ln \left[{{\tau}_{S}}\left(\lambda \right)\right]\,,\end{eqnarray} \tag{ 204 }$$

where C_S are some normalization constants that can be conveniently chosen. As a consequence of (203) we have

$$\begin{eqnarray}&&\left[{{H}^{(S,k)}},{{H}^{\left({{S}^{\prime}},{{k}^{\prime}}\right)}}\right]=0.\end{eqnarray} \tag{ 205 }$$

Apart from the special case S  =  1/2 these conservation laws are quasi-local. This means that (in the infinite volume) their general structure is

$$\begin{eqnarray}&&{{H}^{(S>\frac{1}{2},k)}}=\underset{j=-\infty}{\overset{\infty}{\sum}}\underset{k\geqslant 1}{\sum}\underset{{{\alpha}_{1}},\ldots,{{\alpha}_{k}}}{\sum}f_{{{\alpha}_{1}}{{\alpha}_{2}}\ldots {{\alpha}_{k}}}^{(k)}\sigma _{j}^{{{\alpha}_{1}}}\sigma _{j+1}^{{{\alpha}_{2}}}\ldots \sigma _{j+k-1}^{{{\alpha}_{k}}},\end{eqnarray} \tag{ 206 }$$

where ${{\alpha}_{j}}=0,x,y,z$, and the coefficient functions $f_{{{\alpha}_{1}}\ldots {{\alpha}_{k}}}^{(k)}$ decay sufficiently fast with k so that the conservation laws are extensive. As shown in [64], the initial data of the quantum quench can again be encoded in suitably chosen generating functions, which are generalizations of (189)

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\Omega _{\Psi(0)}^{(S)}\left(\lambda \right)=\underset{L\to \infty}{{\lim}}\frac{\text{i}}{L}\langle \Psi(0)|\tau _{S}^{\prime}\left(\lambda \right)\tau _{S}^{-1}\left(\lambda \right)|\Psi(0)\rangle \\ =\underset{L\to \infty}{{\lim}}\frac{\text{i}}{L}{{\left[\frac{\sinh \left(\frac{\eta}{2}(1-2S)-\text{i}\lambda \right)}{\sinh \left(\frac{\eta}{2}(1+2S)-\text{i}\lambda \right)}\right]}^{L}}\langle \Psi(0)|\tau _{S}^{\prime}\left(\lambda \right){{\tau}_{S}}\left(\lambda -\eta \right)|\Psi(0)\rangle.\end{array}\end{eqnarray} \tag{ 207 }$$

The generating functionals $\Omega _{\Psi(0)}^{(S)}\left(\lambda \right)$ can be evaluated for matrix product states by the same method discussed above (although the computational effort increases with the value of S).

The most convenient description of the stationary state turns out to be in terms of the generalized microcanonical ensemble discussed above. The steady state is characterized by the set $\{{\rho}_{n,p}^{\text{SS}}(\lambda )|n=1,\ldots \}$ of particle densities or the equivalent set of hole densities $\{{\rho}_{n,h}^{\text{SS}}(\lambda )|n=1,\ldots \}$. Ultimately this set must be determined by the initial conditions, which are encoded in (207). We now use that ${{X}_{S}}(\lambda )=\tau _{S}^{\prime}(\lambda )\tau _{S}^{-1}(\lambda )$ can be diagonalized by algebraic Bethe Ansatz [20, 238]. The eigenvalues of ${{X}_{S}}(\lambda )$ for M-particle states with $M\sim L$ are of the form

$$\begin{eqnarray}&&{{\nu}_{S}}(\lambda )=\underset{k=1}{\overset{M}{\sum}}\frac{2\sinh (2S\eta )}{\cos (2\lambda +2{{\lambda}_{k}})-\cosh (2S\eta )}+o(L),\end{eqnarray} \tag{ 208 }$$

where the ${{\lambda}_{k}}$ are solutions to the Bethe Ansatz equations [238]

$$\begin{eqnarray}&&{{\left(\frac{\sin \left({{\lambda}_{j}}+\text{i}\eta S\right)}{\sin \left({{\lambda}_{j}}-\text{i}\eta S\right)}\right)}^{L}}=\underset{k\ne j=1}{\overset{M}{\prod}}\frac{\sin \left({{\lambda}_{j}}-{{\lambda}_{k}}+\text{i}\eta \right)}{\sin \left({{\lambda}_{j}}-{{\lambda}_{k}}-\text{i}\eta \right)},\quad j=1,\ldots,M.\end{eqnarray} \tag{ 209 }$$

In the thermodynamic limit this can be simplified by following through the usual logic of the string hypothesis and the thermodynamic Bethe Ansatz [213, 214, 225]. Rather than with solutions to the Bethe equations (209) one then works with macro-states $|\boldsymbol{\rho}\rangle$, which are described by sets $\left\{{{\rho}_{n,p}}\left(\lambda \right)|n=1,\ldots \right\}$ of particle densities or the equivalent set of hole densities $\left\{{{\rho}_{n,h}}\left(\lambda \right)|n=1,\ldots \right\}$. The eigenvalue equation (208) then becomes [64]

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\frac{1}{L}\langle \boldsymbol{\rho}|{{X}_{S}}\left(\mu \right)|\boldsymbol{\rho}\rangle =\underset{k \in \mathbb{Z}}{\sum}\frac{{{\text{e}}^{-\text{i}2k\mu}}}{\cosh \left(k\eta \right)}\left({\int}_{-\pi /2}^{\pi /2}\text{d}\lambda \,{{\text{e}}^{2\text{i}k\lambda}}{{\rho}_{2S,h}}\left(\lambda \right)-{{\text{e}}^{-2S|k|\eta}}\right).\end{eqnarray} \tag{ 210 }$$

For $|\boldsymbol{\rho}\rangle =|\boldsymbol{\rho}^{{\text{SS}}}\rangle$ the right hand side of (210) must agree with the initial values after the quench (207). This is achieved by setting [64]

$$\begin{eqnarray}&&{\rho}_{2S,h}^{\text{SS}}\left(\lambda \right)={{a}_{2S}}\left(\lambda \right)+\frac{1}{2\pi}\left[\Omega _{\Psi(0)}^{(S)}\left(\lambda +\text{i}\frac{\eta}{2}\right)+\Omega _{\Psi(0)}^{(S)}\left(\lambda -\text{i}\frac{\eta}{2}\right)\right],\end{eqnarray} \tag{ 211 }$$

where ${{a}_{2S}}\left(\lambda \right)$ is a function independent of the initial state that is defined in appendix B.2. This shows that the initial data (207), which involves both ultra-local and quasi-local conservation laws, completely determines the macro-state that defines the generalized microcanonical ensemble. We note that the derivation did not invoke the maximum entropy principle. For the particular case of quenches from the Néel state (211) agrees with the one obtained by the quench action approach [65, 66, 68, 69].

The generalization of the approach discussed above for quenches to particular values of $\Delta$ with $-1<\Delta <1$ in the Heisenberg model (174) was achieved in [239].

We have defined our quench protocol such that it results in initial states $|\Psi(0)\rangle$ that have a cluster decomposition property (12). This requirement is often relaxed in both numerical and analytical investigations, and in some cases $|\Psi(0)\rangle$ is taken to be a Schrödinger cat state, see e.g. [65, 66, 68, 69]. An example is provided by quenches where the system is initialized in a classical Néel state $|\uparrow \downarrow \uparrow \downarrow \ldots \rangle$. This breaks translational invariance and it can be calculationally convenient to work instead with a translationally invariant cat state

$$\begin{eqnarray}&&\frac{1}{\sqrt{2}}\left(|\uparrow \downarrow \uparrow \downarrow \uparrow \ldots \rangle \,+\,|\downarrow \uparrow \downarrow \uparrow \downarrow \ldots \rangle \right).\end{eqnarray} \tag{ A.1 }$$

While for specific calculations such replacements can be useful, they significantly affect the steady state behaviour in general. This can be seen by considering a ${\mathbb{Z}_{2}}$ symmetric pre-quench Hamiltonian H₀ with a ground state that spontaneously breaks the ${\mathbb{Z}_{2}}$ symmetry. An example is provided by the transverse-field Ising chain (104) with h  <  1. In the thermodynamic limit H₀ has two ground states $|{{\Psi}_{\pm}}\rangle$, both of which have a cluster decomposition property. We now consider a general linear combination

$$\begin{eqnarray}&&|\Psi(0)\rangle =\cos \theta |{{\Psi}_{+}}\rangle +{{\text{e}}^{\text{i}\phi}}\sin \theta |{{\Psi}_{-}}\rangle.\end{eqnarray} \tag{ A.2 }$$

As $|{{\Psi}_{\pm}}\rangle$ are macroscopically distinct (as they lead to different order parameters), we conclude that expectation values of local operators $\mathcal{O}$ are given by

$$\begin{eqnarray}&&\langle \Psi(0)|\mathcal{O}|\Psi(0)\rangle ={{\cos}^{2}}\theta \langle {{\Psi}_{+}}|\mathcal{O}|{{\Psi}_{+}}\rangle +{{\sin}^{2}}\theta \langle {{\Psi}_{-}}|\mathcal{O}|{{\Psi}_{-}}\rangle.\end{eqnarray} \tag{ A.3 }$$

As the Hamiltonian of our system is short ranged, this decomposition persists at all finite times, i.e.

$$\begin{eqnarray}&&\langle \Psi(t)|\mathcal{O}|\Psi(t)\rangle ={{\cos}^{2}}\theta \langle {{\Psi}_{+}}(t)|\mathcal{O}|{{\Psi}_{+}}(t)\rangle +{{\sin}^{2}}\theta \langle {{\Psi}_{-}}(t)|\mathcal{O}|{{\Psi}_{-}}(t)\rangle \end{eqnarray} \tag{ A.4 }$$

To see this, we may use a result derived in [180] for the time evolution of local operators with short-ranged Hamiltonians: restricting $\mathcal{O}(t)={{\text{e}}^{\text{i}Ht}}\mathcal{O}{{\text{e}}^{-\text{i}Ht}}$ to a subsystem S of size $|S|$ gives an error that scales as ${{\text{e}}^{\left(2vt-|S|\right)/\zeta}}$, where ζ is a constant and v is the Lieb–Robinson velocity. This means that it is possible to approximate $\mathcal{O}(t)$ to a given accuracy by a local operator of a 'size' that scales as $2vt$. This in turn implies that $\langle {{\Psi}_{+}}|\mathcal{O}(t)|{{\Psi}_{-}}\rangle =0$, because $|{{\Psi}_{\pm}}\rangle$ are macroscopically distinct.

On the other hand, our system relaxes locally by construction if we initialize it in $|{{\Psi}_{\pm}}(0)\rangle$, i.e.

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}|{{\Psi}_{\pm}}(t)\rangle \langle {{\Psi}_{\pm}}(t)|\,={{\,}_{\text{loc}}}{\rho}_{\pm}^{\text{SS}}.\end{eqnarray} \tag{ A.5 }$$

Putting everything together we conclude that

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}|\Psi(t)\rangle \langle \Psi(t)|\,={{\,}_{\text{loc}}}{{\rho}^{\text{SS}}}={{\cos}^{2}}\theta {\rho}_{+}^{\text{SS}}+{{\sin}^{2}}\theta {\rho}_{-}^{\text{SS}}.\end{eqnarray} \tag{ A.6 }$$

The problem is that this form of the stationary state can be different from what one would expect on the basis of local relaxation to (generalized) Gibbs ensembles. To be specific, we consider the example of a quench to an anisotropic spin-1/2 Heisenberg chain

$$\begin{eqnarray}&&H=\frac{J}{4}\underset{\ell}{\sum}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}+\sigma _{\ell}^{z}\sigma _{\ell +1}^{z}+\Delta \sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+g\sigma _{\ell}^{x}\sigma _{\ell +2}^{x}.\end{eqnarray} \tag{ A.7 }$$

Imposing $\Delta,g\ne 0$ renders this model non-integrable, but is has one local conservation law

$$\begin{eqnarray}&&Q=\underset{\ell}{\sum}\sigma _{\ell}^{x},\quad \left[H,Q\right]=0.\end{eqnarray} \tag{ A.8 }$$

Our initial state is of the form (A.2), where $|{{\Psi}_{\pm}}(0)\rangle$ are the two ground states of the TFIC in the thermodynamic limit at h  <  1. We note that the energy density e of (A.7) is the same in both $|{{\Psi}_{+}}\rangle$ and $|{{\Psi}_{-}}\rangle$. As these states have a cluster decomposition property, it follows from our general discussion that the respective stationary states are locally equivalent to grand canonical ensembles

$$\begin{eqnarray}&&\underset{t\to \infty}{{\lim}}|{{\Psi}_{\pm}}(t)\rangle \langle {{\Psi}_{\pm}}(t)|\,={{\,}_{\text{loc}}}{\rho}_{\pm}^{\text{GC}}=\frac{{{\text{e}}^{-{{\beta}_{\pm}}H-{{\mu}_{\pm}}Q}}}{{{Z}_{\pm}}},\end{eqnarray} \tag{ A.9 }$$

where the values for ${{\beta}_{\pm}}$ and ${{\mu}_{\pm}}$ are obtained by fixing the energy density e and charge density ${{q}_{\pm}}$ to agree with their initial values at time t  =  0. Equation (A.6) tells us that the correct stationary state in this example is then

$$\begin{eqnarray}&&{{\rho}^{\text{SS}}}{{=}_{\text{loc}}}{{\cos}^{2}}\theta {\rho}_{+}^{\text{GC}}+{{\sin}^{2}}\theta {\rho}_{-}^{\text{GC}}.\end{eqnarray} \tag{ A.10 }$$

Crucially, while the density matrices ${\rho}_{\pm}^{\text{GC}}$ have the cluster decomposition property, ${{\rho}^{\text{SS}}}$ does not. This can be seen as follows. The local conservation law distinguishes between the two states $|{{\Psi}_{\pm}}\rangle$

$$\begin{eqnarray}&&{{q}_{+}}=\langle {{\Psi}_{+}}|\sigma _{\ell}^{x}|{{\Psi}_{+}}\rangle \ne \langle {{\Psi}_{-}}|\sigma _{\ell}^{x}|{{\Psi}_{-}}\rangle ={{q}_{-}}.\end{eqnarray} \tag{ A.11 }$$

Using the cluster decomposition property of ${\rho}_{\pm}^{\text{GC}}$ we have

$$\begin{eqnarray}&&\underset{|n-\ell |\to \infty}{{\lim}}\text{Tr}({\rho}_{\pm}^{\text{GC}}\sigma _{\ell}^{x}\sigma _{n}^{x})=q_{\pm}^{2},\end{eqnarray} \tag{ A.12 }$$

which in turn establishes that ${{\rho}^{\text{SS}}}$ does not have the cluster decomposition property

$$\begin{eqnarray}&&\underset{|n-\ell |\to \infty}{{\lim}}[\text{Tr}({{\rho}^{\text{SS}}}\sigma _{\ell}^{x}\sigma _{n}^{x})-\text{Tr}({{\rho}^{\text{SS}}}\sigma _{\ell}^{x})\text{Tr}({{\rho}^{\text{SS}}}\sigma _{n}^{x})]=\frac{{{({{q}_{+}}-{{q}_{-}})}^{2}}}{4}{{\sin}^{2}}(2\theta )\ne 0.\end{eqnarray} \tag{ A.13 }$$

On the other hand, if we were to apply our formalism of local relaxation blindly to our Hamiltonian (A.7), we would conclude that the stationary state is locally equivalent to a grand canonical ensemble, which is expected to have a cluster decomposition property [267].

Our discussion can be summarized as follows: If there exists at least one integral of motion that distinguishes $|{{\Psi}_{-}}\rangle$ from $|{{\Psi}_{+}}\rangle$, the stationary state associated with the time evolution of the cat state (A.2) does not possess the cluster decomposition property (12) and hence is not described by a standard generalized Gibbs ensemble.

Let us consider a post-quench Hamiltonian H with a set of local conservation laws I ⁽ⁿ⁾, see (31) and (32). In our basic definition of a quantum quench we initialize the system in a pure state $|\Psi(0)\rangle$. Then the cluster decomposition property implies that the probability distribution of energy and all local conservation laws approach delta-functions in the thermodynamic limit, e.g.

$$\begin{eqnarray}&&{{P}_{n}}(\epsilon )=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}[\rho (0)\delta ({{I}^{(n)}}-L\epsilon )]=\delta (\varepsilon -{{i}^{(n)}}),\quad {{i}^{(n)}}=\underset{L\to \infty}{{\lim}}\frac{1}{L}\text{Tr}[\rho (0){{I}^{(n)}}],\end{eqnarray} \tag{ A.14 }$$

where $\rho (0)=|\Psi(0)\rangle \langle \Psi(0)|$. As we have pointed out, it is sometimes desirable to consider initial density matrices $\rho (0)$ that are not pure states. When doing so one must ensure that (A.14) continues to hold. In cases where it does not it is clearly impossible for the system to locally relax to a GGE, because there the probability distributions of all conservation laws approach delta functions in the thermodynamic limit. Generalizations of GGE ideas to such cases have been explored in [266].

Let us consider a model of free fermions with Hamiltonian

$$\begin{eqnarray}&&H=\underset{k}{\sum}\epsilon (k)n(k),\end{eqnarray} \tag{ B.1 }$$

where $n(k)={{c}^{\dagger}}(k)c(k)$. Imposing periodic boundary conditions quantizes the allowed momenta

$$\begin{eqnarray}&&{{k}_{n}}=\frac{2\pi n}{L},\quad n=-\frac{L}{2}+1,\ldots,\frac{L}{2}.\end{eqnarray} \tag{ B.2 }$$

We now focus on a special class of Fock states $\prod\nolimits_{j=1}^{N}{{c}^{\dagger}}\left({{k}_{j}}\right)|0\rangle$, for which the particle densities

$$\begin{eqnarray}&&{{\rho}_{p}}\left({{k}_{j}}\right)=\frac{1}{L\left({{k}_{j+1}}-{{k}_{j}}\right)}\end{eqnarray} \tag{ B.3 }$$

approach smooth functions in the thermodynamic limit $N,L\to \infty$, n  =  N/L fixed. For such states, the number of particles in the interval $\left[{{k}_{n}},{{k}_{n}}+\Delta k\right]$ for large L is given by

$$\begin{eqnarray}&&{{\rho}_{p}}\left({{k}_{n}}\right)\Delta k.\end{eqnarray} \tag{ B.4 }$$

It is convenient to define a hole density by ${{\rho}_{h}}\left({{k}_{j}}\right)=\frac{1}{2\pi}-{{\rho}_{p}}\left({{k}_{j}}\right)$. In the thermodynamic limit many different choices of {k_j} lead to the same macro-state described by a given particle density. To enumerate them we note that the number of different states in the interval $\left[k,k+\Delta k\right]$, that give rise to a given density in the thermodynamic limit, is obtained by distributing ${{\rho}_{p}}(k)L\Delta k$ particles among $\left({{\rho}_{p}}(k)+{{\rho}_{h}}(k)\right)L\Delta k$ vacancies. This follows a binomial distribution. For large L the latter may be approximated by Stirling's formula, and in the thermodynamic limit one obtains the well-known expression for the entropy per site

$$\begin{eqnarray}&&s\left[{{\rho}_{p}}\right]={\int}_{-\pi}^{\pi}\text{d}k\left[\left({{\rho}_{p}}(k)+{{\rho}_{h}}(k)\right)\ln \left({{\rho}_{p}}(k)+{{\rho}_{h}}(k)\right)-{{\rho}_{p}}(k)\ln \left({{\rho}_{p}}(k)\right)-{{\rho}_{h}}(k)\ln \left({{\rho}_{h}}(k)\right)\right].\end{eqnarray} \tag{ B.5 }$$

Let us now investigate what kind of macro-states exist for a given energy density e. The most likely (maximum entropy) macro-state can be constructed using equilibrium statistical mechanics. Extremizing the free energy per site

$$\begin{eqnarray}&&f\left[{{\rho}_{p}}\right]={\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\epsilon (k){{\rho}_{p}}(k)-Ts\left[{{\rho}_{p}}\right]\end{eqnarray} \tag{ B.6 }$$

with respect to the particle density ${{\rho}_{p}}$ gives

$$\begin{eqnarray}&&{\rho}_{p}^{\text{eq}}(k)=\frac{1}{2\pi}\frac{1}{1+{{\text{e}}^{\epsilon (k)/T}}}.\end{eqnarray} \tag{ B.7 }$$

The 'temperature' T is related to the energy density e by

$$\begin{eqnarray}&&e={\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{\epsilon (k)}{1+{{\text{e}}^{\epsilon (k)/T}}}.\end{eqnarray} \tag{ B.8 }$$

The entropy per site of this equilibrium state is

$$\begin{eqnarray}&&s\left[{\rho}_{p}^{\text{eq}}\right]=\frac{\partial}{\partial T}T{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\ln \left[1+{{\text{e}}^{-\epsilon (k)/T}}\right].\end{eqnarray} \tag{ B.9 }$$

By construction the macro-state (B.7) is the typical state at energy density (B.8): if we randomly pick an energy eigenstate with energy density e for a very large system size L, the probability for this state to have particle density (B.7) is exponentially close (in L) to one. On the other hand, there are atypical finite entropy macro-states characterized by their respective particle densities ${{\rho}_{p}}(k)$. As a particular example we consider our tight-binding model (40) for $\mu =0$, which gives a dispersion $\epsilon (k)=-2J\cos (k)$. We fix the particle density to be 1/2 and the energy density to be $e=-0.405\,838J$, which corresponds to temperature T  =  J in (B.7), (B.8), i.e.

$$\begin{eqnarray}&&{\rho}_{p}^{\text{eq}}(k)=\frac{1}{2\pi}\frac{1}{1+{{\text{e}}^{-2\cos (k)}}}.\end{eqnarray} \tag{ B.10 }$$

The entropy of the equilibrium state is ${{s}_{\text{eq}}}=0.511\,571$. Let us now consider the family of macro-states described by the particle density

$$\begin{eqnarray}&&{\rho}_{p}^{\left(\lambda \right)}(k)=\frac{1}{2\pi}\frac{1}{1+{{\text{e}}^{-4\cos (k)-\lambda \cos (3k)}}}.\end{eqnarray} \tag{ B.11 }$$

Fixing $\lambda =2.430\,96...$ gives us the same particle and energy densities as for the equilibrium state, i.e. n  =  1/2 and $e=-0.405\,838J$. The entropy density $s=0.396\,781$ is of course lower than that of the (maximum entropy) equilibrium state. This means that if we randomly select an energy eigenstate with particle density n  =  1/2 and energy density $e=-0.405\,838J$ for a very large system size L, we are exponentially more likely by a factor ${{\text{e}}^{L\left({{s}_{\text{eq}}}-s\right)}}$ to end up with a state described by (B.10) than one described by (B.11). Unsurprisingly, expectation values of local operators are generally different in the two macro-states. As an example, let us consider

$$\begin{eqnarray}&&{{\mathcal{O}}_{j}}=c_{j}^{\dagger}{{c}_{j+3}}+c_{j+3}^{\dagger}{{c}_{j}}.\end{eqnarray} \tag{ B.12 }$$

Setting again $\lambda =2.430\,96...$ and taking the thermodynamic limit we have

$$\begin{eqnarray}&&\langle {\rho}_{p}^{\text{eq}}(k)|{{\mathcal{O}}_{j}}|{\rho}_{p}^{\text{eq}}(k)\rangle =-0.027\,1229\ne \langle {\rho}_{p}^{\left(\lambda \right)}(k)|{{\mathcal{O}}_{j}}|{\rho}_{p}^{\left(\lambda \right)}(k)\rangle =0.215\,148.\end{eqnarray} \tag{ B.13 }$$

Here the expectation values may be taken with regards to any micro-state that gives rise to the appropriate macro-state in the thermodynamic limit.

The situation in integrable models is analogous to what we just discussed for free fermions. The main difference arises from the more complicated structure and interacting nature of the elementary excitations in integrable models. For the sake of definiteness we consider the particular example of the spin-1/2 Heisenberg model (174) with $\Delta \geqslant 1$. Energy eigenstates $|{{\lambda}_{1}},\ldots,{{\lambda}_{M}}\rangle$ on a ring of length L are parametrized by M rapidity variables ${{\lambda}_{j}}$, which fulfil the quantization conditions [225]

$$\begin{eqnarray}&&{{\left(\frac{\sin ({{\lambda}_{j}}+\text{i}\frac{\eta}{2})}{\sin ({{\lambda}_{j}}-\text{i}\frac{\eta}{2})}\right)}^{L}}=\underset{k\ne j=1}{\overset{M}{\prod}}\frac{\sin ({{\lambda}_{j}}-{{\lambda}_{k}}+\text{i}\eta )}{\sin ({{\lambda}_{j}}-{{\lambda}_{k}}-\text{i}\eta )},\quad j=1,\ldots,M.\end{eqnarray} \tag{ B.14 }$$

Each ${{\lambda}_{j}}$ is associated with an elementary 'magnon' excitation over the ferromagnetic state with all spins up, and the total energy and momentum of the eigenstate are additive $E=\sum\nolimits_{j=1}^{M}e\left({{\lambda}_{j}}\right)$, $P=\sum\nolimits_{j=1}^{M}p\left({{\lambda}_{j}}\right)$. As a result of interactions, magnons form bound states. These correspond to 'string solutions' of the quantization conditions (B.14)

$$\begin{eqnarray}&&{\lambda}_{\alpha}^{n,j}={\lambda}_{\alpha}^{n}+\text{i}\frac{\eta}{2}\left(n+1-2j\right)+\text{i}\delta _{\alpha}^{n,j},\qquad j=1,\ldots,n,\end{eqnarray} \tag{ B.15 }$$

where $\delta _{\alpha}^{n,j}$ are deviations from 'ideal' strings that become negligible when we take the thermodynamic limit at finite densities of magnons and bound states, see 6.2.A of [268]. As a consequence of integrability these bound states are stable excitations. The generalization of particle density description of macro-states in free theories to interacting integrable models is then clear: macro-states are characterized by an (infinite) set of densities $\left\{{{\rho}_{n,p}}\left(\lambda \right)|n=1,2,\ldots \right\}$ for magnons (n  =  1) and bound states of all lengths $\left(n\geqslant 2\right)$. Just as in the case of free fermions, we can define corresponding hole densities ${{\rho}_{n,h}}\left(\lambda \right)$. The relation between particle and hole densities is fixed by the quantization conditions (B.14)

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\rho}_{n,t}}\left(\lambda \right)\equiv {{\rho}_{n,p}}\left(\lambda \right)+{{\rho}_{n,h}}\left(\lambda \right) & = & {{a}_{n}}\left(\lambda \right)-\underset{m=1}{\overset{\infty}{\sum}}{\int}_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\text{d}\mu ~{{A}_{nm}}\left(\lambda -\mu \right){{\rho}_{m,p}}\left(\mu \right), \end{array}\end{eqnarray} \tag{ B.16 }$$

where ${{A}_{nm}}\left(\lambda \right)=\left(1-{{\delta}_{n,m}}\right){{a}_{|n-m|}}\left(\lambda \right)+2{{a}_{|n-m|+2}}\left(\lambda \right)+\cdots +2{{a}_{n+m-2}}\left(\lambda \right)+{{a}_{n+m}}\left(\lambda \right)$ and ${{a}_{n}}\left(\lambda \right)=\frac{1}{2\pi}\frac{2\sinh \left(n\eta \right)}{\cosh \left(n\eta \right)-\cos \left(2\lambda \right)}$. Equation (B.16) allows one to express the hole densities in terms of the particle densities and vice versa, but in contrast to the non-interacting case the relationship is non-trivial. The energy and entropy per site of a macro-state are given by

$$\begin{eqnarray}\begin{array}{r l l} e\left[\left\{{{\rho}_{n,p}}\right\}\right] & = & -J\sinh \left(\eta \right)\underset{n=1}{\overset{\infty}{\sum}}{\int}_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\text{d}\lambda ~{{\rho}_{n,p}}\left(\lambda \right)~{{a}_{n}}\left(\lambda \right), \\ s\left[\left\{{{\rho}_{n,p}},{{\rho}_{n,h}}\right\}\right] & = & \underset{n=1}{\overset{\infty}{\sum}}{\int}_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\text{d}\lambda \left[\left({{\rho}_{n,t}}\left(\lambda \right)\right)\ln \left({{\rho}_{n,t}}\left(\lambda \right)\right)-{{\rho}_{n,p}}\left(\lambda \right)\ln \left({{\rho}_{n,p}}\left(\lambda \right)\right)-{{\rho}_{n,h}}\left(\lambda \right)\ln \left({{\rho}_{n,h}}\left(\lambda \right)\right)\right]. \end{array}\end{eqnarray} \tag{ B.17 }$$

The typical state at a given energy density is then obtained by extremizing the free energy per set $f\left[\left\{{{\rho}_{n,p}}\right\}\right]=e\left[\left\{{{\rho}_{n,p}}\right\}\right]-Ts\left[\left\{{{\rho}_{n,p}},{{\rho}_{n,h}}\right\}\right]$ with respect to particle and hole densities under the constraints (B.16). This results in the so-called Thermodynamic Bethe Ansatz equations [225]. By construction this state is thermal, i.e. corresponds to a standard Gibbs distribution.

Atypical finite entropy states can be constructed by specifying a set of particle densities $\left\{{{\rho}_{n,p}}\left(\lambda \right)\right\}$. The corresponding hole densities are then obtained by solving equations (B.16). The entropy per site of the resulting macro-state can be calculated from (B.17), and will by construction be smaller than that of the maximum entropy state.

Here the amplitude is of the form [26]

$$\begin{eqnarray}{{C}^{z}}\left(\ell \right)={{\bar{C}}^{z}}{{\ell}^{-{{\alpha}^{z}}}},\qquad {{\alpha}^{z}}=\left\{\begin{array}{*{35}{l}} 1 & \text{if}~~|\ln h|>|\ln {{h}_{0}}|, \\ 0 & \text{if}~~{{h}_{0}}=1/h, \\ 1/2 & \text{if}~~|\ln h|<|\ln {{h}_{0}}|, \end{array}\right. \end{eqnarray} \tag{ C.1 }$$

where the constant is known exactly

$$\begin{eqnarray}{{C}^{z}}=\left\{\begin{array}{*{35}{l}} \frac{|{{h}_{0}}-1/{{h}_{0}}|}{4\pi}\frac{{{h}_{0}}-h}{h{{h}_{0}}-1} & \text{if}~~|\ln h|>|\ln {{h}_{0}}|, \\ -\frac{{{\left(h-1/h\right)}^{2}}}{2\pi} & \text{if}~~{{h}_{0}}=1/h, \\ \frac{\left(h-1/h\right)\sqrt{|{{h}_{0}}-1/{{h}_{0}}|}\left({{h}_{0}}-h\right)}{8\sqrt{\pi}h}\sqrt{\frac{{{h}_{0}}-h}{{{h}_{0}}\left(h{{h}_{0}}-1\right)}}\frac{{{\text{e}}^{\text{sgn}\left(\ln h\right)|\ln {{h}_{0}}|/2}}}{\sinh \frac{|\ln h|+|\ln {{h}_{0}}|}{2}} & \text{if}~~|\ln h|<|\ln {{h}_{0}}|. \end{array}\right. \end{eqnarray} \tag{ C.2 }$$

Here the large-$\ell$ asymptotics of prefactor ${{C}^{x}}\left(\ell \right)$ are as follows.

• 1.  
  Quench within the ferromagnetic phase (h₀, h  <  1).

  $$\begin{eqnarray}&&{{C}^{x}}\left(\ell \right)\equiv \mathcal{C}_{\text{FF}}^{x}=\frac{1-h{{h}_{0}}+\sqrt{\left(1-{{h}^{2}}\right)\left(1-h_{0}^{2}\right)}}{2\sqrt{1-h{{h}_{0}}}\sqrt[4]{1-h_{0}^{2}}}.\end{eqnarray} \tag{ C.3 }$$
• 2.  
  Quench from the ferromagnetic to the paramagnetic phase (h₀  <  1  <  h).

  $$\begin{eqnarray}&&{{C}^{x}}\left(\ell \right)\equiv \mathcal{C}_{\text{FP}}^{x}=\sqrt{\frac{h\sqrt{1-h_{0}^{2}}}{h+{{h}_{0}}}}.\end{eqnarray} \tag{ C.4 }$$
• 3.  
  Quench from the paramagnetic to the ferromagnetic phase (h₀  >  1  >  h).

  $$\begin{eqnarray}&&{{C}^{x}}\left(\ell \right)\equiv \mathcal{C}_{\text{PF}}^{x}\left(\ell \right)=\sqrt{\frac{~{{h}_{0}}-h}{\sqrt{h_{0}^{2}-1}}}\cos \left(\ell \arctan \frac{\sqrt{\left(1-{{h}^{2}}\right)\left(h_{0}^{2}-1\right)}}{1+{{h}_{0}}h}\right).\end{eqnarray} \tag{ C.5 }$$
• 4.  
  Quench within the paramagnetic phase (1  <  h₀, h).

  $$\begin{eqnarray}{{C}^{x}}\left(\ell \right)\equiv \mathcal{C}_{\text{PP}}^{x}\left(\ell \right)=\left\{\begin{array}{*{35}{l}} -\frac{{{h}_{0}}\sqrt{h}{{\left(h{{h}_{0}}-1+\sqrt{\left({{h}^{2}}-1\right)\left(h_{0}^{2}-1\right)}\right)}^{2}}}{4\sqrt{\pi}{{\left(h_{0}^{2}-1\right)}^{3/4}}{{\left({{h}_{0}}h-1\right)}^{3/2}}\left(h-{{h}_{0}}\right)}~{{\ell}^{-3/2}} & \text{if}~1<{{h}_{0}}<h, \\ \sqrt{\frac{h\left({{h}_{0}}-h\right)\sqrt{h_{0}^{2}-1}}{\left(h+{{h}_{0}}\right)\left(h{{h}_{0}}-1\right)}} & \text{if}~1<h<{{h}_{0}}. \end{array}\right. \end{eqnarray} \tag{ C.6 }$$