Pre-relaxation in weakly interacting models

The emergence of stationary behaviour in closed quantum many-body systems is perhaps one of the most striking features of non-equilibrium dynamics. Rather generally subsystems behave as they were in a 'bath' and correlation functions relax to stationary values that can be described in a statistical fashion.

The maturation and refinement of experimental techniques have led to the design of experiments ever more effective in the extraction of information on the long time dynamics [1–10]. A sensitive communication between theory and experiment has then made it possible to identify the most interesting aspects of the non-equilibrium time evolution in quantum many-body systems [11].

In particular, it was realised that integrable models behave differently from generic ones. The stationary properties are indeed affected by the local conservation laws, which in integrable models are infinite in number. This led to the concept of generalised Gibbs ensemble (GGE) [12–49], which is often defined as the mixed state with maximal entropy under the constraints of the local conservation laws. Non-integrable models with no other local conservation laws except for the Hamiltonian itself are supposed to 'thermalise' at some effective temperature [50–56], whereas integrable models retain infinite information about the initial state.

In the same way as relaxation and thermalisation were associated with the late time behaviour in integrable and non-integrable models, pre-thermalisation [3, 57–61] has been recognised as a typical feature of generic models close to integrable points. Essentially, at intermediate times the non-purely-elastic processes typical of non-integrable models are almost absent and the system behaves as if it were integrable. Despite the strenuous efforts to understand the process of thermalisation in the presence of pre-thermalisation plateaux, the picture is still far from being clear and, so far, only the earliest plateau has found satisfactory descriptions [59, 61].

This state of affairs boosted the research into conserved and quasi-conserved operators in non-integrable models [62, 63], on the one hand, and put physicists' ingenuity to the test to propose sufficiently simple models to study pre-thermalisation [60, 61], on the other.

One interesting proposal [60] was to break the integrability of the transverse field Ising chain (TFIC) by adding a highly nonlocal interaction proportional to the global magnetisation squared per unit length. Even though the model possesses infinite local conservation laws (odd under reflection symmetry) [19, 64], it was argued to behave like a non-integrable model in the sector of reflection symmetric states. In particular, [60] developed a perturbation theory that allows one to follow the time evolution of the ground state of a TFIC for sufficiently long times to see a pre-thermalisation plateau.

In fact, the situation seems to be more complicated. Some techniques that will be developed in this paper allow us to analytically study the dynamics of a class of nonlocal Hamiltonians that includes the model introduced in [60]. The method we use is exact in the thermodynamic limit (in which the large-system limit is taken first): we prove some conjectures of [60], but we find that the time evolution does not result in thermalisation. Therefore, in the thermodynamic limit, models like the one of [60] can not be naively related to the physics underlying pre-thermalisation. Nevertheless, we show that similar types of nonlocal Hamiltonians can emerge at intermediate times as effective descriptions of the dynamics generated by local Hamiltonians. Thus, instead of spoiling the interest in such models, our findings give new motivations for their study.

Our main goal is to investigate the time evolution of local observables under Hamiltonians with local interactions in the particular time windows where such effective descriptions can be used. There are indeed interesting cases where the expectation values start moving significantly from a plateau that could have been approximately described by the stationary state of the unperturbed model. We show that the crossover is driven by the presence of infinitely many local conservation laws that do not commute with one another, which will be referred to as non-abelian integrability. We therefore identify two necessary requirements for a nontrivial time evolution. First, the unperturbed model must have a non-abelian set of local conservation laws. Second, the perturbation must break non-abelian integrability.

Since the crossover appears also in the presence of perturbations preserving integrability, we call it 'pre-relaxation' and the limit 'pre-relaxation limit'.

The pre-relaxation limit has been already considered in [65], where a typical crossover behaviour between two plateaux has been identified in noninteracting models like the XY quantum spin chain. Reference [65] also obtained similar results for a particular quench in an interacting model (XXZ spin-$\frac{1}{2}$ chain). However, despite an effective description was proposed that is supposed to capture the relaxation process for quite general interactions, the idea was tested only on simple cases in which the dynamics is essentially noninteracting.

In this paper we start filling this gap by investigating the pre-relaxation behaviour triggered off by more general (interacting) perturbations. This is a highly non trivial generalisation. The first universal picture of the time evolution of correlation functions after a quantum quench was delineated in conformal field theories [45, 66], but most of the analytic results have been in fact obtained in models that can be mapped to free fermions or bosons [14, 21, 27, 32–35, 47, 65, 67–69]. For the serious complications introduced by the interactions, there are far less examples [70–72] in which the time evolution of some nontrivial observable has been worked out in interacting models. In the pre-relaxation limit some obstacles can be overcome. In particular we show that, at the leading order in the perturbation strength, the dynamics generated by (local) weakly interacting Hamiltonians are equivalent to those generated by time-dependent (quasi-)local (mean-field) Hamiltonians, which can be solved in a self-consistent way. Differently from the common situations, the mean-field mapping presented here is not an uncontrolled approximation, but arises naturally in the timescale investigated under few reasonable assumptions. The possibility to write a compact system of nonlinear differential equations for the time evolution of local observables can therefore be used to investigate the essential aspects of the pre-relaxation limit even analytically.

1.1. Organisation of the paper

We consider the time evolution of some initial state $\mid {{\Psi}_{0}}\rangle$ with cluster decomposition properties³ under translation invariant Hamiltonians of the form

$$\begin{eqnarray}&&H={{H}_{0}}+gV,\end{eqnarray} \tag{ 2.1 }$$

where V is a global perturbation and g is a small coupling constant.

We focus on perturbations V that break some symmetries of H₀ in such a way that the limit of infinite time of the expectation value of a local observable $\mathcal{O}$ does not commute with the limit of infinitesimal g

$$\begin{eqnarray}&&\underset{g\to 0}{{\lim}}\,\underset{t\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}\left({{H}_{0}}+gV\right)t}}\mathcal{O}{{\text{e}}^{-\text{i}\left({{H}_{0}}+gV\right)t}}\mid {{\Psi}_{0}}\rangle \ne \underset{t\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}{{H}_{0}}t}}\mathcal{O}{{\text{e}}^{-\text{i}{{H}_{0}}t}}\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 2.2 }$$

This is the typical situation in which local degrees of freedom experience a pre-thermalisation/pre-relaxation behaviour. Indeed, at not too large times the effect of the perturbation is negligible and the expectation value has time to settle at the stationary value of the unperturbed Hamiltonian. On the other hand, at later times the perturbation can not be ignored anymore and the expectation value varies with a typical timescale that depends on the perturbation strength. Importantly, the amplitude of the variation is $\mathcal{O}\left({{g}^{0}}\right)$, therefore the pre-relaxation behaviour can be understood in the limit of infinitesimal g.

In principle, there could be many pre-relaxation plateaux, depending on how the time t scales with the small parameter g. Here we focus on the limit $g\ll 1$ and large time in such a way that $T=gt\sim O\left({{g}^{0}}\right)$. Figure 1 summarises the various steps of the formalism that will be developed in the next three sections to investigate the pre-relaxation limit. The analysis of explicit examples will be instead carried out in the last two sections.

Zoom In Zoom Out Reset image size

We now present a comprehensive summary of the main results.

• [Sections 3 and 4] Under some assumptions, the pre-relaxation limit of weakly interacting local Hamiltonians can be described by effective nonlocal Hamiltonians of the form
  $$\begin{eqnarray}&&\bar{H}=\frac{1}{{{L}^{{{n}_{1}}-1}}}Q_{1}^{(1)}\ldots Q_{{{n}_{1}}}^{(1)}+\cdots +\frac{1}{{{L}^{{{n}_{m}}-1}}}Q_{1}^{(m)}\ldots Q_{{{n}_{m}}}^{(m)},\end{eqnarray} \tag{ 2.3 }$$
  where $Q_{j}^{(\ell )}$ are (quasi-)local conservation laws⁴ of the (integrable) unperturbed model. In order to be a nontrivial limit, $Q_{j}^{(\ell )}$ should not commute with one another.
• [Section 5 and appendix C] In the thermodynamic limit, the time evolution of a state with cluster decomposition properties under Hamiltonians of the form (2.3) is completely equivalent to the time evolution under the time-dependent mean-field Hamiltonian
  $$\begin{eqnarray}&&\bar{H}_{\text{MF}}^{{{\Psi}_{0}}}(gt)=\underset{i=1}{\overset{m}{\sum}}\,\underset{j=1}{\overset{{{n}_{i}}}{\sum}}\,c_{j}^{(i)}\left(gt;{{\Psi}_{0}}\right)Q_{j}^{(i)},\end{eqnarray} \tag{ 2.4 }$$
  where the coefficients $c_{j}^{(i)}\left(gt;{{\Psi}_{0}}\right)$ are obtained self-consistently as follows:
  $$\begin{eqnarray}&&c_{j}^{(i)}\left(gt;{{\Psi}_{0}}\right)=\underset{n\ne j}{\prod}\,\text{Tr}\,\left[{{\rho}_{\text{GGE}}}\bar{U}_{\text{MF}}^{\dagger}\left(gt;{{\Psi}_{0}}\right)\frac{Q_{n}^{(i)}}{L}{{\bar{U}}_{\text{MF}}}\left(gt;{{\Psi}_{0}}\right)\right].\end{eqnarray} \tag{ 2.5 }$$
  Here ${{\rho}_{\text{GGE}}}$ is the generalised Gibbs ensemble that emerges in the time evolution under the unperturbed Hamiltonian and ${{\bar{U}}_{\text{MF}}}\left(gt;{{\Psi}_{0}}\right)$ is the time evolution operator of $\bar{H}_{\text{MF}}^{{{\Psi}_{0}}}(gt)$
  $$\begin{eqnarray}&&{{\bar{U}}_{\text{MF}}}\left(gt;{{\Psi}_{0}}\right)=\text{T}\exp\left(-\text{i}{\int}_{0}^{gt}\text{d}\tau \bar{H}_{\text{MF}}^{{{\Psi}_{0}}}(\tau )\right);\end{eqnarray} \tag{ 2.6 }$$
  $\text{T}$ is the time-ordering operator, which formally orders operators depending on $\tau$ in such a way that those on the left are associated with larger (or equivalent) values of $\tau$.Despite mean-field approximations being extremely common also in the field of non-equilibrium physics [73, 74], in the thermodynamic limit the mean-field mapping presented here is actually exact under a mild assumption.
• [Section 6] We consider Slater determinant initial states evolving under
  $$\begin{eqnarray}&&H=J\underset{\ell}{\sum}\,\left(\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}+\frac{g}{4}\sigma _{\ell}^{z}\sigma _{\ell +1}^{z}+\frac{gU}{4}\sigma _{\ell}^{z}\sigma _{\ell +2}^{z}\right)+\frac{gh}{2}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}\end{eqnarray} \tag{ 2.7 }$$
  in the the pre-relaxation limit $gt\sim O\left({{g}^{0}}\right)$ with $g\ll 1$. We identify three different behaviours:
  • Emergence of a second plateau that can be described by a GGE constructed with the local conservation laws of the unperturbed model. The information about the initial state is not encoded in a finite number of parameters because the effective Hamiltonian commutes with the unperturbed Hamiltonian and with infinitely many other conservation laws in involution. If broken, one-site shift invariance is generally not restored.
  • Oscillatory behaviour: the expectation values of local observables keep oscillating with frequency proportional to g. Nevertheless, they can be described by a time-dependent GGE (that is not one-site shift invariant).
  • In the pre-relaxation limit the expectation values are independent of time. This happens whenever the initial state is one-site shift invariant, or, more generally, when it is an excited state of $\bar{H}_{\text{MF}}^{{{\Psi}_{0}}}(0)$ (2.4) (notice that the very definition of $\bar{H}_{\text{MF}}^{{{\Psi}_{0}}}$ depends on the state, so there are implicit self-consistent conditions to be satisfied). In the latter case, one-site shift invariance is generally broken.
  We stress that, rather unexpectedly, a genuine one-site shift invariant interaction is not always sufficient to induce the restoration of one-site shift invariance. This is an indication that non-abelian integrability can survive interacting perturbations (at least at the lowest orders of perturbation theory).We will report explicit examples of the aforementioned behaviours.
• [Section 7] The equivalence with the mean-field description does not rely on the fact that the operators $Q_{j}^{(\ell )}$ of (2.3) commute with the unperturbed Hamiltonian; relaxing such hypothesis allows us to construct more general models. Reference [60] has recently proposed the model with Hamiltonian
  $$\begin{eqnarray}&&H(g,\lambda )=-\underset{\ell}{\overset{L}{\sum}}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+g\sigma _{\ell}^{z}\right)+\frac{\lambda}{L}{{\left({\sum}_{\ell}^{L}\sigma _{\ell}^{z}-\overline{{\sum}_{\ell}^{L}\sigma _{\ell}^{z}}\right)}^{2}}\end{eqnarray} \tag{ 2.8 }$$
  as a convenient framework for studying pre-thermalisation/thermalisation issues. Here $\overset{\overline{}}{{\cdots}}\,$ denotes the time average with respect to $H(g,0)$. In the thermodynamic limit $L\to \infty$, we show that the time evolution does not result in thermalisation.The apparent conflict between our results and [60] can be traced back to the different order of limits. We indeed investigate the expectation values of local observables in the limit⁵
  $$\begin{equation} (\lim_{t\rightarrow\infty})\lim_{L\rightarrow\infty}\langle{\Psi_0\vert {\rm e}^{{\rm i}H(\tilde g,\lambda) t}\mathcal O{\rm e}^{-{\rm i}H(\tilde g,\lambda) t}\vert \Psi_0}\rangle, \end{equation} \tag{ 2.9 }$$
  while, for what concerns the stationary properties, [60] considered the time average in finite systems
  $$\begin{eqnarray}&&\underset{T\to \infty}{{\lim}}\,\frac{1}{T}{\int}_{0}^{T}\text{d}t\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}H\left(\tilde{g},\lambda \right)t}}\mathcal{O}{{\text{e}}^{-\text{i}H\left(\tilde{g},\lambda \right)t}}\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 2.10 }$$
  We find that the stationary behaviour of local observables in the limit (2.9) is not characterised by a finite number of parameters. More generally, we argue that thermal-like behaviour can be excluded when every linear combination of the operators $Q_{j}^{(\ell )}$ is the Hamiltonian of an integrable model.
• [Section 7] For Hamiltonians like (2.12) one can generally identify 'critical regions' for which the late time dynamics is unstable under a small change of the parameters: an infinitesimal variation can lead both to relaxation and to persistent oscillatory behaviour (see e.g. [73] for similar discussions in quantum field theories). In addition, the variance of the expectation value of some local operator in an arbitrarily large time window
  $$\begin{eqnarray}&&\Delta \mathcal{O}=\underset{T\to \infty}{{\lim}}\,{{\left(\frac{1}{T}{\int}_{T}^{2T}\text{d}t{{\langle \mathcal{O}(t)\rangle}^{2}}-{{\left(\frac{1}{T}{\int}_{T}^{2T}\text{d}t\langle \mathcal{O}(t)\rangle \right)}^{2}}\right)}^{1/2}}\end{eqnarray} \tag{ 2.11 }$$
  behaves like an 'order parameter' for the transition, indeed it is not analytic at the boundaries of the relaxation region and vanishes inside.The crossover between relaxation and persistent oscillations can be illustrated by means of diagrams that depict the relaxation properties as a function of the Hamiltonian parameters; in particular here we focus on the limit of small quench (see also [21] for some clarifications about the meaning of 'small quench'). To highlight the close connection with the dephasing mechanisms [18] that allow local relaxation we call the diagrams 'quench dephasing diagrams'⁶. We point out that restricting ourselves to small quenches makes it easier to interpret the results in terms of low lying excitations.We investigated the quench dephasing diagram of a simplified version of (2.8), i.e.
  $$\begin{eqnarray}&&H\left(\tilde{g},\lambda \right)=-\underset{\ell}{\overset{L}{\sum}}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\tilde{g}\sigma _{\ell}^{z}\right)+\frac{\lambda}{L}{{\left({\sum}_{\ell}^{L}\sigma _{\ell}^{z}\right)}^{2}},\end{eqnarray} \tag{ 2.12 }$$
  in the limit in which the initial state is almost the ground state of the Hamiltonian. In figure 2 we can identify a 'critical' piecewise smooth curve at the boundaries of the relaxation region (the bright area). Generally, as we move towards the critical lines from the inside of the regions with persistent oscillations, the variance (2.11) of the transverse magnetisation approaches zero linearly with the distance in the $\left(\tilde{g},\lambda \right)$ parameter space (see figure 10). Close to the Ising critical point ($\tilde{g}=1$ and $\lambda =0$) the region of persistent oscillatory behaviour degenerates into a line ending at the critical point.Finally we show that persistent oscillatory behaviour can be related to the emergence of localised excitations.

Zoom In Zoom Out Reset image size

The time evolution operator for the Hamiltonian (2.1) can be formally written as follows

$$\begin{eqnarray}&&{{\text{e}}^{-\text{i}Ht}}={{U}_{\text{I}}}(T){{\text{e}}^{-\text{i}{{H}_{0}}t}},\end{eqnarray} \tag{ 3.1 }$$

where T   =   gt,

$$\begin{eqnarray}&&{{U}_{\text{I}}}(T)=\text{T}\dagger \exp\left(-\text{i}{\int}_{0}^{T}\text{d}\tau {{\text{e}}^{-\text{i}\frac{{{H}_{0}}}{g}\tau}}V{{\text{e}}^{\text{i}\frac{{{H}_{0}}}{g}\tau}}\right),\end{eqnarray} \tag{ 3.2 }$$

and $\text{T}\dagger$ is the anti-time-ordering operator. ${{U}_{\text{I}}}(T)$ can be interpreted as the evolution backwards in time under the time-dependent effective Hamiltonian

$$\begin{eqnarray}&&V(\tau )={{\text{e}}^{-\text{i}\frac{{{H}_{0}}}{g}\tau}}V{{\text{e}}^{\text{i}\frac{{{H}_{0}}}{g}\tau}}.\end{eqnarray} \tag{ 3.3 }$$

In the pre-relaxation limit, $\tau$ is finite while g is infinitesimal; it is therefore convenient to isolate the stationary (i.e. diagonal) contributions from $V(\tau )$

$$\begin{eqnarray}&&V(\tau )=\bar{V}+\delta V(\tau ),\end{eqnarray} \tag{ 3.4 }$$

where $\bar{V}$ can be formally written as follows

$$\begin{eqnarray}&&\bar{V}=\underset{t\to \infty}{{\lim}}\,\frac{1}{t}{\int}_{0}^{t}\text{d}\tau V(\tau ).\end{eqnarray} \tag{ 3.5 }$$

Let us rewrite $\delta V(\tau )$ in terms of its Fourier transform $\delta \tilde{V}(\varepsilon )$

$$\begin{eqnarray}&&\delta V(\tau )={\int}_{-\infty}^{\infty}\text{d}\varepsilon {{\text{e}}^{\text{i}\frac{\varepsilon \tau}{g}}}\delta \tilde{V}(\varepsilon )\quad \delta {{\tilde{V}}^{\dagger}}(\varepsilon )=\delta \tilde{V}(-\varepsilon ).\end{eqnarray} \tag{ 3.6 }$$

We notice that for noninteracting H₀ (which is the case we are going to consider), the locality of V (in the fermionic picture) implies that $\delta \tilde{V}(\varepsilon )$ is zero for $\mid \varepsilon \mid >{{\varepsilon}_{\text{max}}}\sim O(1)$. After some formal manipulations we find

$$\begin{eqnarray}&&{{U}_{\text{I}}}(T)\left(\text{I}-\text{i}gA(T)\right)=\text{I}-\text{i}gA(0)-\text{i}{\int}_{0}^{T}\text{d}\tau {{U}_{\text{I}}}(\tau )\left[\bar{V}-\text{i}g\delta V(\tau )A(\tau )\right]\end{eqnarray} \tag{ 3.7 }$$

with

$$\begin{eqnarray}&&A(\tau )=\text{i}{\int}_{-\infty}^{\infty}\text{d}\varepsilon {{\text{e}}^{\text{i}\frac{\varepsilon \tau}{g}}}\frac{\delta \tilde{V}(\varepsilon )}{\varepsilon}.\end{eqnarray} \tag{ 3.8 }$$

Equation (3.7) makes sense as long as $A(\tau )$ does. In particular, if $A(\tau )$ per unit length is a bounded operator, in the limit $g\to 0$ all the terms of (3.7) that are multiplied by g can be neglected, i.e. $\delta V(\tau )$ is negligible. More generally, it is sufficient that the matrix elements of $A(\tau )$ that give a relevant contribution in the expectation values of local observables are bounded. In [65] the operator $\delta V(\tau )$ was worked out for a particular noninteracting perturbation and the previous assumption turned out to be satisfied. Besides the noninteracting case, the irrelevance of $\delta V(\tau )$ was also implicitly assumed in [61]. There, using the continuous unitary transformation (CUT) formalism [57, 77], the authors worked out the dynamics after a quantum quench in a weakly interacting model and checked the results against numerical data. For the readers familiar with CUT, we indeed notice that the formal simplification of the terms proportional to g in (3.7) is equivalent to CUT at O(g⁰): At the lowest order of perturbation theory the CUT unitary transformation can be replaced by the identity; a residual dependence on g remains in the CUT Hamiltonian ${{H}_{\text{CUT}}}$, because we are considering the limit of infinite time with T   =   gt finite, so ${{H}_{\text{CUT}}}$ must be computed at O(g). The latter is simply the time average of the Hamiltonian over its unperturbed part, namely ${{H}_{\text{CUT}}}={{H}_{0}}+g\bar{V}+o(g)$. The excellent agreement of [61] with the numerical simulations suggests that the interaction can be replaced by its time average for rather general perturbations. Therefore, from now on we shall assume that the effect of $\delta V(\tau )$ is negligible in the pre-relaxation limit.

Going back to (3.1), if $\delta V(\tau )$ is negligible we find

$$\begin{eqnarray}&&{{\text{e}}^{-\text{i}Ht}}\to {{\text{e}}^{-\text{i}T\bar{V}}}{{\text{e}}^{-\text{i}{{H}_{0}}T/g}},\end{eqnarray} \tag{ 3.9 }$$

where $\left[\bar{V},{{H}_{0}}\right]=0$. We now consider the time evolution of the expectation value of some local operator $\mathcal{O}$

$$\begin{eqnarray}&&\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}Ht}}\mathcal{O}{{\text{e}}^{-\text{i}Ht}}\mid {{\Psi}_{0}}\rangle \to \langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}{{H}_{0}}T/g}}{{\text{e}}^{\text{i}T\bar{V}}}\mathcal{O}\,{{\text{e}}^{-\text{i}T\bar{V}}}\,{{\text{e}}^{-\text{i}{{H}_{0}}T/g}}\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 3.10 }$$

It is well established that the stationary properties of (quasi-)local observables after quenches in translation invariant noninteracting models from states with cluster decomposition properties can be described by means of a generalised Gibbs ensemble (GGE) of the form

$$\begin{eqnarray}&&{{\rho}_{\text{GGE}}}=\frac{{{\text{e}}^{-{\sum}_{j}{{\lambda}_{j}}{{Q}_{j}}}}}{Z},\end{eqnarray} \tag{ 3.11 }$$

where Q_j are local conservation laws and ${{\lambda}_{j}}$ are real parameters determined by the initial state [19, 34]. We also remind the reader that, in order to avoid an explicit dependence of the charges on the initial state, in some special cases the set of charges {Q_j} could be non-abelian [65].

Let us now assume that the perturbation is sufficiently 'nice' that the late time dynamics of ${{\text{e}}^{\text{i}T\bar{V}}}\mathcal{O}{{\text{e}}^{-\text{i}T\bar{V}}}$ under H₀ can be obtained by replacing the state with the corresponding GGE

$$\begin{eqnarray}&&{{\text{e}}^{-\text{i}{{H}_{0}}T/g}}\,\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid \,{{\text{e}}^{\text{i}{{H}_{0}}T/g}}\to {{\rho}_{\text{GGE}}}=\underset{\mid S\mid \to \infty}{{\lim}}\,\underset{t\to \infty}{{\lim}}\,\text{t}{{\text{r}}_{{\bar{S}}}}\left[{{\text{e}}^{-\text{i}{{H}_{0}}t}}\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}{{H}_{0}}t}}\right].\end{eqnarray} \tag{ 3.12 }$$

This step can be easily justified for $T=gt\ll t$ if $\bar{V}$ is (quasi-)local, however in the next sections we will show that, in the presence of interactions, $\bar{V}$ belongs to a larger class of operators, so it is convenient to postpone the explanation of (3.12) after having clarified the properties of $\bar{V}$.

From (3.12) it follows

$$\begin{eqnarray}&&\underset{g\to 0}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}{{H}_{0}}T/g}}\,{{\text{e}}^{\text{i}T\bar{V}}}\mathcal{O}{{\text{e}}^{-\text{i}T\bar{V}}}\,{{\text{e}}^{-\text{i}{{H}_{0}}T/g}}\mid {{\Psi}_{0}}\rangle =\text{tr}\left[{{\rho}_{\text{GGE}}}{{\text{e}}^{\text{i}T\bar{V}}}\,\mathcal{O}{{\text{e}}^{-\text{i}T\bar{V}}}\right],\end{eqnarray} \tag{ 3.13 }$$

which suggests that the pre-relaxation limit can be described by the time-dependent ensemble

$$\begin{eqnarray}&&{{\rho}_{\text{tGGE}}}(t)={{\text{e}}^{-\text{i}\bar{V}gt}}{{\rho}_{\text{GGE}}}{{\text{e}}^{\text{i}\bar{V}gt}},\end{eqnarray} \tag{ 3.14 }$$

where we re-expressed the rescaled time T   =   gt in terms of the time.

It is important to note that both ${{\rho}_{\text{GGE}}}$ and $\bar{V}$ commute with the Hamiltonian. Consequently, if the two operators can be written in terms of the same set of local conservation laws in involution, the time dependence disappears ${{\rho}_{\text{tGGE}}}(t)={{\rho}_{\text{GGE}}}$. In the next section we will show that in many cases of interest $\bar{V}$ can be approximated by a polynomial of the local conservation laws. Therefore, in order to see some nontrivial pre-relaxation behaviour, the unperturbed Hamiltonian H₀ must have a non-abelian set of local charges. We refer the reader to [65] for an extensive discussion of noninteracting models with that property.

The simplest noninteracting model that possesses local conservation laws that are not mutually commuting is the XY model, whose Hamiltonian is given by

$$\begin{eqnarray}&&{{H}_{\text{XY}}}=J\underset{\ell}{\sum}\,\left(\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}\right),\end{eqnarray} \tag{ 4.1 }$$

where $\sigma _{\ell}^{\alpha}$ act like Pauli matrices on the site $\ell$ and like the identity elsewhere. If the initial state $\mid {{\Psi}_{0}}\rangle$ breaks one-site shift invariance, the latter symmetry is generally not restored in the GGE that describes local observables at infinite time after the quench.

On the other hand, an infinitesimally small one-site shift invariant perturbation that breaks the non-abelian integrability of (4.1) is expected to catalyse symmetry restoration, which may be captured by the pre-relaxation limit. Similar issues of symmetry restoration have been pointed out long ago, e.g. in [78].

A perturbation that preserves the noninteracting character of the Hamiltonian was already considered in [65]. Here we investigate perturbations that have a 4-fermion representation in terms of the noninteracting fermions that diagonalise (4.1), namely

$$\begin{eqnarray}&&V\sim \underset{\ell}{\sum}\,a_{\ell +{{n}_{1}}}^{{{\alpha}_{1}}}a_{\ell +{{n}_{2}}}^{{{\alpha}_{2}}}a_{\ell +{{n}_{3}}}^{{{\alpha}_{3}}}a_{\ell +{{n}_{4}}}^{{{\alpha}_{4}}},\end{eqnarray} \tag{ 4.2 }$$

where $a_{\ell}^{\alpha}$ are the Majorana fermions ($\left\{a_{\ell}^{\alpha},a_{n}^{\beta}\right\}=2{{\delta}_{\alpha \beta}}{{\delta}_{\ell n}}$)

$$\begin{eqnarray}&&a_{\ell}^{\alpha}=\left(\underset{j<\ell}{\prod}\,\sigma _{j}^{z}\right)\sigma _{\ell}^{\alpha}\quad \alpha \in \left\{x,y\right\}.\end{eqnarray} \tag{ 4.3 }$$

From the qualitative argument presented in section 3, the relevant Hamiltonian in the pre-relaxation limit is determined by the time average of the perturbation. The calculation is not difficult but rather lengthy. However, a close inspection of the various contributions reveals a hidden structure that helps simplifying the computation. We indeed find (see appendix B).

Property 4.1. The time average under ${{H}_{\text{XY}}}$ of a one-site shift invariant four fermion operator can be written as follows

$$\begin{eqnarray}&&\frac{1}{L}\overline{\underset{\ell}{\sum}\,a_{\ell +{{n}_{1}}}^{{{\alpha}_{1}}}a_{\ell +{{n}_{2}}}^{{{\alpha}_{2}}}a_{\ell +{{n}_{3}}}^{{{\alpha}_{3}}}a_{\ell +{{n}_{4}}}^{{{\alpha}_{4}}}}=F_{\left\{n\right\}}^{\left\{\alpha \right\}}+A_{\left\{n\right\}}^{\left\{\alpha \right\}},\end{eqnarray} \tag{ 4.4 }$$

where F is a linear combination of factorised terms and A is an anomalous contribution originated by the nontrivial solutions of the energy constraint

$$\begin{eqnarray}&&\varepsilon \left({{k}_{1}}\right)+\varepsilon \left({{k}_{2}}\right)=\varepsilon \left({{k}_{3}}\right)+\varepsilon \left({{k}_{1}}+{{k}_{2}}+{{k}_{3}}\right).\end{eqnarray} \tag{ 4.5 }$$

The latter exists only in the thermodynamic limit and strongly depends on the details of the dispersion relation, whereas F has a structure that is almost model independent:

$$\begin{equation} F_{\{n\}}^{\{\alpha\}}=\sum_{s=0}^1\underbrace{a_{n_1}^{\alpha_1} a_{n_2}^{\alpha_2}}_{s}\underbrace{a_{n_3}^{\alpha_3} a_{n_4}^{\alpha_4}}_{s}-\underbrace{a_{n_1}^{\alpha_1} a_{n_3}^{\alpha_3}}_{s}\underbrace{a_{n_2}^{\alpha_2} a_{n_4}^{\alpha_4}}_{s}+\underbrace{a_{n_1}^{\alpha_1} a_{n_4}^{\alpha_4}}_{s}\underbrace{a_{n_2}^{\alpha_2} a_{n_3}^{\alpha_3}}_{s}, \end{equation} \tag{ 4.6 }$$

with

$$\begin{equation}\underbrace{a_{n_1}^{\alpha} a_{n_2}^{\beta}}_s=\overline{ \frac{1}{L}\sum_\ell (-1)^{s\ell}a_{\ell+n_1}^{\alpha} a_{\ell+n_2}^{\beta} }. \end{equation} \tag{ 4.7 }$$

Essentially, index s appears because the XY model (with zero magnetic field) has local conservation laws with momentum $\pi$, while one-site shift invariance constrains the total momentum to be multiple of $2\pi$. We also notice that the factorised part of the time average can be easily generalised to an arbitrary number of fermions, keeping the same structure of the Wick decomposition.

To the best of our knowledge, the (quasi-)local conservation laws of the XY model (4.1) are noninteracting (for $\mid \gamma \mid \ne 1$) and $A_{\left\{n\right\}}^{\left\{\alpha \right\}}$ of (4.6) seems to be a nonlocal conservation law that can not be (not even approximately) written as a function of the local charges. We then expect $A_{\left\{n\right\}}^{\left\{\alpha \right\}}$ to become important (for local observables) only at times proportional to the chain length, which are far beyond the pre-relaxation limit. This persuaded us to conjecture that the anomalous terms are not relevant to our problem, which is equivalent to assume

$$\begin{eqnarray}&&\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}Ht}}\left[\mathcal{O},A\right]{{\text{e}}^{-\text{i}Ht}}\mid {{\Psi}_{0}}\rangle =0\end{eqnarray} \tag{ 4.8 }$$

for any local observable $\mathcal{O}$. In appendix D we check the self-consistency of our approximation, showing that it is compatible with (4.8). We leave further investigations to future works.

Proposition 4.1 suggests that in many cases of interest the effective Hamiltonian describing the pre-relaxation limit takes the form

$$\begin{eqnarray}&&{{H}_{\text{eff}}}=\frac{1}{{{L}^{{{n}_{1}}-1}}}H_{1}^{(1)}\ldots H_{{{n}_{1}}}^{(1)}+\cdots +\frac{1}{{{L}^{{{n}_{m}}-1}}}H_{1}^{(m)}\ldots H_{{{n}_{m}}}^{(m)},\end{eqnarray} \tag{ 4.9 }$$

where $H_{j}^{(\ell )}$ are (quasi-)local (i.e. their density is (quasi-)local, see e.g. [63, 79]) translation invariant operators (i.e. n-site shift invariant for some $n\in \mathbb{N}$). Indeed, provided that the anomalous terms in (4.4) can be disregarded, similar factorisations appear whenever the unperturbed Hamiltonian is noninteracting (e.g., in the model considered in [61]).

Hamiltonians of the form (4.9) are therefore the perfect workbench for pre-relaxation or pre-thermalisation issues.

We notice that the non-equilibrium dynamics generated by a subclass of Hamiltonians of the form (4.9) have been already worked out in [75]. The authors considered 'completely connected quantum models', in which the Hamiltonian is symmetric under any permutation of the sites, and exhibited a mapping onto an effective classical Hamiltonian dynamics.

We also point out that the simplest models of the form (4.9) (e.g. Curie-Weiss quantum Heisenberg models) have often been used as toy models to investigate the statistical properties in the presence of long range interactions [80].

The rest of the paper will be focussed on the following points:

• (i)  
  Solution of the non-equilibrium problem for Hamiltonians of the form (4.9);
• (ii)  
  Characterisation of the pre-relaxation limit in an interacting model, also in the presence of perturbations that break integrability;
• (iii)  
  Non-equilibrium time evolution under (2.12).

For the sake of clarity, we stress again that (ii) relies on two assumptions:

• (a)  
  In the limit $g\to 0$ with gt finite, the time evolution under H   =   H₀ + gV can be split in two steps:

  • 1.  
    infinite time evolution under the unperturbed Hamiltonian H₀, which is supposed to give rise to a generalised Gibbs ensemble $\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid \to {{\rho}_{\text{GGE}}}$;
  • 2.  
    time evolution with rescaled time T   =   gt under the effective Hamiltonian given by the perturbation V averaged with respect to H₀ (3.5);

    $$\begin{eqnarray}&&{{\text{e}}^{-\text{i}\left({{H}_{0}}+gV\right)t}}\,\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid \,{{\text{e}}^{\text{i}\left({{H}_{0}}+gV\right)t}}\sim {{\text{e}}^{-\text{i}gt\bar{V}}}\,{{\rho}_{\text{GGE}}}{{\text{e}}^{\text{i}gt\bar{V}}}.\end{eqnarray} \tag{ 4.10 }$$
• (b)  
  The 'anomalous terms' that appear in the time average of V give a negligible contribution (see appendix D, property 4.1 and discussion below).

On the other hand, (i) and (iii) will be treated as ab initio problems.

In this section we work out problem (i). We are going to show that, despite the nonlocal appearance, operators of the form (4.9) generate a dynamics which is equivalent to that of a (quasi-)local time-dependent mean-field Hamiltonian.

Here we only report some results and three useful corollaries, the details of the derivation can be found in appendix C.

For the sake of simplicity we only consider cases in which $H_{j}^{(\ell )}$ have local densities, however, as far as we can see, all the results can be generalised to quasi-local operators with tails that decay exponentially with the distance.

In the light of (4.9), we define a class of operator $\mathcal{E}$ as follows:

Definition 1. We say that an operator acting on a spin-$\frac{1}{2}$ chain belongs to the class $\mathcal{E}$ if it is written as in (4.9), namely as a finite linear combination of operators of the form

$$\begin{eqnarray}&&\frac{1}{{{L}^{n-1}}}{{H}_{1}}\cdots {{H}_{n}},\end{eqnarray} \tag{ 5.1 }$$

where n is finite, H_j are local translation invariant operators, and L is the chain length.

We consider spin chains so that the local Hilbert space is finite dimensional. This turns out to be a fundamental assumption for most of our results.

One of our goals is to show that the time evolution preserves cluster decomposition properties, which is the key element that allows us to simplify the calculation of expectation values. For example we have

Lemma 5.1. Let $\mathcal{O}\in \mathcal{E}$ and $\mid \Psi\rangle$ a state with cluster decomposition properties. The expectation value of $\mathcal{O}/L$ in $\mid \Psi\rangle$ can be reduced to the expectation values of the local translation invariant operators it consists of:

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\,\langle \Psi\mid \frac{{{H}_{1}}}{L}\cdots \frac{{{H}_{n}}}{L}\mid \Psi\rangle =\underset{L\to \infty}{{\lim}}\,\underset{j}{\prod}\,\frac{\langle \Psi\mid {{H}_{j}}\mid \Psi\rangle}{L}.\end{eqnarray} \tag{ 5.2 }$$

Using this lemma it is rather natural to put in relation the dynamics under (4.9) to that under the mean-field Hamiltonian defined as follows:

Definition 2. Mean-field effective Hamiltonian. Let $H\in \mathcal{E}$. We define the time-dependent mean-field Hamiltonian $H_{\text{MF}}^{{{\Psi}_{0}}}(t)$ as the operator resulting from mapping any generic term (5.1) of the Hamiltonian (4.9) to an operator with local density, as follows

$$\begin{eqnarray}&&\frac{1}{{{L}^{n-1}}}{{H}_{1}}\cdots {{H}_{n}}\to \underset{j=1}{\overset{n}{\sum}}\,\underset{\ell \ne j}{\prod}\,\frac{\langle {{\Psi}_{0}}\mid {{{\bar{U}}}^{\dagger}}(t){{H}_{\ell}}\bar{U}(t)\mid {{\Psi}_{0}}\rangle}{L}{{H}_{j}},\end{eqnarray} \tag{ 5.3 }$$

where $\bar{U}(t)$ is the time evolution under $H_{\text{MF}}^{{{\Psi}_{0}}}(t)$

$$\begin{eqnarray}&&\bar{U}(t)=\text{T}\exp\left(-\text{i}{\int}_{0}^{t}\text{d}\tau H_{\text{MF}}^{{{\Psi}_{0}}}(\tau )\right).\end{eqnarray} \tag{ 5.4 }$$

Thus, generally $H_{\text{MF}}^{{{\Psi}_{0}}}(t)$ must be computed in a self-consistent way.

For example, the Hamiltonian

$$\begin{eqnarray}&&H=-\frac{1}{4}\underset{\ell}{\overset{L}{\sum}}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}\right)+\frac{\lambda}{L}{{\left({\sum}_{\ell}^{L}\sigma _{\ell}^{z}\right)}^{2}}\end{eqnarray} \tag{ 5.5 }$$

belongs to $\mathcal{E}$. In this trivial case $\underset{\ell}{\sum}\,\sigma _{\ell}^{z}$ commutes with H, so the mean-field Hamiltonian is independent of time and it is given by

$$\begin{eqnarray}&&H_{\text{MF}}^{{{\Psi}_{0}}}(t)=-\frac{1}{4}\underset{\ell}{\sum}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}\right)+2\lambda \langle {{\Psi}_{0}}\mid \frac{1}{L}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}\mid {{\Psi}_{0}}\rangle \underset{\ell}{\sum}\,\sigma _{\ell}^{z}.\end{eqnarray} \tag{ 5.6 }$$

We point out that the expectation value (per unit length) of $H\in \mathcal{E}$ in the state $\bar{U}(t)\mid {{\Psi}_{0}}\rangle$ is generally different from that of $H_{\text{MF}}^{{{\Psi}_{0}}}$:

$$\begin{eqnarray}&&\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\frac{1}{{{L}^{n}}}{{H}_{1}}\cdots {{H}_{n}}\bar{U}(t)\mid {{\Psi}_{0}}\rangle =\underset{j=1}{\overset{n}{\prod}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\frac{{{H}_{j}}}{L}\bar{U}(t)\mid {{\Psi}_{0}}\rangle\\ &&\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\underset{j=1}{\overset{n}{\sum}}\,\underset{\ell \ne j}{\prod}\,\frac{\langle {{\Psi}_{0}}\mid {{{\bar{U}}}^{\dagger}}(t){{H}_{\ell}}\bar{U}(t)\mid {{\Psi}_{0}}\rangle}{L}{{H}_{\ell}}\bar{U}(t)\mid {{\Psi}_{0}}\rangle =n\underset{j=1}{\overset{n}{\prod}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\frac{{{H}_{j}}}{L}\bar{U}(t)\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 5.7 }$$

The main property that is proved in appendix C is the exactness of the mean-field description in the thermodynamic limit:

Lemma 5.2. Let $\mid {{\Psi}_{0}}\rangle$ be a translation invariant state with cluster decomposition properties and $H,\mathcal{O}\in \mathcal{E}$. Let the expectation value of $\mathcal{O}$ in the state that time evolves with $H_{\text{MF}}^{{{\Psi}_{0}}}(t)$ be an analytic function of t in the strip $\mid \text{Im}\left[t\right]\mid <r$, with r a nonzero constant. In the thermodynamic limit, the time evolution with H can be replaced by the time evolution with the mean-field Hamiltonian:

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}Ht}}\frac{\mathcal{O}}{L}{{\text{e}}^{-\text{i}Ht}}\mid {{\Psi}_{0}}\rangle =\underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\frac{\mathcal{O}}{L}\bar{U}(t)\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 5.8 }$$

Remark 1. The validity of the hypothesis of analyticity on a strip can be verified a posteriori. The idea is the following. The self-consistent mean-field problem can be generally recast into an infinite nonlinear system of ordinary differential equations. The finiteness of n in (5.1) implies that the system can be written as $\dot{\vec u}=\vec F(\vec u, t)$, with $\vec{F}$ a polynomial. If the system was finite, the solution would have been analytic. This is not always the case for an infinite system but, in practice, the numerical solution is obtained by introducing a cutoff parameter N that makes the system finite. If the mean-field time evolution had a point of non-analyticity, the solution of the system of equations should display a non-trivial dependence of the mean-field parameters on the cutoff as $N\to \infty$.

Corollary 5.3. Lemma 5.2 holds true in particular for local operators.

The local equivalence with the mean-field time evolution can also be expressed in terms of reduced density matrices:

Corollary 5.4. Let $\mid {{\Psi}_{0}}\rangle$ a translation invariant state with cluster decomposition properties and $H\in \mathcal{E}$. In the thermodynamic limit, the time evolution of the reduced density matrix (RDM) of some spin block S is equal to the RDM in the state that time evolves with the mean-field Hamiltonian:

$$\begin{eqnarray}&&{{\rho}_{S}}(t)=\text{t}{{\text{r}}_{{\bar{S}}}}\left[{{\text{e}}^{-\text{i}Ht}}\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}Ht}}\right]=\text{t}{{\text{r}}_{{\bar{S}}}}\left[\bar{U}(t)\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\right].\end{eqnarray} \tag{ 5.9 }$$

The previous lemmas and corollaries are sufficient to reduce the time evolution under $H\in \mathcal{E}$ to the time evolution under a local time-dependent Hamiltonian. There is however another simple corollary to lemma 5.2 that will be useful to assess whether or not at large times it is possible to encode the entire information about the initial state in a finite number of parameters ('thermal-like behaviour').

Corollary 5.5. Let $H\in \mathcal{E}$ and $\mid \Psi\rangle$ a state with cluster decomposition properties. If $\mid \Psi\rangle$ is an excited state of the corresponding mean-field Hamiltonian $H_{\text{MF}}^{\Psi}$

$$\begin{eqnarray}&&H_{\text{MF}}^{\Psi}\mid \Psi\rangle ={{E}_{\Psi}}\mid \Psi\rangle ,\end{eqnarray} \tag{ 5.10 }$$

the expectation value of local observables in ${{\text{e}}^{-\text{i}Ht}}\mid \Psi\rangle$ is independent of time. Therefore, $\mid \Psi\rangle$ behaves locally as an excited state of H.

The reverse is also true. If an excited state of H is locally equivalent to a state with cluster decomposition properties, then the latter is (equivalent to) an excited state of the corresponding mean-field Hamiltonian.

5.1. Time-dependent GGE

We are now in a position to justify (3.12), and in turn (3.13) and (3.14).

In the limit of small g the expectation value of a local observable $\mathcal{O}$ reads as (3.10)

$$\begin{eqnarray}&&\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}{{H}_{0}}T/g}}{{\text{e}}^{\text{i}T\bar{V}}}\mathcal{O}{{\text{e}}^{-\text{i}T\bar{V}}}\,{{\text{e}}^{-\text{i}{{H}_{0}}T/g}}\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 5.11 }$$

Since H₀ is local, the state ${{\text{e}}^{-\text{i}{{H}_{0}}T/g}}\mid {{\Psi}_{0}}\rangle$ has cluster decomposition properties beyond some typical distance proportional to T/g (in order to be outside of the light cone). From corollary 5.3 it follows that the time evolution under $\bar{V}$ is equivalent to that under the corresponding mean-field operator. We indeed only need $JT\ll gL$ (see (C.8)), which is trivially satisfied in the thermodynamic limit. Thus we obtain

$$\begin{eqnarray}&&\langle {{\Psi}_{T}}\mid U_{{\bar{V}}}^{\dagger}(T)\mathcal{O}U_{{\bar{V}}}^{\dagger}(T)\mid {{\Psi}_{T}}\rangle ,\end{eqnarray} \tag{ 5.12 }$$

with

$$\begin{eqnarray}&&{{U}_{{\bar{V}}}}(t)=\text{T}\exp\left(-\text{i}{\int}_{0}^{t}\text{d}\tau \bar{V}_{\text{MF}}^{{{\Psi}_{T}}}(\tau )\right)\end{eqnarray} \tag{ 5.13 }$$

and

$$\begin{eqnarray}&&\mid {{\Psi}_{T}}\rangle ={{\text{e}}^{-\text{i}{{H}_{0}}T/g}}\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ 5.14 }$$

Incidentally, we notice that the time-ordering in (5.13) can not be simplified because $\bar{V}_{\text{MF}}^{{{\Psi}_{T}}}$ is generally written in terms of conservation laws that are not in involution with one another.

For the sake of simplicity we assume that the time-dependent coupling constants of $\bar{V}_{\text{MF}}^{{{\Psi}_{T}}}$ are bounded. The operator $U_{{\bar{V}}}^{\dagger}(T)\mathcal{O}U_{{\bar{V}}}^{\dagger}(T)$ is then quasi-local with a typical range $\xi$ proportional to T [81]. On the other hand $\mid {{\Psi}_{T}}\rangle$ is the time evolution of $\mid {{\Psi}_{0}}\rangle$ at the time $(\infty \leftarrow )T/g\gg T\sim \xi$, which is the limit in which it is reasonable to expect that the state can be replaced by the corresponding generalised Gibbs ensemble (of the unperturbed Hamiltonian)

$$\begin{eqnarray}&&\langle {{\Psi}_{T}}\mid U_{{\bar{V}}}^{\dagger}(T)\mathcal{O}U_{{\bar{V}}}^{\dagger}(T)\mid {{\Psi}_{T}}\rangle \sim \left[{{\rho}_{\text{GGE}}}U_{{\bar{V}}}^{\dagger}(T)\mathcal{O}U_{{\bar{V}}}^{\dagger}(T)\right].\end{eqnarray} \tag{ 5.15 }$$

The operator $\bar{V}_{\text{MF}}^{{{\Psi}_{T}}}$ is obtained self-consistently by computing the expectation values of (quasi-)local conservation laws, which can be obtained from (5.15). Therefore, in the definition (5.3) of the mean-field Hamiltonian we can replace $\mid {{\Psi}_{0}}\rangle$ by ${{\rho}_{\text{GGE}}}$

$$\begin{eqnarray}&&\frac{1}{{{L}^{n-1}}}{{H}_{1}}\cdots {{H}_{n}}\to \underset{j=1}{\overset{n}{\sum}}\,\underset{\ell \ne j}{\prod}\,\frac{\text{tr}\left[{{\rho}_{\text{GGE}}}U_{{\bar{V}}}^{\dagger}(T){{H}_{\ell}}{{U}_{{\bar{V}}}}(T)\right]}{L}{{H}_{\ell}},\end{eqnarray} \tag{ 5.16 }$$

which is consistent with (3.13). We denote by ${{\bar{H}}_{\text{MF}}}(T)$ the mean-field Hamiltonian with the expectation values computed in the GGE.

In this section we investigate the (integrable) XYZ spin-$\frac{1}{2}$ chain in the limit of small anisotropy in the z direction and also the effect of a small perturbation that breaks integrability. The Hamiltonian

$$\begin{eqnarray}&&\quad H=J\underset{\ell}{\sum}\,\left(\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}+\frac{g}{4}\sigma _{\ell}^{z}\sigma _{\ell +1}^{z}+\frac{gU}{4}\sigma _{\ell}^{z}\sigma _{\ell +2}^{z}\right)+\frac{gh}{2}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}\end{eqnarray} \tag{ 6.1 }$$

has the form (2.1) with ${{H}_{0}}={{H}_{\text{XY}}}$ (4.1) and

$$\begin{eqnarray}&&V=\frac{J}{4}\underset{\ell}{\sum}\,\left(\sigma _{\ell}^{z}\sigma _{\ell +1}^{z}+U\sigma _{\ell}^{z}\sigma _{\ell +2}^{z}\right)+\frac{h}{2}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}.\end{eqnarray} \tag{ 6.2 }$$

For a fixed g   ≠   0 the model is integrable for JU   =   h   =   0, corresponding to the spin-$\frac{1}{2}$ XYZ model, and for $JU=\gamma =0$, corresponding to the XXZ spin-$\frac{1}{2}$ chain; otherwise it is non-integrable.

Following sections 3 and 4, in the pre-relaxation limit $g\ll 1$ with $gt\sim O\left({{g}^{0}}\right)$, the initial state can be replaced by the corresponding GGE of the unperturbed Hamiltonian

$$\begin{eqnarray}&&\mid {{\Psi}_{0}}\rangle \to {{\rho}_{\text{GGE}}}=\underset{\mid S\mid \to \infty}{{\lim}}\,\underset{t\to \infty}{{\lim}}\,\text{t}{{\text{r}}_{{\bar{S}}}}\left[{{\text{e}}^{-\text{i}{{H}_{\text{XY}}}t}}\,\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid \,{{\text{e}}^{\text{i}{{H}_{\text{XY}}}t}}\right],\end{eqnarray} \tag{ 6.3 }$$

and V by the time averaged perturbation (3.5). We notice that the free Hamiltonian ${{H}_{\text{XY}}}$ does not play any role in the pre-relaxation limit, because it commutes with ${{\rho}_{\text{GGE}}}$. The mapping into a mean-field problem can be decomposed in the following steps:

–Compute the time averaged perturbation $\bar{V}$;

–Construct the mean-field Hamiltonian ${{\bar{H}}_{\text{MF}}}$;

–Solve the time evolution under ${{\bar{H}}_{\text{MF}}}$ for any local observable.

Some properties of the unperturbed Hamiltonian ${{H}_{\text{XY}}}$ dramatically simplify the first step. ${{H}_{\text{XY}}}$ is mapped to noninteracting fermions by a Jordan–Wigner transformation. Up to irrelevant (to our purposes) boundary terms, it can be written as follows

$$\begin{eqnarray}{{H}_{XY}}\sim \frac{1}{4}\underset{\ell ,n}{\sum}\,\left(\begin{array}{c} a_{2\ell -1}^{x} & a_{2\ell -1}^{y} & a_{2\ell}^{x} & a_{2\ell}^{y} \\ \end{array}\right)\left[\mathcal{H}\right]_{\ell n}^{(2)}\left(\begin{array}{c} a_{2n-1}^{x} \\ a_{2n-1}^{y} \\ a_{2n}^{x} \\ a_{2n}^{y} \\ \end{array}\right)\end{eqnarray} \tag{ 6.4 }$$

with $a_{\ell}^{\alpha}$ the Majorana fermions (4.3); $\mathcal{H}$ is the block-circulant matrix

$$\begin{eqnarray}&&\left[\mathcal{H}\right]_{\ell n}^{(2)}\sim {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}{{\text{e}}^{-\text{i}(n-\ell )k}}{{\mathcal{H}}^{(2)}}(k),\end{eqnarray} \tag{ 6.5 }$$

where the 4-by-4 matrix ${{\mathcal{H}}^{(2)}}(k)$ is usually called symbol (see also appendix A) and it is given by

$$\begin{eqnarray}&&{{\mathcal{H}}^{(2)}}(k)=-{{\varepsilon}_{k}}{{\sigma}^{x}}{{\text{e}}^{\text{i}k/2{{\sigma}^{z}}}}\otimes {{\sigma}^{y}}{{\text{e}}^{\text{i}{{\theta}_{k}}{{\sigma}^{z}}}};\end{eqnarray} \tag{ 6.6 }$$

${{\varepsilon}_{k}}$ and ${{\theta}_{k}}$ are the dispersion relation and the Bogoliubov angle, respectively

$$\begin{eqnarray}&&{{\varepsilon}_{k}}=J\sqrt{{{\cos}^{2}}k/2+{{\gamma}^{2}}{{\sin}^{2}}k/2}\quad {{\text{e}}^{\text{i}{{\theta}_{k}}}}=\frac{\cos k/2+\text{i}\gamma \sin k/2}{\sqrt{{{\cos}^{2}}k/2+{{\gamma}^{2}}{{\sin}^{2}}k/2}}.\end{eqnarray} \tag{ 6.7 }$$

Here we have chosen the two-site shift invariant representation of the Hamiltonian (i.e. we gathered together the fermionic degrees of freedom of pairs of adjacent sites) in order to be able to treat a larger class of initial states.

In translation invariant noninteracting models almost any calculation can be traced back to operations on the symbol associated with the operator, which is the Fourier transform of a block-row of the block-circulant matrix that appears in the fermionic representation of the operator as in (6.4) (see also appendix A). More generally the two-site representation of the symbol is a 4-by-4 Hermitian matrix, function of the momentum and odd under simultaneous transposition and reversion of the momentum. A 2n-by-2n symbol completely identifies a noninteracting operator that is translation invariant by k sites, with k a divisor of n, by the same kind of relations that we wrote for the Hamiltonian (i.e. (6.4)–(6.6)). Thus, we will often report the symbols instead of writing the operators explicitly.

Coming back to the calculation of $\bar{V}$, we find that the three constituents of the interaction term in (6.2) have the following fermionic representation

$$\begin{eqnarray}&&\frac{1}{4}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}\sigma _{\ell +j}^{z}=\frac{1}{4}\underset{\ell}{\sum}\,\text{i}a_{\ell}^{y}a_{\ell}^{x}\text{i}a_{\ell +j}^{y}a_{\ell +j}^{x},\quad j=1,2,\end{eqnarray} \tag{ 6.8 }$$

$$\begin{eqnarray}&&\frac{1}{2}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}=\frac{1}{2}\underset{\ell}{\sum}\,\text{i}a_{\ell}^{y}a_{\ell}^{x}.\end{eqnarray} \tag{ 6.9 }$$

Therefore, on the basis of our assumptions and decomposition (4.6), in the limit $g\ll 1$ with $gt\sim O\left({{g}^{0}}\right)$, we expect the local Hamiltonian (6.1) to be dynamically equivalent to the following nonlocal one

$$\begin{eqnarray}&&H\to \bar{H}={{H}_{\text{XY}}}+g\bar{V}(U,h)\end{eqnarray} \tag{ 6.10 }$$

where ${{H}_{\text{XY}}}$ is the XY Hamiltonian (4.1) and the nonlocal perturbation is given by

$$\begin{eqnarray}&&\bar{V}(U,h)=\frac{J}{L}\underset{s=0}{\overset{1}{\sum}}\,\underset{j=1}{\overset{2}{\sum}}\,{{U}^{j-1}}\left({{(-1)}^{sj}}\bar{H}_{s}^{z}\bar{H}_{s}^{z}+\bar{H}_{s,j}^{xy}\bar{H}_{s,j}^{yx}-\bar{H}_{s,j}^{xx}\bar{H}_{s,j}^{yy}\right)+h\bar{H}_{0}^{z}.\end{eqnarray} \tag{ 6.11 }$$

The time averaged quadratic operators appearing on the right hand side of (6.11) are the fundamental blocks of (4.6) and read as

$$\begin{eqnarray}\bar{H}_{s}^{z}&=&\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\sigma _{\ell}^{z}}=\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\text{i}a_{\ell}^{y}a_{\ell}^{x}}\\ \bar{H}_{s,j}^{xy}&=&\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\sigma _{\ell}^{x}{{\left(\sigma _{\ell +1}^{z}\right)}^{j-1}}\sigma _{\ell +j}^{y}}=\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\left(-\text{i}\right)a_{\ell}^{y}a_{\ell +j}^{y}}\\ \bar{H}_{s,j}^{yx}&=&\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\sigma _{\ell}^{y}{{\left(\sigma _{\ell +1}^{z}\right)}^{j-1}}\sigma _{\ell +j}^{x}}=\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\text{i}a_{\ell}^{x}a_{\ell +j}^{x}}\\ \bar{H}_{s,j}^{xx}&=&\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\sigma _{\ell}^{x}{{\left(\sigma _{\ell +1}^{z}\right)}^{j-1}}\sigma _{\ell +j}^{x}}=\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\left(-\text{i}\right)a_{\ell}^{y}a_{\ell +j}^{x}}\\ \bar{H}_{s,j}^{yy}&=&\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\sigma _{\ell}^{y}{{\left(\sigma _{\ell +1}^{z}\right)}^{j-1}}\sigma _{\ell +j}^{y}}=\overline{\frac{1}{2}\underset{\ell}{\sum}\,{{(-1)}^{s\ell}}\text{i}a_{\ell}^{x}a_{\ell +j}^{y}}.\end{eqnarray} \tag{ 6.12 }$$

Since we have to compute the time average of quadratic operators evolving according to a noninteracting Hamiltonian (see (3.5) and (6.4)), we can use (A.6). This allows to find the following exact result

$$\begin{eqnarray}\bar{\mathcal{O}}(k)&=&\underset{T\to +\infty}{{\lim}}\,\frac{1}{T}{\int}_{0}^{T}\text{d}t\mathcal{O}(k,t)=\underset{T\to +\infty}{{\lim}}\,\frac{1}{T}{\int}_{0}^{T}\text{d}t{{\text{e}}^{\text{i}{{\mathcal{H}}^{(2)}}(k)t}}\mathcal{O}(k){{\text{e}}^{-\text{i}{{\mathcal{H}}^{(2)}}(k)t}}\\&&\qquad=\frac{1}{2}\mathcal{O}(k,0)+\frac{1}{2}\left[{{\sigma}^{x}}{{\text{e}}^{\text{i}\frac{k}{2}{{\sigma}^{z}}}}\otimes {{\sigma}^{y}}{{\text{e}}^{\text{i}{{\theta}_{k}}{{\sigma}^{z}}}}\right]\,\mathcal{O}(k,0)\left[{{\sigma}^{x}}{{\text{e}}^{\text{i}\frac{k}{2}{{\sigma}^{z}}}}\otimes {{\sigma}^{y}}{{\text{e}}^{\text{i}{{\theta}_{k}}{{\sigma}^{z}}}}\right],\end{eqnarray} \tag{ 6.13 }$$

where $\mathcal{O}(k)$ is the symbol of a quadratic operator and $\overset{\overline{}}{{\cdots}}\,$ denotes the time average. The symbols of the operators (6.12) read

$$\begin{eqnarray}\bar{H}_{s}^{z}(k)&=&\frac{J}{\varepsilon _{k}^{2}}\left({{\delta}_{s,0}}{{\mathcal{Q}}_{2}}(k)-\gamma {{\delta}_{s,1}}{{\mathcal{Q}}_{8}}(k)\right) \\ \bar{H}_{s,1}^{xy}(k)&=&{{\delta}_{s,0}}{{\mathcal{Q}}_{4}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{6}}(k) \\ \bar{H}_{s,2}^{xy}(k)&=&{{\delta}_{s,0}}{{\mathcal{Q}}_{3}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{7}}(k) \\ \bar{H}_{s,1}^{yx}(k)&=&-{{\delta}_{s,0}}{{\mathcal{Q}}_{4}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{6}}(k) \\ \bar{H}_{s,2}^{yx}(k)&=&-{{\delta}_{s,0}}{{\mathcal{Q}}_{3}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{7}}(k) \\ \bar{H}_{s,1}^{xx}(k)&=&-\frac{J}{2\varepsilon _{k}^{2}}\left((1+\gamma )+(1-\gamma )\cos k\right)\left({{\delta}_{s,0}}{{\mathcal{Q}}_{1}}(k)-{{\delta}_{s,1}}{{\mathcal{Q}}_{5}}(k)\right) \\ \bar{H}_{s,2}^{xx}(k)&=&-\frac{J}{\varepsilon _{k}^{2}}\left[\gamma +(1-\gamma )(s+\cos k)\right]\left({{\delta}_{s,0}}{{\mathcal{Q}}_{2}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{8}}(k)\right) \\ \bar{H}_{s,1}^{yy}(k)&=&-\frac{J}{2\varepsilon _{k}^{2}}\left((1-\gamma )+(1+\gamma )\cos k\right)\left({{\delta}_{s,0}}{{\mathcal{Q}}_{1}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{5}}(k)\right) \\ \bar{H}_{s,2}^{yy}(k)&=&\frac{J}{\varepsilon _{k}^{2}}\left(\gamma +s(1-\gamma )-(1+\gamma ){{(-1)}^{s}}\cos k\right)\left({{\delta}_{s,0}}{{\mathcal{Q}}_{2}}(k)+{{\delta}_{s,1}}{{\mathcal{Q}}_{8}}(k)\right).\end{eqnarray} \tag{ 6.14 }$$

Here we expressed the results in terms of the symbols of the local charges of H_XY [65]

$$\begin{eqnarray}{{\mathcal{Q}}_{1}}(k)&=&\mathcal{I}_{1}^{+\left(\text{e}\right)}(k)={{\varepsilon}_{k}}\left[{{\sigma}^{x}}{{\text{e}}^{\text{i}\frac{k}{2}{{\sigma}^{z}}}}\right]\,\otimes \,\left[{{\sigma}^{y}}{{\text{e}}^{\text{i}{{\theta}_{k}}_{{}}{{\sigma}^{z}}}}\right] \\ {{\mathcal{Q}}_{2}}(k)&=&\mathcal{I}_{1}^{+\left(\text{o}\right)}(k)=\cos(k/2){{\varepsilon}_{k}}1\,\otimes \,\left[{{\sigma}^{y}}{{\text{e}}^{\text{i}{{\theta}_{k}}{{\sigma}^{z}}}}\right] \\ {{\mathcal{Q}}_{3}}(k)&=&\mathcal{I}_{1}^{-\left(\text{e}\right)}(k)=\sin(k)1\,\otimes \,1 \\ {{\mathcal{Q}}_{4}}(k)&=&\mathcal{I}_{1}^{-\left(\text{o}\right)}(k)=\sin(k/2)\left[{{\sigma}^{x}}{{\text{e}}^{\text{i}\frac{k}{2}{{\sigma}^{z}}}}\right]\,\otimes \,1 \\ {{\mathcal{Q}}_{5}}(k)&=&\mathcal{J}_{1}^{+\left(\text{e}\right)}(k)={{\varepsilon}_{k}}\left[{{\sigma}^{y}}{{\text{e}}^{\text{i}\frac{k}{2}{{\sigma}^{z}}}}\right]\,\otimes \,\left[{{\sigma}^{x}}{{\text{e}}^{\text{i}{{\theta}_{k}}{{\sigma}^{z}}}}\right] \\ {{\mathcal{Q}}_{6}}(k)&=&\mathcal{J}_{1}^{+\left(\text{o}\right)}(k)=\cos(k/2)\left[{{\sigma}^{y}}{{\text{e}}^{\text{i}\frac{k}{2}{{\sigma}^{z}}}}\right]\,\otimes \,{{\sigma}^{z}} \\ {{\mathcal{Q}}_{7}}(k)&=&\mathcal{J}_{1}^{-\left(\text{e}\right)}(k)=\sin(k){{\sigma}^{z}}\otimes \,{{\sigma}^{z}} \\ {{\mathcal{Q}}_{8}}(k)&=&\mathcal{J}_{1}^{-\left(\text{o}\right)}(k)=\sin(k/2){{\varepsilon}_{k}}{{\sigma}^{z}}\,\otimes \,\left[{{\sigma}^{x}}{{\text{e}}^{\text{i}{{\theta}_{k}}{{\sigma}^{z}}}}\right].\end{eqnarray} \tag{ 6.15 }$$

The first four symbols correspond to one-site shift invariant operators (the standard conservation laws of the quantum XY model), while the others change sign under a shift by one site.

We remind the reader that from the symbol of an operator it is possible to infer its locality properties [65]. In particular, a smooth symbol is associated with a quasi-local operator. If in addition the symbol has a finite number of nonzero Fourier components, as in (6.15), the associated operator is local. Equations (6.14) imply that $\bar{H}_{s}^{xy},\bar{H}_{s}^{yx}$ are local while $\bar{H}_{s}^{xx},\bar{H}_{s}^{yy},\bar{H}_{s}^{z}$ are quasi-local, thus the Hamiltonian (6.10) is a member of the quasi-local extension of the class $\mathcal{E}$ studied in section 5. As pointed out in section 5, we expect all the theorems of section 5, in particular corollary 5.4, to remain valid also for quasi-local operators. This guarantees the time evolution generated by (6.10) to be locally equivalent to the one generated by the following mean-field Hamiltonian

$$\begin{eqnarray}{{\bar{H}}_{\text{MF}}}(T)&=&{{H}_{XY}}+2Jg\underset{s}{\sum}\,\left({{(-1)}^{s}}+U\right)\frac{{{\langle \bar{H}_{s}^{z}\rangle}_{T}}}{L}\bar{H}_{s}^{z}+Jg\underset{s,j}{\sum}\,{{U}^{j-1}}\left(\frac{{{\langle \bar{H}_{sj}^{xy}\rangle}_{T}}}{L}\bar{H}_{sj}^{yx}+\frac{{{\langle \bar{H}_{sj}^{yx}\rangle}_{T}}}{L}\bar{H}_{sj}^{xy}\right)\\&&-Jg\underset{s,j}{\sum}\,{{U}^{j-1}}\left(\frac{{{\langle \bar{H}_{sj}^{xx}\rangle}_{T}}}{L}\bar{H}_{sj}^{yy}+\frac{{{\langle \bar{H}_{sj}^{yy}\rangle}_{T}}}{L}\bar{H}_{sj}^{xx}\right)+hg\bar{H}_{0}^{z},\end{eqnarray} \tag{ 6.16 }$$

where ${{\langle \mathcal{O}\rangle}_{T}}$ is the expectation value of the operator $\mathcal{O}$ in the mean-field description (see (5.15))

$$\begin{eqnarray}&&{{\langle \mathcal{O}\rangle}_{T}}=\text{Tr}\left[{{U}_{{\bar{V}}}}(T){{\rho}_{\text{GGE}}}U_{{\bar{V}}}^{\dagger}(T)\mathcal{O}\right].\end{eqnarray} \tag{ 6.17 }$$

To determine the time evolution generated by ${{\bar{H}}_{\text{MF}}}(T)$ we need to solve the self-consistency conditions encoded in (6.16) and (6.17). To this end, it is again convenient to exploit the representation in terms of symbols. Using (6.14), the symbol ${{\mathcal{H}}_{\text{MF}}}(k,T)$ of the time-dependent mean-field Hamiltonian can be written in terms of the symbols $\left\{{{\mathcal{Q}}_{\alpha}}(k),\alpha =1,\ldots ,8\right\}$, as follows⁷

$$\begin{eqnarray}&&{{\mathcal{H}}_{\text{MF}}}(k,T)=-{{\mathcal{Q}}_{1}}(k)+g{{\mathcal{V}}_{\text{MF}}}(k,T),\end{eqnarray} \tag{ 6.18 }$$

$$\begin{eqnarray}&&{{\mathcal{V}}_{\text{MF}}}(k,T)=\frac{h}{\varepsilon _{k}^{2}}{{\mathcal{Q}}_{2}}(k)+\underset{\alpha =1}{\overset{8}{\sum}}\,{{c}_{\alpha}}\left(k;{{\tilde{y}}_{\alpha}}\right){{\mathcal{Q}}_{\alpha}}(k).\end{eqnarray} \tag{ 6.19 }$$

The coefficients are given by

$$\begin{equation} \begin{array}{l} c_{1}(k; \tilde y_1)&=&-\frac{1+\cos k}{2\varepsilon^2_k}(\tilde y_{1}^{(0)}+\tilde y_{1}^{(1)})+\gamma^2\frac{1-\cos k}{2\varepsilon^2_k}(\tilde y_{1}^{(0)}-\tilde y_{1}^{(1)}),\nonumber\\ c_{2}(k; \tilde y_2)&=&2\frac{1+U}{ \varepsilon_k^2} \tilde y_{2}^{(0)}-2 U \frac{\cos k}{\varepsilon_k^2}\tilde y_{2}^{(1)} +2 U \gamma^2 \frac{1-\cos k}{\varepsilon_k^2}(\tilde y_{2}^{(0)}-\tilde y_{2}^{(1)}),\nonumber\\ c_{3}(k;\tilde y_{3})&=&-U(1+ \gamma^2) \tilde y_{3}^{(0)}-U(1-\gamma^2) \tilde y_{3}^{(1)},\nonumber\\ c_{4}(k;\tilde y_{4})&=&-(1+ \gamma^2) \tilde y_{4}^{(0)}-(1-\gamma^2) \tilde y_{4}^{(1)},\nonumber\\ c_{5}(k;\tilde y_5)&=& \frac{ 1+\cos k}{2\varepsilon^2_k} (\tilde y_{5}^{(0)}+\tilde y_{5}^{(1)})-\gamma^2 \frac{1-\cos k}{2\varepsilon^2_k} (\tilde y_{5}^{(0)}-\tilde y_{5}^{(1)}),\nonumber\\ c_{6}(k; \tilde y_{6})&=& (1+ \gamma^2) \tilde y_{6}^{(0)}+(1-\gamma^2) \tilde y_{6}^{(1)},\nonumber\\ c_{7}(k; \tilde y_{7})&=& U(1+ \gamma^2) \tilde y_{7}^{(0)}+U(1-\gamma^2) \tilde y_{7}^{(1)},\nonumber\\ c_{8}(k; \tilde y_{8})&=&2 \gamma^2\frac{U-1}{\varepsilon_k^2} \tilde y_{8}^{(0)}-2 \gamma^2 U \frac{\cos k}{ \varepsilon_k^2} \tilde y_{8}^{(1)}+2 U \frac{1+ \cos k}{\varepsilon_k^2} (\tilde y_{8}^{(0)}+\tilde y_{8}^{(1)}), \end{array} \end{equation} \tag{ 6.20 }$$

where we defined

$$\begin{eqnarray}&&\tilde{y}_{\alpha}^{(\ell )}(T)={\int}_{-\pi}^{\pi}\frac{\text{d}p}{2\pi}\frac{\cos(\ell p)}{\varepsilon _{k}^{2}}{{y}_{\alpha}}(p,T),\end{eqnarray} \tag{ 6.21 }$$

and

$$\begin{eqnarray}&&{{y}_{\alpha}}(k,T)=\frac{1}{8}\text{Tr}\left[{{U}_{{{\mathcal{H}}_{\text{MF}}}}}(k,T){{\Gamma}_{\text{GGE}}}(k)U_{{{\mathcal{H}}_{\text{MF}}}}^{\dagger}(k,T){{\mathcal{Q}}_{\alpha}}(k)\right] \\&&{{U}_{{{\mathcal{H}}_{\text{MF}}}}}(k,T)=\text{T}\exp\left[-\text{i}{\int}_{0}^{T}\text{d}s{{\mathcal{V}}_{\text{MF}}}(k,s)\right].\end{eqnarray} \tag{ 6.22 }$$

By taking the first derivative of (6.22) with respect to T and using the (closed) commutator algebra of ${{\mathcal{Q}}_{\alpha}}(k)$ we get

$$\begin{eqnarray}&&{{\dot{y}}_{\alpha}}(k,T)=\frac{h}{\varepsilon _{k}^{2}}\underset{\gamma =1}{\overset{8}{\sum}}\,f_{k}^{2\alpha \gamma}{{y}_{\gamma}}(k,T)+\underset{\beta ,\gamma =1}{\overset{8}{\sum}}\,{{c}_{\beta}}\left(k;{{\tilde{y}}_{\gamma}}\right)f_{k}^{\beta \alpha \gamma}{{y}_{\gamma}}(k,T).\end{eqnarray} \tag{ 6.23 }$$

The nonzero structure constants $f_{k}^{\alpha \beta \gamma}$ that are not connected to one another by symmetry are given by

$$\begin{eqnarray} &&f_{k}^{562}=f_{k}^{548}=2f_{k}^{647}=2 && \qquad f_{k}^{782}=f_{k}^{746}=2f_{k}^{845}=-2(1-\cos k)\\ && f_{k}^{584}=f_{k}^{526}=2f_{k}^{827}=-2\varepsilon _{k}^{2} && \qquad f_{k}^{728}=f_{k}^{764}=2f_{k}^{625}=2(1+\cos k)\end{eqnarray} \tag{ 6.24 }$$

The others follow from $f_{k}^{\beta \alpha \gamma}=-f_{k}^{\alpha \beta \gamma}$. In particular ${{\mathcal{Q}}_{1}}(k)$ and ${{\mathcal{Q}}_{3}}(k)$ commute with all the other charges, so y₁(k) and y₃(k) are conserved and the system (6.23) is reduced to 6 first order integro-differential equations that depend on a continuous variable k.

The solution of (6.23) entirely determines the time evolution generated by ${{\bar{H}}_{\text{MF}}}$. Indeed, the expectation value of any local observable in the pre-relaxation limit can be computed using the Wick theorem with the correlation matrix

$$\begin{eqnarray}&&\langle \left(\begin{array}{c} a_{2n-1}^{x} \\ \begin{array}{c} a_{2n-1}^{y} \\ \begin{array}{c} a_{2n-1}^{y} \\ a_{2n-1}^{y} \\ \end{array} \\ \end{array} \\ \end{array}\right)\left(a_{2\ell -1}^{x}a_{2\ell -1}^{y}a_{2\ell}^{x}a_{2\ell}^{y}\right)\,\rangle ={{\delta}_{\ell n}}{{1}_{4}}+{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}{{\text{e}}^{-\text{i}(n-\ell )k}}\underset{i=1}{\overset{8}{\sum}}\,\frac{8{{y}_{i}}(k,T)}{\text{tr}\left({{\mathcal{Q}}_{i}}{{(k)}^{2}}\right)}{{\mathcal{Q}}_{i}}(k).\end{eqnarray} \tag{ 6.25 }$$

This also means that the reduced density matrices of subsystems are gaussian at any time, so the two assumptions (a) and (b) could be also reformulated as a single hypothesis of RDMs being gaussian.

Equations (6.23) and (6.25) are the main results of this section: they allow us to compute the expectation values of local observables in the pre-relaxation limit of a weakly interacting model by solving a nonlinear system of differential equations, which is rather easy from a numerical point of view.

Reflection symmetry. The Hamiltonian (6.1) is reflection symmetric, that is to say it is invariant under the transformation

$$\begin{eqnarray}&&\sigma _{\ell}^{\alpha}\to \sigma _{s+L-\ell}^{\alpha}\quad \alpha \in \left\{x,y,z\right\},\end{eqnarray} \tag{ 6.26 }$$

where s is odd for reflections about a bond and even for those about a site.

The reflection operator acts on the Majorana fermions as follows

$$\begin{eqnarray}&&a_{\ell}^{x}\to \text{i}\left(\underset{j}{\prod}\,\sigma _{j}^{z}\right)a_{s+L-\ell}^{y} a_{\ell}^{y}\to -\text{i}\left(\underset{j}{\prod}\,\sigma _{j}^{z}\right)a_{s+L-\ell}^{x}.\end{eqnarray} \tag{ 6.27 }$$

Therefore the symbol $\mathcal{H}$ of a one-site shift invariant operator transforms as

$$\begin{eqnarray}&&{{\mathcal{H}}^{(1)}}(k)\to {{\sigma}^{y}}{{\mathcal{H}}^{(1)}}(-k){{\sigma}^{y}},\end{eqnarray} \tag{ 6.28 }$$

while for two-site shift invariant operators we find

$$\begin{equation}\mathcal H^{(2)}(k)\rightarrow\left\{ \begin{array}{ll} \sigma^x\otimes \sigma^y\ \mathcal H^{(2)}(-k)\ \sigma^x\otimes \sigma^y&s \rm{\ odd}\\ {\rm e}^{-{\rm i}\frac{k}{2}\sigma^z}\otimes \sigma^y\ \mathcal H^{(2)}(-k)\ {\rm e}^{{\rm i}\frac{k}{2}\sigma^z}\otimes \sigma^y&s \rm{\ even}. \end{array}\right. \end{equation} \tag{ 6.29 }$$

The symbols (6.15) of the conservation laws have the simple transformation rules

$$\begin{eqnarray}&&{{\mathcal{Q}}_{1,2}}(k){{\mathcal{Q}}_{1,2}}(k) \\&&{{\mathcal{Q}}_{3,4}}(k)-{{\mathcal{Q}}_{3,4}}(k) \\&&{{\mathcal{Q}}_{5,6}}(k)-{{(-1)}^{s}}{{Q}_{5,6}}(k) \\&&{{\mathcal{Q}}_{7,8}}(k){{(-1)}^{s}}{{\mathcal{Q}}_{7,8}}(k).\end{eqnarray} \tag{ 6.30 }$$

Since a shift by one site is equivalent to a reflection about a bond followed by a reflection about a site, we recover the transformation rules pointed out below (6.15).

If the initial state is reflection symmetric about a bond, ${{\mathcal{Q}}_{j}}(k)=0$ for j   =   3,4,7,8. Thus, the system of equations (6.23) can be reduced to

$$\begin{eqnarray}&&{{\dot{y}}_{2}}(k,T)=-2{{c}_{5}}\left(k;{{\tilde{y}}_{5}}\right)\varepsilon _{k}^{2}{{y}_{6}}(k,T)+{{c}_{6}}\left(k;{{\tilde{y}}_{6}}\right)(1+\cos k){{y}_{5}}(k,T) \\&&{{\dot{y}}_{5}}(k,T)=-2{{c}_{6}}\left(k;{{\tilde{y}}_{6}}\right){{y}_{2}}(k,T)+2\left(\frac{h}{\varepsilon _{k}^{2}}+{{c}_{2}}\left(k;{{\tilde{y}}_{2}}\right)\right)\varepsilon _{k}^{2}{{y}_{6}}(k,T) \\&&{{\dot{y}}_{6}}(k,T)=2{{c}_{5}}\left(k;{{\tilde{y}}_{5}}\right){{y}_{2}}(k,T)-\left(\frac{h}{\varepsilon _{k}^{2}}+{{c}_{2}}\left(k;{{\tilde{y}}_{2}}\right)\right)(1+\cos k){{y}_{5}}(k,T).\end{eqnarray} \tag{ 6.31 }$$

We numerically identified three different behaviours:

• Stationarity: The expectation values of the observables remain equal to the initial values given by the unperturbed GGE (figure 3).
• Local relaxation: The observables relax to a different stationary value: one-site shift invariance is restored in some cases (figure 4) while remaining broken in others (figure 5).
• Persistent oscillations: The amplitude of the oscillations of the expectation values of the observables does not approach zero (figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We point out that, even when there is relaxation (at some intermediate times with $Jt\gg {{g}^{-1}}$), the stationary state is not thermal, being the local conservation laws of ${{H}_{\text{XY}}}$ with symbol proportional to ${{\mathcal{Q}}_{1}}(k)$ and ${{\mathcal{Q}}_{3}}(k)$ (namely the charges that preserve non-abelian integrability) conserved in the pre-relaxation limit.

6.1. Perturbations preserving integrability

In this section we consider the case $U=\gamma =0$, in which H (6.1) is the Hamiltonian of the XXZ spin-$\frac{1}{2}$ chain. Because of the U(1) symmetry of rotations around z, there are many simplifications and the system of equations (6.31) can be rewritten as follows:

$$\begin{eqnarray}&&\frac{1}{2}\dot{y}_{2}^{\left[n\right]}(T)=y_{5}^{\left[n\right]}(T)y_{6}^{\left[0\right]}(T)-y_{6}^{\left[n\right]}(T)y_{5}^{\left[0\right]}(T) \\&&\frac{1}{2}\dot{y}_{5}^{\left[n\right]}(T)=-\left(2y_{2}^{\left[n\right]}(T)+y_{2}^{\left[n-1\right]}(T)+y_{2}^{\left[n+1\right]}(T)\right)y_{6}^{\left[0\right]}(T)+\left(h+4y_{2}^{\left[0\right]}(T)\right)y_{6}^{\left[n\right]}(T) \\&&\frac{1}{2}\dot{y}_{6}^{\left[n\right]}(T)=\left(2y_{2}^{\left[n\right]}(T)+y_{2}^{\left[n-1\right]}(T)+y_{2}^{\left[n+1\right]}(T)\right)y_{5}^{\left[0\right]}(T)-\left(h+4y_{2}^{\left[0\right]}(T)\right)y_{5}^{\left[n\right]}(T),\end{eqnarray} \tag{ 6.32 }$$

where we defined

$$\begin{eqnarray}&&\quad y_{2}^{\left[n\right]}(T)={\int}_{-\pi}^{\pi}\frac{\text{d}p}{2\pi}\frac{\cos(np)}{1+\cos p}{{y}_{2}}(p,T),\quad y_{5,6}^{\left[n\right]}(T)={\int}_{-\pi}^{\pi}\frac{\text{d}p}{2\pi}\cos(np){{y}_{5,6}}(p,T).\end{eqnarray} \tag{ 6.33 }$$

We notice that, despite the denominator, $y_{2}^{\left[n\right]}$ are expectation values of local operators, as well as $y_{5,6}^{\left[n\right]}$. Since ${{S}^{z}}=\frac{1}{2}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}$ commutes with the Hamiltonian, the dependence on h is simple and, in particular, the expectation value of the one-site shift invariant conservation laws is independent of the magnetic field. This means that the functions $y_{2}^{\left[n\right]}$ are independent of h.

It is useful to rewrite the system for n   =   0. We find

$$\begin{eqnarray}&&\dot{y}_{2}^{\left[0\right]}(T)=0 \\&&\dot{y}_{5}^{\left[0\right]}(T)=2\left(h+2y_{2}^{\left[0\right]}(T)-2y_{2}^{\left[1\right]}(T)\right)y_{6}^{\left[0\right]}(T) \\&&\dot{y}_{6}^{\left[0\right]}(T)=-2\left(h+2y_{2}^{\left[0\right]}(T)-2y_{2}^{\left[1\right]}(T)\right)y_{5}^{\left[0\right]}(T).\end{eqnarray} \tag{ 6.34 }$$

Inspecting the system we conclude that $y_{2}^{\left[0\right]}$ and ${{\left(y_{5}^{\left[0\right]}\right)}^{2}}+{{\left(y_{6}^{\left[0\right]}\right)}^{2}}$ are conserved. Moreover, the system (6.34) can be directly solved, it yields

$$\begin{eqnarray}&& y_{2}^{\left[0\right]}(T)=\frac{\langle {{S}^{z}}\rangle}{2L}\equiv \frac{{{s}^{z}}}{2} \\&&y_{5}^{\left[0\right]}(T)=y_{5}^{\left[0\right]}(0)\cos\left({\int}_{0}^{T}\text{d}\tau \left(2h+4m(\tau )\right)\right)+y_{6}^{\left[0\right]}(0)\sin\left({\int}_{0}^{T}\text{d}\tau \left(2h+4m(\tau )\right)\right) \\&&y_{6}^{\left[0\right]}(T)=y_{6}^{\left[0\right]}(0)\cos\left({\int}_{0}^{T}\text{d}\tau \left(2h+4m(\tau )\right)\right)-y_{5}^{\left[0\right]}(0)\sin\left({\int}_{0}^{T}\text{d}\tau \left(2h+4m(\tau )\right)\right).\end{eqnarray} \tag{ 6.35 }$$

Here we defined

$$\begin{eqnarray}&&\quad m(T)\equiv \frac{{{s}^{z}}}{2}-y_{2}^{\left[1\right]}(T)={{s}^{z}}-\frac{{{\langle {{Q}_{2}}\rangle}_{T}}}{L}=\frac{1}{4}\langle \sigma _{\ell}^{z}\rangle +\frac{1}{8}{{\langle \sigma _{\ell -1}^{x}\sigma _{\ell}^{z}\sigma _{\ell +1}^{x}+\sigma _{\ell -1}^{y}\sigma _{\ell}^{z}\sigma _{\ell +1}^{y}\rangle}_{T}}.\end{eqnarray} \tag{ 6.36 }$$

If both $y_{5}^{\left[0\right]}(T)$ and $y_{6}^{\left[0\right]}(T)$ are zero, (6.32) has the solution

$$\begin{eqnarray}&&y_{2}^{\left[n\right]}(T)=y_{2}^{\left[n\right]}(0) \\&&y_{5}^{\left[n\right]}(T)=y_{5}^{\left[n\right]}(0)\cos\left(2\left(h+2{{s}^{z}}\right)T\right)+y_{6}^{\left[n\right]}(0)\sin\left(2\left(h+2{{s}^{z}}\right)T\right) \\&&y_{6}^{\left[n\right]}(T)=y_{6}^{\left[n\right]}(0)\cos\left(2\left(h+2{{s}^{z}}\right)T\right)-y_{5}^{\left[n\right]}(0)\sin\left(2\left(h+2{{s}^{z}}\right)T\right).\end{eqnarray} \tag{ 6.37 }$$

For $h\ne -2{{s}^{z}}$, local observables keep oscillating in time, otherwise, on the pre-relaxation timescale, the expectation values of local observables do not move from the values reached at times $1\ll Jt\ll {{g}^{-1}}$.

More generally (${{\left(y_{5}^{\left[0\right]}\right)}^{2}}+{{\left(y_{6}^{\left[0\right]}\right)}^{2}}\ne 0$), from (6.35) and (6.36) we immediately infer that relaxation is possible only if

$$\begin{eqnarray}&&\exists \underset{T\to \infty}{{\lim}}\,m(T)=-\frac{h}{2},\end{eqnarray} \tag{ 6.38 }$$

$$\begin{eqnarray}&&\exists \underset{T\to \infty}{{\lim}}\,\mid {\int}_{0}^{T}\text{d}\tau \left(m(\tau )+\frac{h}{2}\right)\mid <\infty .\end{eqnarray} \tag{ 6.39 }$$

We see that m(T) could be interpreted as a sort of 'induced magnetisation' that h must compete with.

The trivial dependence on h is manifest choosing the variables

$$\begin{eqnarray}&&{{Y}_{n}}=y_{2}^{\left[n\right]}(T) \\&&{{\Phi}_{n}}=y_{5}^{\left[n\right]}(T)y_{5}^{\left[0\right]}(T)+y_{6}^{\left[n\right]}(T)y_{6}^{\left[0\right]}(T),\end{eqnarray} \tag{ 6.40 }$$

which satisfy the following system of equations independent of h:

$$\begin{eqnarray}&&m(T)\equiv \frac{{{s}^{z}}}{2}-{{Y}_{1}}(T)\\&& {{{\ddot {Y}}\,}_{n}}(T)=-4{{\Phi}_{0}}\left(2{{Y}_{n}}(T)+{{Y}_{n+1}}(T)+{{Y}_{n-1}}(T)\right)+8\left({{s}^{z}}-m(T)\right){{\Phi}_{n}}(T) \\&&{{{\dot \Phi}}_{n}}(T)=-2\left({{s}^{z}}-m(T)\right){{\dot{Y}}_{n}}(T).\end{eqnarray} \tag{ 6.41 }$$

Here we omitted the time dependence in the conserved quantity ${{\Phi}_{0}}$. Since ${{\Phi}_{0}}\ne 0$ by assumption, the original variables are obtained from the inverse transformation

$$\begin{eqnarray}&&y_{5}^{\left[n\right]}(T)=\frac{2y_{5}^{\left[0\right]}(T){{\Phi}_{n}}(T)+y_{6}^{\left[0\right]}(T){{{\dot{Y}}}_{n}}(T)}{2{{\Phi}_{0}}} \\&&y_{6}^{\left[n\right]}(T)=\frac{2y_{6}^{\left[0\right]}(T){{\Phi}_{n}}(T)-y_{5}^{\left[0\right]}(T){{{\dot{Y}}}_{n}}(T)}{2{{\Phi}_{0}}},\end{eqnarray} \tag{ 6.42 }$$

and (6.35). Performing a qualitative analysis of the system (6.41) we conclude that condition (6.38) and $m(T)+h/2$ approaching 0 faster than 1/T imply relaxation of local degrees of freedom. Therefore the variance (2.11) of m(T) is what in section 2 we called an 'order parameter' for the transition between relaxation and oscillatory behaviour.

As we will show in section 7, some aspects of the solutions of nonlinear systems like (6.41) can be worked out analytically. In the present context this would involve the study of quantum quenches from rather artificial initial states. Therefore we prefer to leave the entire discussion to section 7, where we will obtain a system of equations extremely similar to (6.41), with the advantage that the qualitative analysis can be carried out for more conventional initial states.

6.2. Perturbations breaking integrability: linearisation

In order to gain some insights into the time evolution under the Hamiltonian (6.1) in the non-integrable case we focus on quantum quenches starting from the dimer product state

$$\begin{eqnarray}&&\mid \text{MG}\rangle =\frac{\mid \uparrow \downarrow \rangle -\mid \downarrow \uparrow \rangle}{2}\otimes \cdots \otimes \frac{\mid \uparrow \downarrow \rangle -\mid \downarrow \uparrow \rangle}{2},\end{eqnarray} \tag{ 6.43 }$$

which is the ground state of the Majumdar–Ghosh Hamiltonian

$$\begin{eqnarray}&&{{H}_{0}}=\frac{J}{4}\underset{\ell =1}{\overset{L}{\sum}}\,{{\vec{\sigma}}_{\ell}}\cdot {{\vec{\sigma}}_{\ell +1}}+\frac{1}{2}{{\vec{\sigma}}_{\ell}}\cdot {{\vec{\sigma}}_{\ell +2}}.\end{eqnarray} \tag{ 6.44 }$$

Despite the model being interacting, (6.43) is a two-site shift invariant Slater determinant, whose correlation matrix has the following symbol

$$\begin{eqnarray}&&{{\Gamma}_{\text{MG}}}(k)={{\sigma}^{x}}\otimes {{\sigma}^{y}}.\end{eqnarray} \tag{ 6.45 }$$

The initial conditions for $\left\{{{y}_{\alpha}}(k)\right\}$ (6.22) are determined by the GGE correlation matrix for g   =   0. They can be obtained by expanding ${{\Gamma}_{\text{MG}}}(k)$ in the base of the symbols (6.15) of the conserved charges of H_XY (the remaining space is zeroed by the time evolution, as $1\ll Jt$, see section 5.1)

$$\begin{eqnarray}&&\Gamma(k;0)=\underset{i=1}{\overset{8}{\sum}}\,\frac{\text{tr}\left[{{\Gamma}_{\text{MG}}}(k){{\mathcal{Q}}_{i}}(k)\right]}{\text{tr}\left[{{\left({{\mathcal{Q}}_{i}}(k)\right)}^{2}}\right]}{{\mathcal{Q}}_{i}}(k).\end{eqnarray} \tag{ 6.46 }$$

We find

$$\begin{eqnarray}&&\quad \Gamma(k;0)=\frac{1+\cos k}{1+\cos k+{{\gamma}^{2}}(1-\cos k)}{{\mathcal{Q}}_{1}}(k)-\gamma \frac{1-\cos k}{1+\cos k+{{\gamma}^{2}}(1-\cos k)}{{\mathcal{Q}}_{5}}(k).\end{eqnarray} \tag{ 6.47 }$$

The only nonzero initial conditions are given by (see (6.22))

$$\begin{eqnarray}&&{{y}_{1}}(k,0)=\frac{1+\cos k}{4}, \\&&{{y}_{5}}(k,0)=-\gamma \frac{1-\cos k}{4}.\end{eqnarray} \tag{ 6.48 }$$

The initial state is reflection symmetric about a bond, so we can use the reduced system (6.31). Since y₅ is the only nonzero initial condition that appears in (6.31), for $\gamma =0$ (see previous section) the solution of the system of equations is independent of time, namely the pre-relaxation limit is trivial.

It is easy to see that also for h   =   0 system (6.31) and (6.48) has a stationary solution. We therefore assume $\gamma ,h\ne 0$. Since c_j(k;y_j) are linear homogeneous functions of y_j, the magnetic field h enters into the equations essentially as a scale factor. We rescale the variables as follows

$$\begin{eqnarray}&&\tau =2hT=2hgt\quad \epsilon =\frac{\gamma}{2h}\quad {{z}_{j}}=\frac{2{{y}_{j}}}{\gamma}\quad {{\gamma}_{j}}=\frac{2{{c}_{j}}}{\gamma}.\end{eqnarray} \tag{ 6.49 }$$

From (6.31) we then obtain

$$\begin{eqnarray}&&{{\partial}_{\tau}}{{z}_{2}}(k,\tau )=-\epsilon \varepsilon _{k}^{2}{{\gamma}_{5}}(k,\tau ){{z}_{6}}(k,\tau )+\epsilon {{\cos}^{2}}\frac{k}{2}{{\gamma}_{6}}(k,\tau ){{z}_{5}}(k,\tau ) \\&&{{\partial}_{\tau}}{{z}_{5}}(k,\tau )=-\epsilon \,{{\gamma}_{6}}(k,\tau ){{z}_{2}}(k,\tau )+\left(1+\epsilon \varepsilon _{k}^{2}{{\gamma}_{2}}(k,\tau )\right){{z}_{6}}(k,\tau ) \\&&{{\partial}_{\tau}}{{z}_{6}}(k,\tau )=\epsilon \,{{\gamma}_{5}}(k,\tau ){{z}_{2}}(k,\tau )-\left(1+\epsilon \varepsilon _{k}^{2}{{\gamma}_{2}}(k,\tau )\right)\frac{{{\cos}^{2}}\frac{k}{2}}{\varepsilon _{k}^{2}}{{z}_{5}}(k,\tau ),\end{eqnarray} \tag{ 6.50 }$$

with the initial conditions

$$\begin{eqnarray}&&{{z}_{2}}(k,0)={{z}_{6}}(k,0)=0\quad {{z}_{5}}(k,0)=-{{\sin}^{2}}\frac{k}{2}.\end{eqnarray} \tag{ 6.51 }$$

For generic $\epsilon$ the system of equations is not exactly solvable, but the limit of small $\epsilon$ allows a linear approximation. For not too large rescaled times (we will come back to this point later) the terms that are multiplied by $\epsilon$ in the last two equations can be neglected, while the functions that appear on the right hand side of the first equation can be computed at $O\left({{\epsilon}^{0}}\right)$.

For z₅ and z₆ we obtain the simple solution

$$\begin{eqnarray}&&{{z}_{5}}(k,\tau )\approx -{{\sin}^{2}}\frac{k}{2}\cos\left(\frac{\cos\frac{k}{2}}{{{\varepsilon}_{k}}}\tau \right)\\&& {{z}_{6}}(k,\tau )\approx {{\sin}^{2}}\frac{k}{2}\frac{\cos\frac{k}{2}}{{{\varepsilon}_{k}}}\sin\left(\frac{\cos\frac{k}{2}}{{{\varepsilon}_{k}}}\tau \right),\end{eqnarray} \tag{ 6.52 }$$

while z₂ is a slightly more complicated function that involves integrals over the momentum of z_5,6, namely

$$\begin{eqnarray}&&{{z}_{2}}(k,\tau )\approx \epsilon \frac{\sin k \sin\frac{k}{2}}{2}{\int}_{-\pi}^{\pi}\frac{\text{d}p}{2\pi}\underset{\sigma =\pm}{\sum}\,{{g}_{\sigma}}(k,p)f\left(\frac{\cos\frac{k}{2}}{{{\varepsilon}_{k}}}+\sigma \frac{\cos\frac{p}{2}}{{{\varepsilon}_{p}}};\tau \right).\end{eqnarray} \tag{ 6.53 }$$

Here we defined

$$\begin{eqnarray}&&\,\,\,{{g}_{\sigma}}(k,p)\equiv \frac{{{\sin}^{2}}\frac{p}{2}}{{{\varepsilon}_{p}}}\left[\frac{{{\cos}^{2}}\frac{k}{2}{{\cos}^{2}}\frac{p}{2}-{{\gamma}^{2}}{{\sin}^{2}}\frac{k}{2}{{\sin}^{2}}\frac{p}{2}}{{{\varepsilon}_{k}}{{\varepsilon}_{p}}}+\sigma \cos\frac{k}{2}\cos\frac{p}{2}\right],\end{eqnarray} \tag{ 6.54 }$$

$$\begin{eqnarray}&&\quad \,f\left(x;\tau \right)\equiv \frac{1-\cos(x\tau )}{x}.\end{eqnarray} \tag{ 6.55 }$$

For small $\epsilon$ and given $\gamma$, ${{z}_{2}}(k,\tau )$ relaxes to the stationary value

$$\begin{eqnarray}&&{{z}_{2}}(k,\infty )=\frac{\epsilon {{\sin}^{2}}k}{8{{\gamma}^{2}}}\left[1-3{{\gamma}^{2}}+\left(2+4{{\gamma}^{2}}\right)\cos k+\left(1-{{\gamma}^{2}}\right)\cos 2k\right].\end{eqnarray} \tag{ 6.56 }$$

Let us now estimate the time window in which the linear approximation is applicable. From (6.52) it follows that ${{\gamma}_{5}}$ and ${{\gamma}_{6}}$ decay to zero as ${{\tau}^{-\frac{3}{2}}}$. Instead, since z₂ approaches a nonzero stationary value, ${{\gamma}_{2}}$ is of the same order of z₂. This means that, as the time increases, the first term on the right hand side of the last two equations of (6.50) becomes more and more negligible with respect to the other term multiplied by $\epsilon$. By neglecting the former we obtain essentially the same solution (6.52) as before, with the replacement

$$\begin{eqnarray}&&\tau \to \tau +\epsilon \varepsilon _{k}^{2}{\int}_{0}^{\tau}\text{d}s{{\gamma}_{2}}(k,s)=\tau \left(1+\epsilon \varepsilon _{k}^{2}{{\gamma}_{2}}(k,\infty )\right)+\ldots .\end{eqnarray} \tag{ 6.57 }$$

Being ${{\gamma}_{2}}\sim O(\epsilon )$, after a rescaled time $\tau \sim \frac{1}{{{\epsilon}^{2}}}$, the correction to z_5,6 becomes comparable with the function itself. Assuming that the relevant part of the time evolution occurs within this time scale, the linear approximation is justified only if $\mid {{z}_{2}}\mid \ll 1$ (and $\epsilon \ll 1$). For $\gamma <1/2$ we find $\mid {{z}_{2}}(k,\infty )\mid <\frac{\epsilon}{6{{\gamma}^{2}}}$, so we obtain the consistency condition

$$\begin{eqnarray}&&\frac{1}{12h}\ll \gamma \ll 2h.\end{eqnarray} \tag{ 6.58 }$$

Figures 7 and 8 report a comparison between the solution of the linearised problem and the full numerical solution of system (6.31) for a set of parameters fulfilling (6.58).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the expressions (6.52) of ${{z}_{5}}(k,\tau ),{{z}_{6}}(k,\tau )$ and (6.56) of ${{z}_{2}}(k,\tau )$, we can directly compute the time evolution of the expectation value of any local observable in the pre-relaxation limit. Indeed corollary 5.4 allows us to apply the Wick theorem at any time (in the limit under examination) and the correlation matrix is given by (6.25).

Any integral involving z₅ and z₆ approaches zero and z₂ becomes independent of time even if not integrated. Therefore, in the limit (6.58) and for large times the expectation value of any local observable relaxes to a stationary value that can be described by the correlation matrix with symbol

$$\begin{eqnarray}&&\underset{\tau \to {{\infty}^{\left(*\right)}}}{{\lim}}\,\Gamma\left(k;\tau \right)=\frac{\gamma {{z}_{1}}(k,\infty )}{\varepsilon _{k}^{2}}{{\mathcal{Q}}_{1}}(k)+\frac{\gamma {{z}_{2}}(k,\infty )}{\varepsilon _{k}^{2}{{\cos}^{2}}(k/2)}{{\mathcal{Q}}_{2}}(k),\end{eqnarray} \tag{ 6.59 }$$

where the infinite time limit $\tau \to {{\infty}^{\left(*\right)}}$ must be understood within the limits of validity of the linear approximation.

We point out that one-site shift invariance is restored, indeed the only contributions to the correlation matrix at infinite times arise from ${{\mathcal{Q}}_{1}}(k)$ and ${{\mathcal{Q}}_{2}}(k)$, which are symbols of one-site shift invariant operators. The manifestly one-site shift invariant expression of the correlation matrix in the limit (6.59) reads

$$\begin{eqnarray}&&\underset{\tau \to {{\infty}^{\left(*\right)}}}{{\lim}}\,\langle \left(\begin{array}{c} a_{n}^{x} \\ a_{n}^{y} \\ \end{array}\right)\left(\begin{array}{c} a_{\ell}^{x} & a_{\ell}^{y} \\ \end{array}\right)\rangle ={{\delta}_{\ell n}}{{1}_{2}}+{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}{{\text{e}}^{-\text{i}(n-\ell )k}}{{\Gamma}^{(1)}}(k),\end{eqnarray} \tag{ 6.60 }$$

with

$$\begin{eqnarray}{{\Gamma}^{(1)}}(k)&=&\frac{2h(1+\cos 2k)+{{\sin}^{2}}k\cos k\left(1-3{{\gamma}^{2}}+\left(2+4{{\gamma}^{2}}\right)\cos 2k+\left(1-{{\gamma}^{2}}\right)\cos 4k\right)}{2h\left(1+\cos 2k+{{\gamma}^{2}}(1-\cos 2k)\right)}\\&&\times \left(\cos k\text{}\text{}{{\sigma}^{y}}-\gamma \sin k\text{}\text{}{{\sigma}^{x}}\right).\end{eqnarray} \tag{ 6.61 }$$

6.3. Remarks on the late time dynamics

We notice that, despite one-site shift invariance being restored in (6.59), the asymptotic value is not given by the average over a shift of the expectation value of the operator in the GGE of the unperturbed model. Indeed the one-site shift average of (6.47) is proportional to ${{\mathcal{Q}}_{1}}(k)$ (see (6.30)) but the symbol of the large time correlation matrix (6.59) has also a term proportional to ${{\mathcal{Q}}_{2}}(k)$. Consequently, the shift-averaged stationary values can not be recovered from those in the limit of small perturbation $g\to 0$. For example we have

$$\begin{eqnarray}&&\underset{1\ll t\ll \frac{1}{g}}{{\lim}}\,\langle \frac{\sigma _{2\ell -1}^{z}+\sigma _{2\ell}^{z}}{2}\rangle =O(g)\end{eqnarray} \tag{ 6.62 }$$

$$\begin{eqnarray}&&\underset{2hgt\to {{\infty}^{\left(*\right)}}}{{\lim}}\,\langle \sigma _{\ell}^{z}\rangle =-\frac{{{\gamma}^{2}}}{16h}\frac{3+\gamma}{{{(1+\gamma )}^{3}}}+O(g),\end{eqnarray} \tag{ 6.63 }$$

where we highlighted that there are O(g) corrections. Besides this particular quench in a non-integrable model, similar issues arise also in the integrable case, where it is generally believed that at infinite time after the quench the expectation values can be computed in a GGE constructed with the (quasi-)local conservation laws of the model. In the rest of the section we propose a reasoning that relates this kind of deviations to possible atypical properties of the model.

We notice that at times $1\ll t\ll \frac{1}{g}$, shift-symmetrised expectation values of quadratic operators, e.g. (6.62), can be obtained from the GGE of the unperturbed model constructed with only the local translation invariant conservation laws. This is because the limit $1\ll t\ll \frac{1}{g}$ with $g\to 0$ can be described by the GGE of the unperturbed model, but the conservation laws that break translation invariance are odd under a shift by one site, so they can in fact be neglected (this equivalence breaks down for operators that consist of the product of more than two Jordan–Wigner fermions).

For nonzero g non-abelian integrability is supposed to break down and the relevant charges are generally assumed to be one-site shift invariant and in involution with one another. We now speculate about the stationary state in the limit of small g if the perturbation does not break integrability. Let the (quasi-)local conservation laws of the interacting (integrable) model be in a smooth one-to-one correspondence with the local one-site shift invariant conservation laws for g   =   0. In the limit of small g, the stationary state should locally approach the GGE constructed with the local one-site shift invariant conservation laws of the unperturbed model. The expectation values of shift-symmetrised quadratic operators at times $1\ll t\ll \frac{1}{g}$ are compatible with such a one-site shift invariant GGE. However, discrepancies like that between (6.62) and (6.63) show that at larger times there is a time window in which the expectation values approach a different value. Our assumption of regularity of the conservation laws as a function of g requires that at even larger times the expectation values should eventually relax to the same values they had in the earliest plateau. This is clearly possible, but an infinite number of operators displaying a similar behaviour is rather surprising. This makes us wonder whether the hypothesis of regularity could break down, that is to say there are (quasi-)local conservation laws for $g\ne 0$ that become nonlocal when g   =   0 (e.g. their typical range could be singular as $g\to 0$) or, vice versa, some one-site shift invariant conservation laws of the unperturbed Hamiltonian do not have analogues at nonzero g.

From this point of view, discrepancies like that between (6.62) and (6.63) could be indications that in the XYZ model there might be quasi-local conservation laws that do not behave well in the limit $g\to 0$. This scenario becomes even more plausible if one takes into account the issues in the construction of the GGE in XXZ spin-$\frac{1}{2}$ chains that were recently pointed out [42–44].

6.4. Summary

The results of section 4 are a compelling motivation for the study of nonlocal Hamiltonians of the form (4.9). We now go beyond that rigid derivation: we skip the formal steps that relate (4.9) to a pre-relaxation limit and start off directly with a Hamiltonian of the form (4.9). We then query whether such models with (apparently) non-integrable long-range interactions could display thermal-like behaviour at late times after a quench. Specifically, we consider the Hamiltonian

$$\begin{eqnarray}&&H\left(\tilde{g},\lambda \right)=-\underset{\ell}{\overset{L}{\sum}}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\tilde{g}\sigma _{\ell}^{z}\right)+\frac{\lambda}{L}{{\left({\sum}_{\ell}^{L}\sigma _{\ell}^{z}\right)}^{2}}.\end{eqnarray} \tag{ 7.1 }$$

This has been recently proposed as a convenient model to investigate pre-thermalisation issues [60]. In fact, in order to recover some temporal cluster decomposition properties, the authors of [60] considered a slightly modified version of (7.1), where, in the term proportional to $\lambda$, $\underset{\ell}{\overset{L}{\sum}}\,\sigma _{\ell}^{z}$ was replaced by $\underset{\ell}{\overset{L}{\sum}}\,\sigma _{\ell}^{z}-\overline{\underset{\ell}{\overset{L}{\sum}}\,\sigma _{\ell}^{z}}$, the latter being its time average for $\lambda =0$. The time average is a simple quadratic quasi-local operator (see the first equation of (6.14)), so our formalism could be readily applied. However, for the sake of simplicity, we consider (4.9) and show later that a redefinition of $\tilde{g}$ is sufficient to recover the results shown in [60].

From corollary 5.4, in the thermodynamic limit $L\to \infty$ the time evolution under (7.1) is locally equivalent to the time evolution under the time-dependent mean-field Hamiltonian

$$\begin{eqnarray}&&H_{\text{MF}}^{{{\Psi}_{0}}}(t)=-\underset{\ell}{\sum}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+h(t)\sigma _{\ell}^{z}\right),\end{eqnarray} \tag{ 7.2 }$$

with h(t) the solution of the self-consistent equation

$$\begin{equation} h(t)=\tilde g-2\lambda \langle{\Psi_0\vert {\rm T^\dagger}\exp\left({\rm i}\int_0^t{\rm d}\tau H_{\rm MF}^{\Psi_0}(\tau)\right)\sigma_\ell^z\ {\rm T}\exp\left(-{\rm i}\int_0^t{\rm d}\tau H_{\rm MF}^{\Psi_0}(\tau)\right)\vert \Psi_0}\rangle. \end{equation} \tag{ 7.3 }$$

Here we used translation invariance to replace $\frac{1}{L}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}$ by the local operator $\sigma _{\ell}^{z}$ (which removes the nasty dependence on the chain length L). For Slater determinant initial states the expectation value can be conveniently written in terms of the symbol $\Gamma(k)$ of the initial correlation matrix (see appendix A). In particular we find

$$\begin{eqnarray}&&h(t)=\tilde{g}-\lambda {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}{{y}_{k}}(t),\end{eqnarray} \tag{ 7.4 }$$

with

$$\begin{eqnarray}&&{{y}_{k}}(t)=\text{tr}\left[{{U}_{{{\mathcal{H}}_{\text{MF}}}}}(k,t)\Gamma(k)U_{{{\mathcal{H}}_{\text{MF}}}}^{\dagger}(k,t){{\sigma}^{y}}\right],\end{eqnarray} \tag{ 7.5 }$$

$$\begin{eqnarray}&&{{U}_{{{\mathcal{H}}_{\text{MF}}}}}(k,t)=\text{T}\exp\left[2\text{i}{\int}_{0}^{t}\text{d}\tau {{\sigma}^{y}}\left(h(\tau )-{{\text{e}}^{\text{i}k{{\sigma}^{z}}}}\right)\right].\end{eqnarray} \tag{ 7.6 }$$

One can easily verify that y_k is the solution of

$$\begin{eqnarray}&&y_{k}^{\prime \prime}(t)=4\left(h(t)-\cos k\right){{\phi}_{k}}(t)+16\left(h(t)\cos k-1\right){{y}_{k}}(t) \\&&\phi _{k}^{\prime}(t)=-4h(t)y_{k}^{\prime}(t),\end{eqnarray} \tag{ 7.7 }$$

with

$$\begin{eqnarray}&&{{\phi}_{k}}(t)=-4\text{tr}\left[{{U}_{{{\mathcal{H}}_{\text{MF}}}}}(k,t)\Gamma(k)U_{{{\mathcal{H}}_{\text{MF}}}}^{\dagger}(k,t){{\sigma}^{y}}{{\text{e}}^{\text{i}k{{\sigma}^{z}}}}\right].\end{eqnarray} \tag{ 7.8 }$$

Equations (7.4) and (7.7) could be solved by discretising the momenta k; however working in real space is more transparent. We therefore introduce the (real space) Fourier coefficients

$$\begin{equation} \tilde y_n=\int_{-\pi}^\pi\frac{{\rm d} k}{2\pi}\cos(n k)y_k(t)\qquad \tilde\phi_n=\int_{-\pi}^\pi\frac{{\rm d} k}{2\pi}\cos(n k)\phi_k(t), \end{equation} \tag{ 7.9 }$$

which are the expectation values of local operators with range n and n + 1 respectively. The system of equations (7.7) can then be recast as follows

$$\begin{equation} \begin{array}{l} h&=&g-\lambda \tilde y_0\\ \tilde y_0^{\prime\prime}&=&4 h \tilde\phi_0-4\tilde\phi_{1}+16 h\tilde y_{1}-16\tilde y_{0}\\ \tilde y_n^{\prime\prime}&=&4 h \tilde\phi_n-2(\tilde\phi_{n+1}+\tilde\phi_{n-1})+8 h (\tilde y_{n+1}+\tilde y_{n-1})-16\tilde y_{n}& (n>0)\\ \tilde \phi_n^{\prime}&=&-4h \tilde y_n^{\prime}& (n\geqslant 0). \end{array} \end{equation} \tag{ 7.10 }$$

The similarity with the system of equations (6.41) for the pre-relaxation limit in XXZ spin chains is remarkable, although the meaning of the variables is different.

For sufficiently large n (in order to be outside of the (deformed) light cone), both ${{\tilde{y}}_{n}}$ and $\tilde\phi_n$ are small (in the noncritical case they are exponentially small in n). Thus, the error originated from truncating the system of equations to $n\leqslant N$ decays very fast to zero with N and (7.10) can be conveniently reduced to a finite system of differential equations. This is the regularisation that we used in our numerical investigations⁸.

The system of equations has at least one integral of motion, namely the energy per unit length $\varepsilon =\langle {{\Psi}_{0}}\mid H\mid {{\Psi}_{0}}\rangle /L$. This can be written as follows

$$\begin{equation} \varepsilon=\frac{h^2-\tilde g^2}{4\lambda}-\frac{1}{8}\tilde\phi_0. \end{equation} \tag{ 7.11 }$$

In addition, (7.1) has infinite noninteracting conservation laws that are odd under reflection symmetry [19, 64]. The Dzyaloshinskii–Moriya interaction

$$\begin{eqnarray}&&{{H}_{\text{D}-\text{M}}}=\underset{\ell}{\sum}\,\sigma _{\ell}^{x}\sigma _{\ell +1}^{y}-\sigma _{\ell}^{y}\sigma _{\ell +1}^{x}\end{eqnarray} \tag{ 7.12 }$$

is one of them. For generic initial states, this is sufficient to rule out thermalisation. However, we embrace the point of view of [60] and wonder whether at infinite time after a quench from a reflection symmetric state some form of thermalisation arises.

7.1. To relax or not to relax

At late times the mean-field Hamiltonian (7.2) can result in two distinct behaviours: either the dynamics is governed by a (time independent) TFIC Hamiltonian (viz. $\exists {{\lim}_{t\to \infty}}h(t)$), or there is no relaxation (viz. $\nexists {{\lim}_{t\to \infty}}h(t)$).

In the following we will focus on quantum quenches from the ground state of the TFIC Hamiltonian

$$\begin{eqnarray}&&{{H}_{0}}=-\underset{\ell}{\sum}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+{{g}_{0}}\sigma _{\ell}^{z}\right).\end{eqnarray} \tag{ 7.13 }$$

A numerical analysis suggests that for generic values of g₀, $\tilde{g}$ and $\lambda$, the system of equation (7.7) does not always describe a relaxation process. We indeed found cases in which h(t) oscillates with almost constant amplitude (see figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If there are (finite) regions of the parameter space associated with relaxation and regions that are instead characterised by persistent oscillations, some quantities shall behave non-analytically at the boundaries of the regions. For example, the relaxation parameter

$$\begin{eqnarray}&&\Delta {{h}_{T}}=\sqrt{\frac{1}{T}{\int}_{T}^{2T}\text{d}\tau {{h}^{2}}(\tau )-{{\left(\frac{1}{T}{\int}_{T}^{2T}\text{d}\tau h(\tau )\right)}^{2}}}\end{eqnarray} \tag{ 7.14 }$$

approaches zero as $T\to \infty$ in the regions of (local) relaxation. On the other hand, if there are (sufficiently regular) persistent oscillations, $\Delta {{h}_{T}}$ remains nonzero for arbitrarily large T.

We numerically analysed the region of the parameter space in which there is no relaxation (at least apparently). This generally happens for sufficiently large $\lambda$ (see figures 2 and 9). In the vicinity of (a numerical estimation of) ${{\lambda}_{c}}$, $\Delta {{h}_{T}}$ is nicely fitted by a line (see figure 10)

$$\begin{eqnarray}&&\Delta {{h}_{T}}=\alpha \left(\lambda -{{\lambda}_{c}}\right)+O\left({{\left(\lambda -{{\lambda}_{c}}\right)}^{2}}\right)\quad \lambda >{{\lambda}_{c}}.\end{eqnarray} \tag{ 7.15 }$$

Since $\Delta {{h}_{T}}$ is positive, linear behaviour is indicative of discontinuous derivative at $\lambda ={{\lambda}_{c}}$. Indeed, for $\lambda <{{\lambda}_{c}}$, $\Delta {{h}_{T}}$ is compatible with zero (see e.g. the three solutions with $\lambda <{{\lambda}_{c}}$ in figure 10).

7.1.1. Small quench.

In the case g₀   =   h(0) the solution of (7.10) is independent of time (as a consequence of corollary 5.5). This can be turned into a condition on the parameter $\tilde{g}$ of the Hamiltonian as follows

$$\begin{eqnarray}&&\tilde{g}=\bar{g}\left({{g}_{0}},\lambda \right),\end{eqnarray} \tag{ 7.16 }$$

where we defined

$$\begin{eqnarray}&&\bar{g}\left({{g}_{0}},\lambda \right)={{g}_{0}}+2\lambda {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{{{g}_{0}}-\cos k}{\sqrt{1+g_{0}^{2}-2{{g}_{0}}\cos k}}.\end{eqnarray} \tag{ 7.17 }$$

The initial state $\mid {{\psi}_{{{g}_{0}}}}\rangle$ corresponding to this quench is then an effective eigenstate of $H\left(\tilde{g},\lambda \right)$ (7.1). In addition, in appendix [a:GS]E we show that one of the solutions of (7.16) corresponds to the state $\mid {{\psi}_{{{g}_{0}}}}\rangle$ with minimal energy among Slater determinants. As a matter of fact, the numerical analysis indicates that $\mid {{\psi}_{{{g}_{0}}}}\rangle$ is the true ground state of (7.1). We can therefore use $\mid {{\psi}_{{{g}_{0}}}}\rangle$ as a reference state to define the limit of small quench.

In this limit, the transition relaxation/no-relaxation can be understood more clearly. Indeed, choosing the parameters such that (7.16) is approximately satisfied, both $\mid {{\tilde{y}}_{n}}(t)-{{\tilde{y}}_{n}}(0)\mid$ and $\vert \tilde \phi_n(t)-\tilde \phi_n(0)\vert$ turn out to be small at any time. We can therefore linearise the system of equations (7.10) isolating a time independent contribution from ${{\tilde{y}}_{n}}$ and ${{\tilde{\phi}}_{n}}$:

$$\begin{equation} \begin{array}{l} h=g_0+(\tilde g- \bar g(g_0,\lambda)- \lambda\delta y_0)\\ \tilde y_n=\bar y_n+\delta y_n\\ \tilde \phi_n=\bar \phi_n+\delta \phi_n,\\ \end{array} \end{equation} \tag{ 7.18 }$$

where variables with a bar on top are expectation values calculated in $\mid {{\psi}_{{{g}_{0}}}}\rangle$. We then obtain

$$\begin{eqnarray}&&\delta y_{n}^{\prime \prime}\approx -16\left(1+g_{0}^{2}\right)\delta {{y}_{n}}+16{{g}_{0}}\left(\delta {{y}_{n+1}}+\delta {{y}_{n-1}}\right)-\lambda {{\bar{v}}_{n}}\delta {{y}_{0}}+\left(\tilde{g}-\bar{g}\right){{\bar{v}}_{n}},\end{eqnarray} \tag{ 7.19 }$$

where

$$\begin{eqnarray}&&{{\bar{v}}_{n}}=4\left({{\bar{\phi}}_{n}}+2{{\bar{y}}_{n-1}}+2{{\bar{y}}_{n+1}}\right).\end{eqnarray} \tag{ 7.20 }$$

The system of equations (7.19) can be readily solved. Since for quenches from the ground states of TFIC Hamiltonians $\tilde{y}_{n}^{\prime}(0)=0$, we find

$$\begin{equation} \delta y_n(t)\approx (\tilde g-\bar g)\sum_j k_j \left[1-\cos(\sqrt{a_j}t)\right] w_{n, j} \end{equation} \tag{ 7.21 }$$

where a_j and w_{n, j} are the eigenvalues and the components of the (right) eigenvectors (at fixed j) of

$$\begin{eqnarray}&&{{A}_{\ell n}}=16\left(1+g_{0}^{2}\right){{\delta}_{\ell n}}-16{{g}_{0}}{{\delta}_{\mid \ell -n\mid ,1}}-16{{g}_{0}}{{\delta}_{\ell 0}}{{\delta}_{n1}}+\lambda {{\bar{v}}_{\ell}}{{\delta}_{n0}}\end{eqnarray} \tag{ 7.22 }$$

and k_j are given by

$$\begin{eqnarray}&&\vec{k}={{(AW)}^{-1}}\overrightarrow{{\bar{v}}}\quad \left({{\left[\vec{k}\right]}_{j}}={{k}_{j}},\quad {{\left[W\right]}_{nj}}={{w}_{n,j}},\quad {{\left[\overrightarrow{{\bar{v}}}\right]}_{j}}={{\bar{v}}_{j}}\right).\end{eqnarray} \tag{ 7.23 }$$

For $\lambda =0$, A can be diagonalised in momentum space; the eigenvalues are given by

$$\begin{eqnarray}&&{{a}_{j}}=16\left(1+g_{0}^{2}-2{{g}_{0}}\cos{{k}_{j}}\right),\end{eqnarray} \tag{ 7.24 }$$

and, in the limit $N\to \infty$ (N is our regularisation parameter), the spectrum becomes continuous (the density of k_j is uniform in $(0,\pi )$) and the eigenvectors unbounded. It is not difficult to show that the rank-1 perturbation $\lambda {{\bar{v}}_{\ell}}{{\delta}_{n0}}$ does not change the continuous part of the spectrum, which is still described by (7.24) and, as $N\to \infty$, the density of k_j remains uniform.

For sufficiently small $\lambda$ the spectrum is continuous. Given that k_j w_n,j is a smooth function of a_j for a fix n (as we numerically checked), we can apply the Riemann–Lebesgue lemma to extract the large time behaviour of (7.21). We find that it relaxes to

$$\begin{eqnarray}&&{{\tilde{y}}_{n}}(\infty )\approx {{\bar{y}}_{n}}+\left(\tilde{g}-\bar{g}\right)\underset{j}{\sum}\,{{\left[{{A}^{-1}}\right]}_{nj}}{{\bar{v}}_{j}}.\end{eqnarray} \tag{ 7.25 }$$

The power-law corrections to this result can be obtained by stationary phase approximations.

On the other hand, when $\pm \lambda$ exceeds some critical value $\pm \lambda _{c}^{\pm}$, A's spectrum develops an isolated eigenvalue and the corresponding eigenvector is bounded, whereas the continuous part of the spectrum is still described by (7.24). The main consequence is that the oscillations associated with the isolated eigenvalue do not cancel by dephasing mechanisms and local degrees of freedom keep oscillating at arbitrarily large times.

7.1.2. The bound state.

The isolated eigenvalue a₀ can be easily worked out in momentum space

$$\begin{eqnarray}&&\underset{n}{\sum}\,{{A}_{\ell n}}{{w}_{n,0}}={{a}_{0}}{{w}_{\ell ,0}}\Rightarrow \left[{{a}_{0}}-16\left(1+g_{0}^{2}-2{{g}_{0}}\cos k\right)\right]w(k;0)=\lambda \bar{v}(k){{w}_{0,0}},\end{eqnarray} \tag{ 7.26 }$$

with $w(k;j)={{w}_{0,j}}+2\underset{n>0}{\sum}\,\cos(nk){{w}_{n,j}}$ and $\bar{v}(k)={{\bar{v}}_{0}}+2\underset{n>0}{\sum}\,\cos(nk){{\bar{v}}_{n}}$. Isolated eigenvalues are such that $\omega =\sqrt{{{a}_{0}}}/4$ is outside of the continuous band, given by the image of $\sqrt{1+g_{0}^{2}-2{{g}_{0}}\cos k}$. We are therefore allowed to write

$$\begin{eqnarray}&&\frac{w(k;0)}{{{w}_{0,0}}}=\frac{\lambda}{16}\frac{\bar{v}(k)}{{{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k},\end{eqnarray} \tag{ 7.27 }$$

which, integrated over k, gives the condition

$$\begin{eqnarray}&&1=\frac{\lambda}{16}{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{\bar{v}(k)}{{{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k}.\end{eqnarray} \tag{ 7.28 }$$

The critical values $\lambda _{c}^{\pm}$ at the boundaries of the relaxation region correspond to the extrema of the continuous band where $\omega =1+\text{sgn}(\lambda )\mid {{g}_{0}}\mid$. Using the explicit form of ${{\bar{v}}_{n}}$ we then find

$$\begin{eqnarray}&&\lambda _{c}^{\pm}={{\left({\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{{{\sin}^{2}}k}{\left(\text{sgn}(\lambda )\mid {{g}_{0}}\mid \,+{{g}_{0}}\cos k\right)\sqrt{1+g_{0}^{2}-2{{g}_{0}}\cos k}}\right)}^{-1}}\end{eqnarray} \tag{ 7.29 }$$

and for $\lambda >\lambda _{c}^{+}$ or $\lambda <\lambda _{c}^{-}$ equation (7.28) can be rewritten as

$$\begin{eqnarray}&&{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{{{\sin}^{2}}k}{\left({{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k\right)\sqrt{1+g_{0}^{2}-2{{g}_{0}}\cos k}}=\frac{1}{2\lambda}.\end{eqnarray} \tag{ 7.30 }$$

Taking the derivative of this expression with respect to $\omega$ we obtain

$$\begin{eqnarray}&&\frac{\text{d}\lambda}{\text{d}\omega}>0\quad 0\leqslant \omega <\mid 1-\mid {{g}_{0}}\mid \mid \vee \omega >1+\mid {{g}_{0}}\mid ;\end{eqnarray} \tag{ 7.31 }$$

this fact, together with the observation that $\lambda >0$ if $\omega >1\,+\,\mid {{g}_{0}}\mid$ and $\lambda <0$ if $0\leqslant \omega <\mid 1\,-\,\mid {{g}_{0}}\mid \mid$, implies that $\lambda$ is an increasing function of $\omega$ (in the allowed dominion), therefore it is injective. This means that there can not be more than one isolated eigenvalue for a fixed $\lambda$.

The right eigenvector corresponding to the isolated eigenvalue is given by

$$\begin{eqnarray}&&{{w}_{n,0}}\propto {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{{{\sin}^{2}}k\cos(nk)}{\left({{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k\right)\sqrt{1+g_{0}^{2}-2{{g}_{0}}\cos k}},\end{eqnarray} \tag{ 7.32 }$$

which decays exponentially with n, confirming that it is a bound state. Analogously, the left eigenvector $w_{0,n}^{L}$ reads as (the factor in front of the integral is due to A's asymmetry)

$$\begin{eqnarray}&&w_{0,n}^{L}\propto \left(2-{{\delta}_{n0}}\right){\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\frac{\cos(nk)}{{{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k}=\frac{\left(2-{{\delta}_{n0}}\right){{(-1)}^{n}}{{\text{e}}^{-n\theta}}}{\sqrt{{{\left(1+g_{0}^{2}-{{\omega}^{2}}\right)}^{2}}-4g_{0}^{2}}}\quad \left(\lambda >0\right),\end{eqnarray} \tag{ 7.33 }$$

where

$$\begin{eqnarray}&&\theta =\text{arccosh}\left(\frac{{{\omega}^{2}}-1-g_{0}^{2}}{2{{g}_{0}}}\right).\end{eqnarray} \tag{ 7.34 }$$

Using (7.21) and (7.23) we can compute the entire contribution of the isolated eigenvalue to the solution of (7.19):

$$\begin{equation} (\tilde g-\bar g)\frac{\vec {\bar v}\cdot \vec w_0^L\ \vec w_0}{\vec w_0^L\cdot \vec w_0}\frac{1-\cos(4\omega t)}{16\omega^2}. \end{equation} \tag{ 7.35 }$$

We stress that, within the linear approximation, the oscillation frequency is independent of $\tilde{g}$. Figure 11 shows that the most important correction to (7.35) lies precisely in the frequency, essentially because it is multiplied by the time, which has to be large for the subleading (time dependent) contributions to be negligible. However, in not-too-large time windows, the numerical data are in excellent agreement with (7.35) (the expression must be modified including a corrective phase shift if $\left(\tilde{g}-\bar{g}\right)Jt$ is not small).

Zoom In Zoom Out Reset image size

We point out that there is a third relevant point ${{\lambda}_{*}}<\lambda _{c}^{-}\left(<0\right)$ at which $\omega =0$:

$$\begin{eqnarray}&&{{\lambda}_{*}}=-{{\left({\int}_{-\pi}^{\pi}\frac{\text{d}k}{\pi}\frac{{{\sin}^{2}}k}{{{\left(1+g_{0}^{2}-2{{g}_{0}}\cos k\right)}^{3/2}}}\right)}^{-1}}.\end{eqnarray} \tag{ 7.36 }$$

For $\lambda <{{\lambda}_{*}}$, the isolated eigenvalue of A becomes negative (i.e. $\omega$ becomes purely imaginary). This would result in an exponential growth of (7.21), which after some finite time would no longer be consistent with the linearisation procedure. However, we point out that, for small quenches in which the energy is close to the ground state one, $\lambda$ is always larger than ${{\lambda}_{*}}$. Indeed, if we assume $\lambda \leqslant {{\lambda}_{*}}\left({{g}_{0}}\right)$, it turns out that the ground state of (7.13) is not equivalent to the ground state of $H\left(\bar{g}\left({{g}_{0}},\lambda \right),\lambda \right)$, with $\bar{g}\left({{g}_{0}},\lambda \right)$ given by (7.17). The latter state is instead associated with the ground state of the TFIC Hamiltonian (7.13) with magnetic field $g_{0}^{\prime}\ne {{g}_{0}}$ such that $\bar{g}\left(g_{0}^{\prime},\lambda \right)=\bar{g}\left({{g}_{0}},\lambda \right)$ and ${{\lambda}_{*}}\left(g_{0}^{\prime}\right)<\lambda$. In the quench dephasing diagram of figure 2 we can indeed easily identify $\lambda _{c}^{\pm}$ (the 'critical' curves for $\lambda$ positive and negative, respectively) but there is no trace of ${{\lambda}_{*}}$.

Interpretation. The bound state of the matrix A may be put in relation with the existence of localised excitations in (7.1). We emphasise that this is not an ab initio calculation but rather a physical picture that explains the observations.

Since the time evolution of the expectation value of any local observable in $\mid {{\psi}_{{{g}_{0}}}}\rangle$ is stationary, we can assume that, to all intents and purposes, $\mid {{\psi}_{{{g}_{0}}}}\rangle$ is an eigenstate of $H\left(\bar{g}\left({{g}_{0}},\lambda \right),\lambda \right)$. We now consider the limit of small quench and assume $\mid \lambda \mid >\mid {{\lambda}_{c}}\mid$. From (7.21) it follows that the projection on the bound state of $y_{n}^{\prime \prime}(t)+4\text{i}\omega y_{n}^{\prime}(t)$ is proportional to an oscillating phase

$$\begin{eqnarray}&&\underset{n}{\sum}\,w_{0,n}^{L}\left(y_{n}^{\prime \prime}(t)+4\text{i}\omega y_{n}^{\prime}(t)\right)\propto {{\text{e}}^{4\text{i}\omega t}}.\end{eqnarray} \tag{ 7.37 }$$

The left hand side can be written as the expectation value $\langle {{\psi}_{{{{\tilde{g}}}_{0}}}}(t)\mid B_{{{g}_{0}}}^{\dagger}\mid {{\psi}_{{{{\tilde{g}}}_{0}}}}(t)\rangle$ of a noninteracting operator with symbol

$$\begin{eqnarray}&&\sin k\frac{\left(h(t)-\cos(k)\right){{\sigma}^{x}}-\sin k{{\sigma}^{y}}+\text{i}\omega {{\sigma}^{z}}}{{{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k},\end{eqnarray} \tag{ 7.38 }$$

where ${{\tilde{g}}_{0}}$ approximately satisfyes (7.16). In the no-quench limit the mean-field parameter is constant $h(t)\to {{g}_{0}}$ (and ${{\tilde{g}}_{0}}\to {{g}_{0}}$), so such operator can be written as

$$\begin{eqnarray}B_{{{g}_{0}}}^{\dagger}\sim \underset{\ell ,n}{\sum}\,\left(\begin{array}{c} a_{\ell}^{x} & a_{\ell}^{y} \\ \end{array}\right){{\left[\mathcal{B}_{{{g}_{0}}}^{\dagger}\right]}_{\ell -n}}\left(\begin{array}{c} a_{n}^{x} \\ a_{n}^{y} \\ \end{array}\right),\end{eqnarray} \tag{ 7.39 }$$

where

$$\begin{eqnarray}&&{{\left[\mathcal{B}_{{{g}_{0}}}^{\dagger}\right]}_{\ell}}={\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}{{\text{e}}^{\text{i}\ell k}}\sin k\frac{\left({{g}_{0}}-\cos(k)\right){{\sigma}^{x}}-\sin k{{\sigma}^{y}}+\text{i}\omega {{\sigma}^{z}}}{{{\omega}^{2}}-1-g_{0}^{2}+2{{g}_{0}}\cos k}\end{eqnarray} \tag{ 7.40 }$$

and $a_{\ell}^{\alpha}$ are the Majorana fermions (4.3). In (7.39) we left out the normalisation. Importantly, ${{\left[\mathcal{B}_{{{g}_{0}}}^{\dagger}\right]}_{\ell}}\sim {{\text{e}}^{-\mid \ell \mid \theta}}$, with $\theta$ defined in (7.34); we have therefore found a quasi-local operator whose expectation value approximately oscillates in time as in (7.37). As a consequence, $\mathcal{B}_{{{g}_{0}}}^{\dagger}$ acts like an excitation over the initial state. Indeed we have

$$\begin{eqnarray}&&\langle {{\psi}_{{{{\tilde{g}}}_{0}}}}(t)\mid B_{{{g}_{0}}}^{\dagger}\mid {{\psi}_{{{{\tilde{g}}}_{0}}}}(t)\rangle =\underset{\ell ,n}{\sum}\,\langle {{\psi}_{{{{\tilde{g}}}_{0}}}}\mid \ell \rangle \langle n\mid {{\psi}_{{{{\tilde{g}}}_{0}}}}\rangle {{\text{e}}^{\text{i}\left({{E}_{\ell}}-{{E}_{n}}\right)}}{{\left[B_{{{g}_{0}}}^{\dagger}\right]}_{\ell n}}\sim {{\text{e}}^{4\text{i}\omega t}}\end{eqnarray} \tag{ 7.41 }$$

and if the state $\mid {{\psi}_{{{{\tilde{g}}}_{0}}}}\rangle$ has a sufficiently general representation (i.e. the overlaps with the eigenstates of $H\left(\tilde{g},\lambda \right)$ are generally nonzero), (7.41) tells us that $B_{{{g}_{0}}}^{\dagger}$ connects only states with energy difference equal to $4\omega$.

In conclusion, the bound state of (7.22) seems to be a manifestation of a localised excitation of $H\left(\tilde{g},\lambda \right)$.

7.1.3. Remark on quenches from the ordered phase.

In the previous section we ignored a subtlety that in principle could have invalidated part of the discussion (and part of the diagram in figure 2). The mapping to a mean-field Hamiltonian relies on cluster decomposition properties but the ground state of the TFIC Hamiltonian in the ferromagnetic phase is the superposition of two Slater determinants [21], that separately do not possess cluster decomposition properties. Nevertheless, the discussion (and in turn figure 2) remains correct also in this problematic case, at least for the operators commuting with $\underset{\ell}{\prod}\,\sigma _{\ell}^{z}$. This can be seen as follows:

• (i)  
  The mean-field mapping is exact for the true ground state, which breaks the spin-flip symmetry realised by $\underset{\ell}{\prod}\,\sigma _{\ell}^{z}$. Thus (7.3) is valid.
• (ii)  
  Using that ${{\sigma}^{z}}$ is a quadratic operator in the Jordan–Wigner fermions, in (7.3) the ground state can be replaced by one of the two Slater determinants, therefore h(t) is still solution of (7.7). Analogously, the expectation value of any operator commuting with $\underset{\ell}{\prod}\,\sigma _{\ell}^{z}$ can be found replacing the initial state with one of the two Slater determinants and then using Wick theorem.
• (iii)  
  For any given t, the expectation value of operators ${{\mathcal{O}}^{o}}(\ell )$ that anticommute with $\underset{\ell}{\prod}\,\sigma _{\ell}^{z}$, like the order parameter, can be obtained from the large r limit of $\langle {{\Psi}_{0}}\mid U_{\text{MF}}^{\dagger}(t){{\mathcal{O}}^{o}}(\ell ){{\mathcal{O}}^{o}}(\ell +r){{U}_{\text{MF}}}(t)\mid {{\Psi}_{0}}\rangle$, using cluster decomposition properties.

However, using similar general arguments we are not able to exclude that the expectation values of the odd operators might keep oscillating also when all the even operators relax. There is indeed a subtle problem of limits that comes from the trick of Point (iii). Nevertheless, in the cases considered we have never encountered this situation, suggesting that for the model under examination such complications do not arise. Here we provide a heuristic argument. When the limit of infinite time for h(t) exists, at sufficiently long times the dynamics is essentially the same as for a quantum quench in the TFIC from a certain state with cluster decomposition properties. In the latter situation we can apply (a direct generalisation of) the results of [19], which showed that the expectation values of odd operators decay to zero. This suggests that the diagram of figure 2 could be valid also taking into account the odd operators.

7.2. Relaxation properties

We now focus on the regions (of the parameter space) in which the limit ${{\lim}_{t\to \infty}}h(t)\equiv {{h}_{\infty}}$ exists and investigate more closely the relaxation properties. Even beyond the linear approximation (7.21), we can still guess an asymptotic form for the solution y_k of (7.7)

$$\begin{eqnarray}&&{{y}_{k}}\to {{\zeta}_{k}}+\left({{A}_{k}}{{\text{e}}^{2\text{i}{{\varepsilon}_{k}}t}}+\text{h}.\text{c}.\right),\end{eqnarray} \tag{ 7.42 }$$

where ${{\varepsilon}_{k}}\sim 2\sqrt{1+h_{\infty}^{2}-2{{h}_{\infty}}\cos k}$. Equation (7.42) is compatible with the relaxation of h(t) with corrections O(t^−1/2−j) (stationary phase approximation), where j is an integer that depends on the behaviour of A_k around the extremal points of the dispersion relation ($k=0\vee \pi$). Our numerical analysis for several quenches from TFIC initial states is compatible with j   =   1. This is not surprising since the same exponent governs the late time behaviour of $\langle {{\sigma}^{z}}\rangle$ after quenches in the TFIC [21]. Thus, we conjecture the large time expansion

$$\begin{eqnarray}&&h(t)\sim {{h}_{\infty}}+\frac{{{A}_{0}}\cos\left(4\left(1-{{h}_{\infty}}\right)t+{{\varphi}_{0}}\right)+{{A}_{\pi}}\cos\left(4\left(1+{{h}_{\infty}}\right)t+{{\varphi}_{\pi}}\right)}{{{t}^{3/2}}},\end{eqnarray} \tag{ 7.43 }$$

and similar behaviours for ${{\tilde{y}}_{n}}(t)$. Remarkably, the (leading) oscillatory frequency is only determined by ${{h}_{\infty}}$. In figure 12, (7.43) is compared against numerical data for a quench that leads to relaxation.

Zoom In Zoom Out Reset image size

Since for asymptotically large times the time evolution is equivalent to the one generated by the TFIC Hamiltonian

$$\begin{eqnarray}&&{{H}_{f}}=-\underset{\ell}{\sum}\,\left(\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+{{h}_{\infty}}\sigma _{\ell}^{z}\right),\end{eqnarray} \tag{ 7.44 }$$

at late times local observables are described by a generalised Gibbs ensemble constructed with the local conservation laws of H_f. However, the Lagrange multipliers can not be simply fixed by computing the corresponding integrals of motion at the initial time, as they are in fact conserved only at asymptotically large times. This is an example of a stationary state written in terms of operators commuting with one another but not with the Hamiltonian. In fact, this is not the first time that such an unusual description emerges: in [61] the pre-thermalisation plateau was described by a GGE constructed with operators in involution that are however not conserved at the perturbative order that was worked out.

As it will be clarified in the next section, such a stationary state coincides with the pre-thermalisation plateau of [60]. In the thermodynamic limit this is just the stationary state that emerges at infinite time after the quench. Indeed, in the regions in which there is relaxation, we have not found any indication of pre-relaxation/pre-thermalisation behaviour. Furthermore, 'relaxation' is not synonym of 'thermalisation'. Indeed it is not difficult to show that at late times the system still retains infinite information about the initial state. To this aim, as initial state we choose the ground state of the local Hamiltonian

$$\begin{eqnarray}&&H\left[\left\{a\right\}\right]=\underset{j}{\overset{n}{\sum}}\,{{a}_{j}}{{H}_{j}},\quad {{a}_{j}}\in \mathbb{R},\end{eqnarray} \tag{ 7.45 }$$

where n is finite and H_j are the most local reflection symmetric conservation laws of a TFIC model with Hamiltonian H₀. Such state can be easily constructed [82] and, using the notations of appendix E, is completely characterised by the function

$$\begin{eqnarray}&&m(k)=-\text{sgn}\left(\underset{j}{\overset{n}{\sum}}\,{{a}_{j}}\cos(jk)\right),\end{eqnarray} \tag{ 7.46 }$$

which is equal to 1 if and only if the excitation $\alpha _{k}^{\dagger}$ of H₀ is present in the state. It is important to note that different characteristic functions m(k) correspond to locally inequivalent states. Thermal-like behaviour would imply that the only information about m(k) that is retained at late times is the corresponding energy and magnetisation. We now consider the special cases in which the initial magnetisation is such that H₀ is the mean-field Hamiltonian at the initial time. In this way, by corollary 5.5, the expectation value of local observables is independent of time and the late time stationary state is equivalent to the initial state. The only scenario compatible with thermal-like behaviour is that each distinct function m(k) of the form (7.46) corresponds to a distinct pair $\left\{{{h}_{m}},{{\varepsilon}_{m}}\right\}$, where h_m is the parameter of the late time mean-field Hamiltonian and ${{\varepsilon}_{m}}$ the energy. The self-consistent conditions behind the no-quench limit are worked out in appendix E and are given by

$$\begin{eqnarray}&&{{h}_{m}}=g+2\lambda {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}m(k)\frac{{{h}_{m}}-\cos k}{\sqrt{1+h_{m}^{2}-2{{h}_{m}}\cos k}}.\end{eqnarray} \tag{ 7.47 }$$

$$\begin{eqnarray}&&{{\varepsilon}_{m}}=\frac{h_{m}^{2}-{{g}^{2}}}{4\lambda}-{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}m(k)\frac{{{h}_{m}}\cos k-1}{\sqrt{1+h_{m}^{2}-2{{h}_{m}}\cos k}}.\end{eqnarray} \tag{ 7.48 }$$

It is easy to see that the same pair $\left\{{{h}_{m}},{{\epsilon}_{m}}\right\}$ is associated with infinitely many functions m(k) of the form (7.46) (it is enough to choose n   =   3 in (7.46) to find (infinite) examples). Thus, the non-equilibrium time evolution under (7.1) does not generally result in thermalisation.

We conclude the analysis of (7.1) considering the expectation value ${{\langle {{Q}_{n}}\rangle}_{t}}$ of the local conservation laws of H_f (7.44), which can be written as follows [19]

$$\begin{eqnarray}&&\frac{{{Q}_{n}}}{L}={\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\cos(nk){{\varepsilon}_{k}}\left(\alpha _{k}^{\dagger}\alpha _{k}^{\dagger}-\frac{1}{2}\right),\end{eqnarray} \tag{ 7.49 }$$

where ${{\varepsilon}_{k}}$ is the dispersion relation of H_f and ${{\alpha}_{k}}$ are noninteracting fermions that diagonalise H_f. Using free-fermion techniques we obtain

$$\begin{equation} \langle{\frac{Q_n}{L}}\rangle_t=-\frac{h_\infty}{2}\tilde y_n-\frac{1}{8}\tilde\phi_n. \end{equation} \tag{ 7.50 }$$

In particular, the expectation value of ${{H}_{f}}\equiv {{Q}_{0}}$ can be written as follows

$$\begin{eqnarray}&&{{\langle \frac{{{H}_{f}}}{L}\rangle}_{t}}=\varepsilon -\frac{{{\left(h(t)-{{h}_{\infty}}\right)}^{2}}-{{\left({{h}_{\infty}}-\tilde{g}\right)}^{2}}}{4\lambda},\end{eqnarray} \tag{ 7.51 }$$

from which it is clear that the relaxation exponent of ${{\langle {{H}_{f}}\rangle}_{t}}$ is twice the exponent of h (which in the cases that we investigated is 3/2, see (7.43)). This result can be easily generalised to any local conservation law of H_f. By taking the time derivative of (7.50), the last equation of (7.10) implies

$$\begin{eqnarray}&&\langle \frac{{{Q}_{n}}}{L}\rangle _{t}^{\prime}=\frac{h-{{h}_{\infty}}}{2}y_{n}^{\prime}.\end{eqnarray} \tag{ 7.52 }$$

Since both $h-{{h}_{\infty}}$ and $y_{n}^{\prime}$ are expected to decay as t^−3/2 (with oscillatory factors like in (7.43)), we immediately obtain

$$\begin{eqnarray}&&{{\langle \frac{{{Q}_{n}}}{L}\rangle}_{\infty}}-{{\langle \frac{{{Q}_{n}}}{L}\rangle}_{t}}\sim O\left({{t}^{-3}}\right).\end{eqnarray} \tag{ 7.53 }$$

7.2.1. Comparison with [60].

It is not a coincidence that the same relaxation exponents were found in [60] in a perturbative framework. In order to understand the relation between the two models we must come back to the modified version of Hamiltonian considered in [60], i.e. (2.8). The mean-field Hamiltonian for that precise model reads

$$\begin{equation} H_{\rm MF}(t)=-\sum_\ell\sigma_\ell^x\sigma_{\ell+1}^x-\texttt{g}\sum_\ell\sigma_\ell^z+2\lambda\langle{\sigma_\ell^z-\bar\sigma_\ell^z}\rangle_{t, \rm MF}\left(\sum_\ell\sigma_\ell^z-\overline{\sum_\ell\sigma_\ell^z}\right), \end{equation} \tag{ 7.54 }$$

where the time average $\overset{\overline{}}{{\cdots}}\,$ is taken with respect to the Hamiltonian with $\lambda =0$. It is convenient to introduce the auxiliary Hamiltonian

$$\begin{equation} \tilde{H}_{\rm MF}(t)=-\sum_\ell\sigma_\ell^x\sigma_{\ell+1}^x-\texttt{g}\sum_\ell\sigma_\ell^z+2\lambda\langle{\sigma_\ell^z-\tilde\sigma_\ell^z}\rangle_{t, \rm MF}\left(\sum_\ell\sigma_\ell^z-\widetilde{\sum_\ell\sigma_\ell^z}\right), \end{equation} \tag{ 7.55 }$$

where the time average $\widetilde{\cdots}$ is now taken with respect to $\tilde{H}_{f}\equiv\lim_{t\rightarrow\infty}\tilde{H}_{\rm MF}(t)$, limit which is assumed to exist. It is now simple to prove that $\tilde H_f=H(\text{g},0)$, where $H\left(\text{g},\lambda \right)$ is given in equation (2.8). To this aim let us consider the expectation value of ${{\sigma}^{z}}$ evolving via ${{\tilde{H}}_{\text{MF}}}(t)$; it fulfils

$$\begin{eqnarray}\underset{t\to \infty}{{\lim}}\,{{\langle \tilde{\sigma}_{\ell}^{z}\rangle}_{t,\text{MF}}}&=&\underset{t\to \infty}{{\lim}}\,\underset{T\to \infty}{{\lim}}\,\frac{1}{T}{\int}_{0}^{T}\text{d}s\langle {{\Psi}_{0}}\mid \tilde{U}_{\text{MF}}^{\dagger}(t){{\text{e}}^{\text{i}{{{\tilde{H}}}_{f}}s}}\sigma _{\ell}^{z}{{\text{e}}^{-\text{i}{{{\tilde{H}}}_{f}}s}}{{\tilde{U}}_{\text{MF}}}(t)\mid {{\Psi}_{0}}\rangle \\&=&\underset{t\to \infty}{{\lim}}\,\underset{T\to \infty}{{\lim}}\,\frac{1}{T}{\int}_{0}^{T}\text{d}s\langle {{\Psi}_{0}}\mid \tilde{U}_{\text{MF}}^{\dagger}(t+s)\sigma _{\ell}^{z}{{\tilde{U}}_{\text{MF}}}(t+s)\mid {{\Psi}_{0}}\rangle \\&=&\underset{t\to \infty}{{\lim}}\,{{\langle \sigma _{\ell}^{z}\rangle}_{t,\text{MF}}}={{\langle {{\sigma}^{z}}\rangle}_{\infty ,\text{MF}}},\end{eqnarray} \tag{ 7.56 }$$

where ${{\tilde{U}}_{\text{MF}}}(t)$ is the time evolution operator constructed with ${{\tilde{H}}_{\text{MF}}}(t)$ and in the second step we replaced ${{\text{e}}^{-\text{i}{{{\tilde{H}}}_{f}}s}}{{\tilde{U}}_{\text{MF}}}(t)$ with ${{\tilde{U}}_{\text{MF}}}(t+s)$, as it is legitimate at late times. From this it follows $\tilde H_f=H(\text{g},0)$ and, in turn, the equivalence between (7.54) and (7.55) ${{H}_{\text{MF}}}(t)={{\tilde{H}}_{\text{MF}}}(t)$. Importantly, this means that (7.54) and (7.13) have the same infinite time limit if we set ${{h}_{\infty}}=\text{g}$.

As a matter of fact, the equivalence between (7.54) and (7.13) is not restricted to infinite times. The mean-field time evolution operator for (7.54) can indeed be written as follows

$$\begin{eqnarray}&&{{U}_{\text{MF}}}(t)={{\text{e}}^{4\text{i}\lambda {\int}_{0}^{t}\text{d}s{{\langle \delta {{\sigma}^{z}}\rangle}_{s,\text{MF}}}\overline{{{S}^{z}}}}}\text{T}\exp\left(-\text{i}{\int}_{0}^{t}\text{d}s{{H}_{f}}+4\lambda {{\langle \delta {{\sigma}^{z}}\rangle}_{s,\text{MF}}}{{S}^{z}}(s)\right).\end{eqnarray} \tag{ 7.57 }$$

Here ${{S}^{z}}=\frac{1}{2}\underset{\ell}{\sum}\,\sigma _{\ell}^{z}$, $\delta {{\sigma}^{z}}={{\sigma}^{z}}-{{\bar{\sigma}}^{z}}$ and

$$\begin{eqnarray}&&{{S}^{z}}(s)={{\text{e}}^{-4\text{i}\lambda {\int}_{0}^{s}\text{d}{{s}^{\prime}}{{\langle \delta {{\sigma}^{z}}\rangle}_{{{s}^{\prime}},\text{MF}}}\overline{{{S}^{z}}}}}{{S}^{z}}{{\text{e}}^{4\text{i}\lambda {\int}_{0}^{s}\text{d}{{s}^{\prime}}{{\langle \delta {{\sigma}^{z}}\rangle}_{{{s}^{\prime}},\text{MF}}}\overline{{{S}^{z}}}}}.\end{eqnarray} \tag{ 7.58 }$$

If the magnetisation relaxes faster than $1/{{t}^{1+\alpha}}$, with $\alpha >0$, the operator at the exponent of the first term of (7.57) is a bounded function of the time, so that exponential can be safely expanded. In fact the entire term can be neglected (it gives corrections $o(\lambda )$). The same holds true in S^z(s), indeed the finiteness of $\underset{s}{{\sup}}\,\lambda {\int}_{0}^{s}\text{d}{{s}^{\prime}}{{\langle \delta {{\sigma}^{z}}\rangle}_{{{s}^{\prime}},\text{MF}}}$ guarantees that the series expansion of the exponentials in S^z(s) can be truncated for any s with an error that goes to zero as $\lambda \to 0$. By considering the first terms of the expansion one immediately realizes that the correction is $o(\lambda )$ and approaches zero for large s as $1/{{s}^{\alpha}}$. Putting everything together, replacing S^z(s) by S^z in (7.57) introduces an error $o(\lambda )$, independently of the time. Therefore we obtain

$$\begin{eqnarray}&&{{U}_{\text{MF}}}(t)\sim \text{T}\exp\left(-\text{i}{\int}_{0}^{t}\text{d}s{{H}_{f}}+4\lambda {{\langle \delta {{\sigma}^{z}}\rangle}_{s,\text{MF}}}{{S}^{z}}\right)+o(\lambda ).\end{eqnarray} \tag{ 7.59 }$$

We can also replace ${{\langle \delta {{\sigma}^{z}}\rangle}_{s,\text{MF}}}$ with ${{\langle {{\sigma}^{z}}\rangle}_{s,\text{MF}}}-{{\langle {{\sigma}^{z}}\rangle}_{\infty ,\text{MF}}}$: the difference between the two terms is $o\left({{\lambda}^{0}}\right)$ and approaches zero at large times at least as $1/{{t}^{1+\alpha}}$. In this way we have reduced Hamiltonian (7.54) to (7.1), provided that the condition (see first equation of (7.10)) $\tilde g= \texttt{g}+4\lambda m^z_\infty$ is satisfied. With this choice we recover the perturbative results of [60], e.g. the relaxation exponents (7.53). A posteriori we note that the large time behaviour of $\langle \sigma _{\ell}^{z}\rangle$ under (7.1) (see (7.43)) is sufficiently fast ($\alpha =1/2$) to justify the approximation of (7.57) by (7.59). We also checked that the mean-field solution of (7.1) is perfectly compatible with the results shown in figure 1 of [60] (see figure 13).

Zoom In Zoom Out Reset image size

Finally, we point out that [60] introduced the term $\overline{\underset{\ell}{\sum}\,\sigma _{\ell}^{z}}$ to fix some conditions in the long-time limit, where (7.1) and the Hamiltonian of [60] turn out to be equivalent. (7.1) is therefore a sensible replacement for the Hamiltonian of [60].

7.3. Generalisations

The construction of (low-entangled) stationary states that we proposed in the previous section ((7.45) and below) can be applied also to other Hamiltonians of the form (4.9) if the corresponding mean-field Hamiltonian is integrable. In those cases we can rule out thermalisation if the self-consistent problem satisfied by the stationary solutions at fixed energy and mean-field parameters has more than one solution.

When the mean-field Hamiltonian is noninteracting, following the lines of the proof sketched for the nonlocal generalisation of the Ising model in section 7.2, one can generally show that the solution is not unique.

It is also reasonable to expect that, also in the presence of interactions, the finite number of constraints given by the energy conservation and the late time values of the mean-field coupling constants could not reduce the parameter space of the initial Hamiltonian (7.45) to a single point. We indeed believe that thermalisation is unlikely to emerge if the mean-field Hamiltonian describes an integrable model at any time.

Nevertheless, the interacting case exhibits counterintuitive behaviours, for example in the energy level-spacing statistics. Generally integrable models exhibit Poisson statistics, whereas generic models follow a Wigner distribution [83, 84]. There are many exceptions to this rule [85], however the nearest neighbour spacing distribution is probably the most reliable numerical check of integrability.

In figure 14 the level spacing distribution is shown for various chain sizes for the Hamiltonian (7.1) with $\lambda =0.5$ and $\tilde{g}=0.5$ in the reflection symmetric sector of the zero momentum subspace with spin-flip parity $\underset{\ell}{\prod}\,\sigma _{\ell}^{z}$ equal to 1. The numerical data suggest that in the thermodynamic limit the curves collapse to an exponential distribution (Poisson statistics). This is consistent with our observation that at arbitrarily large times after a quantum quench the system keeps retaining infinite information about the initial state.

Zoom In Zoom Out Reset image size

A completely different scenario appears for the Hamiltonian

$$\begin{eqnarray}&&H=-\underset{\ell}{\sum}\,\left(\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}\right)+\frac{\lambda}{4L}{{\left({\sum}_{\ell}\sigma _{\ell}^{z}\sigma _{\ell +1}^{z}\right)}^{2}}.\end{eqnarray} \tag{ 7.60 }$$

The corresponding mean-field Hamiltonian for a given one-site shift invariant initial state $\mid {{\Psi}_{0}}\rangle$ is given by

$$\begin{eqnarray}&&\quad H_{\text{MF}}^{{{\Psi}_{0}}}(t)=-\underset{\ell}{\sum}\,\left(\frac{1+\gamma}{4}\sigma _{\ell}^{x}\sigma _{\ell +1}^{x}+\frac{1-\gamma}{4}\sigma _{\ell}^{y}\sigma _{\ell +1}^{y}\right)+\frac{\lambda}{2}\langle {{\Psi}_{t}}\mid \sigma _{\ell}^{z}\sigma _{\ell +1}^{z}\mid {{\Psi}_{t}}\rangle \underset{\ell}{\sum}\,\sigma _{\ell}^{z}\sigma _{\ell +1}^{z}.\end{eqnarray} \tag{ 7.61 }$$

At fixed time, this describes an XYZ model, which is known to be integrable for any choice of the coupling constants. Therefore, assuming relaxation, the stationary properties of local observables should be described by a GGE constructed with the conservation laws of an XYZ model. We considered $\gamma =0.25$ and $\lambda =0.5$. We found large finite size corrections in the level-spacing statistics (in the same sector as before) and, in particular, a remarkable even-odd parity effect (see figure 15). However, it seems that increasing L the curves collapse to a Wigner distribution, as it commonly happens in non-integrable models.

Zoom In Zoom Out Reset image size

This might appear in contradiction with our conjecture that thermalisation should not be expected when the mean-field Hamiltonian is interacting and integrable. In fact, we have not taken into account that at late times the mean-field parameters are fixed. In the previous section we ruled out thermalisation by constructing an infinite family of stationary states with the same energy and the same mean-field parameters. If this is possible, then we should find a signature of the huge degeneracy in the level-spacing statistics by restricting the space to the excited states that lie in some shell with mean-field parameters almost fixed.

In the restricted space our preliminary analysis is indeed compatible with Poisson statistics also for the Hamiltonian (7.60). However, our data turn out to be compatible with Poisson statistics even if the mean-field Hamiltonian does not describe an integrable model. This is in contrast to our expectations that in generic systems there should not be more than a few parameters that characterise the stationary state, namely the energy and, at worst, the mean-field coupling constants at infinite time after the quench.

Our interpretation of these contradictory results is that we did not investigate sufficiently large chains, so our analysis of the energy-level statistics in the restricted space is not sufficiently indicative. We are confident that a more accurate analysis will show a different behaviour in the non-integrable case.

Finally, we point out that the situation is trickier when there are isolated points in the parameter space of the mean-field coupling constants that correspond to integrable models. For example, it is not clear to us whether or not we should expect thermalisation when at asymptotically large times the coupling constants of the mean-field Hamiltonian match the integrability points.

7.4. Summary

Pre-relaxation is a dynamical phenomenon that arises when small perturbations break symmetries that affect the late time behaviour of local observables. When the perturbation breaks (abelian) integrability, this is usually called pre-thermalisation, which is generally thought as a two-step process in which local observables experience virtual relaxation before approaching thermal-like expectation values. However the relaxation process can also be more complicated, following many steps of quasi-stationary behaviour. This happens in particular when the model is close to a non-abelian integrable point. In order to extract the pre-relaxation behaviour one must therefore identify the correct time scale of the phenomenon.

We have considered the problem of pre-relaxation after quantum quenches in weakly interacting models, starting from initial states with cluster decomposition properties. We focussed on the particular situation in which the unperturbed (one-site shift invariant) Hamiltonian has a non-abelian set of local conservation laws that break one-site shift invariance. In particular we considered interacting perturbations to the XY spin-$\frac{1}{2}$ chain and investigated both integrable extensions, like the Heisenberg XYZ model, and the effects of perturbations that break integrability.

We identified the inverse perturbation strength as a relevant time scale of pre-relaxation and studied the dynamics of local observables at times proportional to it.

Despite the model being interacting, the noninteracting structure, remnant of the unperturbed Hamiltonian and manifested in the Wick theorem, survives the pre-relaxation limit. However interactions do affect the dynamics by introducing a nontrivial time dependence in the effective noninteracting Hamiltonian that generates the time evolution. The most striking effect is probably that, even if local degrees of freedom approach stationary values, these can not be generally predicted without following the entire dynamics.

We have shown how to recast the non-equilibrium problem into a system of nonlinear differential equations involving expectation values of quasi-local operators. The system of equations has qualitatively distinct solutions, which vary from trivial stationarity to persistent oscillatory behaviour over the entire time window considered. We have not found any relevant difference between integrable and non-integrable perturbations, suggesting that the scenario of thermalisation in generic models arises at much larger times.

For the very nature of the local conservation laws of the XY model, in order to have a nontrivial time evolution the initial state must break one-site shift invariance. For a particular initial state of that kind we considered a limit in which the equations can be linearised and exhibited the analytic solution, in which one-site shift invariance is eventually restored. The regime worked out analytically shows quite clearly the importance of cluster decomposition in the non-equilibrium problem. While, as mentioned above, the pre-relaxation limit is trivial for one-site shift invariant states, a shift symmetrisation of the two-site shift invariant initial state has a nontrivial time evolution. This is because cluster decomposition has been lost with the symmetrisation. It is important to take into account such aspect when analytic predictions of late time stationary behaviour are compared with numerical data at times in which one-site shift invariance is not yet restored.

The crossover between oscillatory behaviour and relaxation is quite interesting per se. This has been the main motivation for the analysis of a simplified model that shares most of the formal aspects with the effective description of pre-relaxation in the perturbed XY model, but that, in fact, has not been derived from a pre-relaxation limit. We considered a transverse-field Ising chain with an additional nonlocal interaction proportional to the magnetisation squared per unit length. This model was already studied in [60] in the framework of a perturbation theory. We used some general properties, proven for Hamiltonians of that form, to obtain nonperturbative results and showed that in the thermodynamic limit subsystems retain infinite information about the initial state, whatever large the time is. This is not in disagreement with [60], where thermalisation was conjectured for time averages in finite systems: for this model the diagonal ensemble could not be locally equivalent to the stationary state that emerges in the thermodynamic limit when the quench parameters are compatible with relaxation.

We showed that the late time behaviour in the thermodynamic limit (which corresponds to the 'pre-thermal' behaviour of [60]) can not always be described by a stationary state. In the parameter space there are indeed 'critical lines' that separate relaxation from persistent oscillatory behaviour. We defined a limit of small quench and, in that limit, exhibited the analytic expressions for such critical lines. The appearance of oscillatory behaviour has been interpreted as a consequence of the emergence of localised excitations.

We also discussed the generalisations to other Hamiltonians in which some terms have the form of interactions with macroscopic observables, like the magnetisation squared per unit length of the model above. In particular, we ruled out thermal-like behaviour in a large class of models of that kind.

Finally, we would like to stress that our description of the pre-relaxation limit is based on a few hypotheses. In particular, we neglected some 'anomalous terms', proving only the self-consistency of the conjecture. Some preliminary checks against iTEBD simulations are confirming the validity of the assumptions [89]; however, we leave a more rigorous analysis of the regimes of validity of our approximations to future research.

We thank M Collura, F Essler and A Gambassi for useful discussions and M Marcuzzi for having provided us with the data shown in figure 13. This work was supported by the EPSRC under grants EP/I032487/1 (BB) and EP/J014885/1 (MF) and by the LabEx ENS-ICFP: ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL* (MF).

We briefly summarise some useful relations valid in noninteracting models. Additional details can be found in [65] and [86].

Let us consider a local one-site shift invariant spin-chain Hamiltonian H that is mapped to noninteracting fermions by the Jordan–Wigner transformation ($\left\{a_{\ell}^{\alpha},a_{n}^{\beta}\right\}=2{{\delta}_{\alpha \beta}}{{\delta}_{\ell n}}$)

$$\begin{eqnarray}&&a_{\ell}^{\alpha}=\left(\underset{j<\ell}{\prod}\,\sigma _{j}^{z}\right)\sigma _{\ell}^{\alpha}\quad \alpha \in \left\{x,y\right\}.\end{eqnarray} \tag{ A.1 }$$

Up to boundary terms, H reads as

$$\begin{eqnarray}H\sim \frac{1}{4}\underset{\ell ,m}{\overset{L}{\sum}}\,\left(\begin{array}{c} a_{n\ell -n+1}^{x} & a_{n\ell -n+1}^{y} & \ldots & a_{n\ell}^{x} & a_{n\ell}^{y} \\ \end{array}\right){{\left[{{\mathcal{H}}^{(n)}}\right]}_{\ell m}}\left(\begin{array}{c} a_{nm-n+1}^{x} \\ a_{nm-n+1}^{y} \\ \vdots \\ a_{nm}^{x} \\ a_{nm}^{y} \\ \end{array}\right),\end{eqnarray} \tag{ A.2 }$$

where L is the chain length and $\mathcal{H}$ is a block-circulant⁹ matrix

$$\begin{eqnarray}&&{{\left[{{\mathcal{H}}^{(n)}}\right]}_{\ell m}}=\frac{n}{L}\underset{k}{\sum}\,{{\text{e}}^{-\text{i}(m-\ell )k}}{{\mathcal{H}}^{(n)}}(k)\quad {{\text{e}}^{\text{i}Lk/n}}=1\end{eqnarray} \tag{ A.3 }$$

where

$$\begin{eqnarray}&&{{\mathcal{H}}^{{{(n)}^{\dagger}}}}(k)={{\mathcal{H}}^{{{(n)}^{{}}}\,}}(k),\quad {{\mathcal{H}}^{{{(n)}^{t}}}}(k)=-{{\mathcal{H}}^{{{(n)}^{{}}}}}(-k).\end{eqnarray}$$

The index n is a divisor of L (but in the thermodynamic limit any positive integer is allowed). We call the 2n-by-2n matrix ${{\mathcal{H}}^{(n)}}(k)$ the n-site representation of the symbol, which completely characterises the block-circulant matrix. In the following we will refer to ${{\mathcal{H}}^{(n)}}(k)$ also as the symbol of H.

For given n, the following properties hold:

• (i)  
  Any function of block-circulant matrices is a block circulant matrix, with symbol equal to the function of the respective symbols

  $$\begin{eqnarray}&&f\left(\mathcal{A},\mathcal{B},\ldots \right)\to f\left(\mathcal{A}(k),\mathcal{B}(k),\ldots \right).\end{eqnarray} \tag{ A.4 }$$
• (ii)  
  Let A and B as in (A.2). Their commutator [A,B ] has the form (A.2), with symbol equal to the commutators of the symbols

  $$\begin{eqnarray}&&\left[A,B\right]\to \left[\mathcal{A}(k),\mathcal{B}(k)\right]\end{eqnarray} \tag{ A.5 }$$
• (iii)  
  The time evolution in the Heisenberg picture under (A.2.) of a noninteracting operator A of the same form is noninteracting, with symbol

  $$\begin{eqnarray}&&{{\text{e}}^{\text{i}Ht}}A{{\text{e}}^{-\text{i}Ht}}\to {{\text{e}}^{\text{i}\mathcal{H}(k)t}}A(k){{\text{e}}^{-\text{i}\mathcal{H}(k)t}}\end{eqnarray} \tag{ A.6 }$$
• (iv)  
  (Wick theorem) The expectation value of any operator in a thermal state of (A.2) (and in any Slater determinant state) can be expressed in terms of the correlation matrix

  $$\begin{eqnarray}&&\Gamma_{m\ell}^{(n)}={{\delta}_{\ell m}}\text{I}-\langle \left(\begin{array}{c} a_{nm-n+1}^{x} \\ a_{nm-n+1}^{y} \\ \vdots \\ a_{nm}^{x} \\ a_{nm}^{y} \\ \end{array}\right)\left(\begin{array}{c} a_{n\ell -n+1}^{x} & a_{n\ell -n+1}^{y} & \ldots & a_{n\ell}^{x} & a_{n\ell}^{y} \\ \end{array}\right)\rangle .\end{eqnarray} \tag{ A.7 }$$
• (v)  
  For a thermal state with inverse temperature $\beta$ the correlation matrix is given by

  $$\begin{eqnarray}&&{{\Gamma}_{\beta}}=-\tan h\left(\frac{\beta}{2}\mathcal{H}\right).\end{eqnarray} \tag{ A.8 }$$

  Therefore, the 'thermal' ground state (which, in the presence of degeneracies, is equivalent to the incoherent superposition of the states) has correlation matrix

  $$\begin{eqnarray}&&{{\Gamma}_{\infty}}=-\text{sgn}\left(\mathcal{H}\right).\end{eqnarray} \tag{ A.9 }$$
• (vi)  
  If $\mathcal{H}$ is block-circulant, $\Gamma$ is block-circulant as well, so it is completely characterised by its symbol $\Gamma(k)$. In particular, the symbol of the thermal ground state reads

  $$\begin{eqnarray}&&{{\Gamma}_{\infty}}(k)=-\text{sgn}\left(\mathcal{H}(k)\right).\end{eqnarray} \tag{ A.10 }$$
• (vii)  
  The correlation matrix of a Slater-determinant state that time evolves under the noninteracting Hamiltonian (A.2) has the symbol

  $$\begin{eqnarray}&&\Gamma(k;t)={{\text{e}}^{-\text{i}\mathcal{H}(k)t}}\Gamma(k;0){{\text{e}}^{\text{i}\mathcal{H}(k)t}}.\end{eqnarray} \tag{ A.11 }$$
• (viii)  
  In the thermodynamic limit ($L\to \infty$), the expectation value of an operator A of the form (A.2) in a Slater-determinant state with correlation matrix $\Gamma$ can be written as follows

  $$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\,\frac{1}{L}\langle A\rangle =\frac{1}{4n}{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}\text{tr}\left[\Gamma(k)\mathcal{A}(k)\right].\end{eqnarray} \tag{ A.12 }$$

In this appendix we show the validity of property 4.1. The property can be more easily proven for Jordan–Wigner fermions

$$\begin{eqnarray}&&{{c}_{j}}=\frac{1}{2}\left(a_{j}^{x}-\text{i}a_{j}^{y}\right);\end{eqnarray} \tag{ B.1 }$$

the relation for the Majorana fermions $a_{j}^{x},a_{j}^{y}$ will then follow by linearity.

In order to proceed it is convenient to introduce the following notation

$$\begin{equation}{\bf{c}}_j^\alpha (t) \equiv \left\{ {\begin{array}{*{20}{l}} {c_j^\dagger (t)}&{\alpha = + }\\ {{c_j}(t)}&{\alpha = - .} \end{array}} \right. \end{equation} \tag{ B.2 }$$

The relation between ${{c}^{\dagger}},c$ and the Bogoliubov fermions ${{b}^{\dagger}}(k),b(k)$ that diagonalise the unperturbed (noninteracting) Hamiltonian H_XY (4.1) can be written as

$$\begin{eqnarray}&&\mathbf{c}_{j}^{\alpha}(t)=\frac{1}{\sqrt{L}}\underset{k}{\sum}\,\underset{\beta =\pm}{\sum}\,{{\text{e}}^{\text{i}\alpha jk}}U(k)_{\beta}^{\alpha}{{\mathbf{b}}_{\beta}}(k){{\text{e}}^{\text{i}\beta {{\varepsilon}_{k}}t}}.\end{eqnarray} \tag{ B.3 }$$

Here U(k) is the $2\times 2$ matrix defining the Bogoliubov transformation, ${{\varepsilon}_{k}}$ is the dispersion relation in the 1-site shift invariant representation

$$\begin{eqnarray}&&{{\varepsilon}_{k}}=J\sqrt{{{\cos}^{2}}k+{{\gamma}^{2}}{{\sin}^{2}}k}\end{eqnarray} \tag{ B.4 }$$

and we set

$$\begin{equation}{\bf{c}}_j^\alpha (t) \equiv \left\{ {\begin{array}{*{20}{l}} {c_j^\dagger (t)}&{\alpha = + }\\ {{c_j}(t)}&{\alpha = - .} \end{array}} \right. \end{equation} \tag{ B.5 }$$

The relation (4.6) is then equivalent to

$$\begin{eqnarray}&&\frac{1}{L}\overline{\underset{j}{\sum}\,\mathbf{c}_{j+n1}^{{{\alpha}_{1}}}\mathbf{c}_{j+n2}^{{{\alpha}_{2}}}\mathbf{c}_{j+n3}^{{{\alpha}_{3}}}\mathbf{c}_{j+n4}^{{{\alpha}_{4}}}}=\mathcal{F}_{\left\{n\right\}}^{\left\{\alpha \right\}}+\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}},\end{eqnarray} \tag{ B.6 }$$

where $\mathcal{F}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ is a factorised term

$$\begin{equation} {\mathcal F_{\{n\}}^{\{\alpha\}}}=\sum_{s=0}^1\underbrace{{\bf c}_{n_1}^{\alpha_1} {\bf c}_{n_2}^{\alpha_2}}_{s}\underbrace{{\bf c}_{n_3}^{\alpha_3} {\bf c}_{n_4}^{\alpha_4}}_{s}-\underbrace{{\bf c}_{n_1}^{\alpha_1} {\bf c}_{n_3}^{\alpha_3}}_{s}\underbrace{{\bf c}_{n_2}^{\alpha_2} {\bf c}_{n_4}^{\alpha_4}}_{s}+\underbrace{{\bf c}_{n_2}^{\alpha_2} {\bf c}_{n_3}^{\alpha_3}}_{s}\underbrace{{\bf c}_{n_1}^{\alpha_1} {\bf c}_{n_4}^{\alpha_4}}_{s} \end{equation} \tag{ B.7 }$$

$$\begin{equation} \underbrace{{\bf c}_{n_1}^{\alpha} {\bf c}_{n_2}^{\beta}}_s= \overline{ \frac{1}{L}\sum_\ell (-1)^{s\ell}{\bf c}_{\ell+n_1}^{\alpha} {\bf c}_{\ell+n_2}^{\beta }} \end{equation} \tag{ B.8 }$$

and $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ is the remaining contribution. Using (B.3) we can explicitly carry out the time average and the sum over $\ell$ in (B.8). We obtain

$$\begin{equation} \underbrace{{\bf c}_{n_1}^{\alpha} {\bf c}_{n_2}^{\beta}}_s= \frac{1}{L}\sum_{k} {\rm e}^{-{\rm i} \alpha (n_2-n_1) k}{\rm e}^{{\rm i} n_2 s\pi}U(k)^{\alpha}_{\gamma}U({\alpha\beta}\bar k_s)^{\beta}_{\bar\gamma}{\bf b}_\gamma(k){\bf b}_{\bar\gamma}({\alpha\beta}\bar k_s), \end{equation} \tag{ B.9 }$$

where $\bar{\alpha}=-\alpha$ and we defined

$$\begin{eqnarray}&&{{k}_{s}}=k+\pi s.\end{eqnarray} \tag{ B.10 }$$

Analogously, (B.6) reads as

$$\begin{eqnarray} &&\frac{1}{L}\overline{\sum_j {\bf c}_{j+n_1}^{\alpha_1} {\bf c}_{j+n_2}^{\alpha_2} {\bf c}_{j+n_3}^{\alpha_3} {\bf c}_{j+n_4}^{\alpha_4}}\, =\,\frac{1}{L^2}\sum_{\{k_i\}} \prod_{j=1}^4\left( {\rm e}^{{\rm i} \alpha_j n_j k_j} U(k_j)^{\alpha_j}_{\beta_j}\right){\bf b}_{\beta_1}(k_1){\bf b}_{\beta_2}(k_2){\bf b}_{\beta_3}(k_3){\bf b}_{\beta_4}(k_4)\\ &&\qquad\qquad\times\delta_{\alpha_1 k_1+\alpha_2 k_2+\alpha_3 k_3+\alpha_4 k_4}\delta_{\beta_1\varepsilon_{k_1}+\beta_2\varepsilon_{k_2}+\beta_3\varepsilon_{k_3}+\beta_4\varepsilon_{k_4}}. \end{eqnarray} \tag{ B.11 }$$

In order to compute the sums over the momenta it is necessary to solve the constraints given by the delta functions, i.e.

$$\begin{eqnarray}&&{{\alpha}_{1}}{{k}_{1}}+{{\alpha}_{2}}{{k}_{2}}+{{\alpha}_{3}}{{k}_{3}}+{{\alpha}_{4}}{{k}_{4}}=0 \\&&{{\beta}_{1}}{{\varepsilon}_{{{k}_{1}}}}+{{\beta}_{2}}{{\varepsilon}_{{{k}_{2}}}}+{{\beta}_{3}}{{\varepsilon}_{{{k}_{3}}}}+{{\beta}_{4}}{{\varepsilon}_{{{k}_{4}}}}=0.\end{eqnarray} \tag{ B.12 }$$

Some of the solutions to these equations can be found by requiring the terms of (B.12.) to cancel in pairs. This would give

$$\begin{eqnarray}&&{{\delta}_{{{\alpha}_{1}}{{k}_{1}}+{{\alpha}_{2}}{{k}_{2}}+{{\alpha}_{3}}{{k}_{3}}+{{\alpha}_{4}}{{k}_{4}}}}\,{{\delta}_{{{\beta}_{1}}{{\varepsilon}_{{{k}_{1}}}}+{{\beta}_{2}}{{\varepsilon}_{{{k}_{2}}}}+{{\beta}_{3}}{{\varepsilon}_{{{k}_{3}}}}+{{\beta}_{4}}{{\varepsilon}_{{{k}_{4}}}}}}=\underset{s=0,1}{\sum}\,\Delta _{1}^{s}+\Delta _{2}^{s}+\Delta _{3}^{s},\end{eqnarray} \tag{ B.13 }$$

with

$$\begin{eqnarray}&&\Delta _{1}^{s}={{\delta}_{{{\beta}_{1}},{{{\bar{\beta}}}_{2}}}}\,{{\delta}_{{{\beta}_{3}},{{{\bar{\beta}}}_{4}}\,}}{{\delta}_{{{{\bar{\alpha}}}_{1}}{{k}_{1}},{{\alpha}_{1}}{{k}_{2,s}}}}\,{{\delta}_{{{{\bar{\alpha}}}_{3}}{{k}_{3}},{{\alpha}_{4}}{{k}_{4,s}}}}\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&\Delta _{2}^{s}={{\delta}_{{{\beta}_{1}},{{{\bar{\beta}}}_{3}}}}\,{{\delta}_{{{\beta}_{2}},{{{\bar{\beta}}}_{4}}}}\,{{\delta}_{{{{\bar{\alpha}}}_{1}}{{k}_{1}},{{\alpha}_{3}}{{k}_{3,s}}}}\,{{\delta}_{{{{\bar{\alpha}}}_{2}}{{k}_{2}},{{\alpha}_{4}}{{k}_{4,s}}}}\end{eqnarray} \tag{ B.15 }$$

$$\begin{eqnarray}&&\Delta _{3}^{s}={{\delta}_{{{\beta}_{1}},{{{\bar{\beta}}}_{4}}}}\,{{\delta}_{{{\beta}_{2}},{{{\bar{\beta}}}_{3}}}}\,{{\delta}_{{{{\bar{\alpha}}}_{1}}{{k}_{1}},{{\alpha}_{4}}{{k}_{4,s}}}}\,{{\delta}_{{{{\bar{\alpha}}}_{2}}{{k}_{2}},{{\alpha}_{3}}{{k}_{3,s}}}}.\end{eqnarray} \tag{ B.16 }$$

For a generic dispersion relation it is reasonable to expect these solutions to be the only possible. For the specific dispersion considered, (B.4), equations (B.12) admit other solutions in the thermodynamic limit. We call these anomalous solutions because they strongly depend on the precise form of the dispersion relation. We now show that the term $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ in (B.6) is exactly the contribution arising from these solutions, i.e.

$$\begin{eqnarray}\mathcal{F}_{\left\{n\right\}}^{\left\{\alpha \right\}}&=&\frac{1}{{{L}^{2}}}\underset{{{k}_{1}},{{k}_{2}},{{k}_{3}},{{k}_{4}}}{\sum}\,\underset{s=0,1}{\sum}\,\underset{j=1}{\overset{4}{\prod}}\,\left({{\text{e}}^{\text{i}{{\alpha}_{j}}{{n}_{j}}{{k}_{j}}}}U\left({{k}_{j}}\right)_{{{\beta}_{j}}}^{{{\alpha}_{j}}}\right){{\mathbf{b}}_{{{\beta}_{1}}}}\left({{k}_{1}}\right){{\mathbf{b}}_{{{\beta}_{2}}}}\left({{k}_{2}}\right){{\mathbf{b}}_{{{\beta}_{3}}}}\left({{k}_{3}}\right){{\mathbf{b}}_{{{\beta}_{4}}}}\left({{k}_{4}}\right)\\&&\times \left\{\Delta _{1}^{s}+\Delta _{2}^{s}+\Delta _{3}^{s}\right\}.\end{eqnarray} \tag{ B.17 }$$

We stress that the operator $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ will be nonzero only in the thermodynamic limit.

We consider for example the term containing ${{\Delta }_{2}}$. We have

$$\begin{eqnarray} &&\frac{1}{L^2}\sum_{k_1, k_2, k_3, k_4} \sum_{s=0,1}\prod_{j=1}^4\left( {\rm e}^{{\rm i} \alpha_j n_j k_j} U(k_j)^{\alpha_j}_{\beta_j}\right) {\bf b}_{\beta_1}(k_1){\bf b}_{\beta_2}(k_2){\bf b}_{\beta_3}(k_3){\bf b}_{\beta_4}(k_4)\Delta^s_2\,\\ &&\qquad=\frac{1}{L^2}\sum_{p,q}\sum_{s=0,1}{\rm e}^{{\rm i} \alpha_1 (n_1-n_3) p + {\rm i} \alpha_2 (n_2-n_4) q}{\rm e}^{{\rm i} n_3 s \pi}{\rm e}^{{\rm i} n_4 s \pi}U(p)^{\alpha_1}_{\beta_1}U(q)^{\alpha_2}_{\beta_2}U({\alpha_1\alpha_3}\bar p_{s})^{\alpha_3}_{\bar\beta_1}U({\alpha_2\alpha_4}\bar q_{s})^{\alpha_4}_{\bar\beta_2}\\ &&\qquad\times{\bf b}_{\beta_1}(p){\bf b}_{\beta_2}(q){\bf b}_{\bar\beta_1}({\alpha_1\alpha_3}\bar p_{s}){\bf b}_{\bar\beta_2}({\alpha_2\alpha_4} \bar q_{s})\\ &&\qquad=-\sum_{s=0,1}\underbrace{{\bf c}_{n_1}^{\alpha_1} {\bf c}_{n_3}^{\alpha_3}}_s\underbrace{{\bf c}_{n_2}^{\alpha_2} {\bf c}_{n_4}^{\alpha_4}}_s+\,O\left(L^{-1}\right), \end{eqnarray} \tag{ B.18 }$$

where we used the commutation relations of the $\left\{{{\mathbf{b}}_{\beta}}(k)\right\}$ in the last step. Although the terms $\mathcal{F}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ and $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ are in fact multiplied by L in the time average of (4.2), the possible corrections O(L⁻¹) in (B.18.) (which would result in corrections O(L⁰) in the effective Hamiltonian) are locally irrelevant, because their density approaches zero in the thermodynamic limit.

We obtain analogous results for ${{\Delta }_{1}}$ and ${{\Delta }_{3}}$, that is to say (B.6).

Remark 2. We point out that for other dispersion relations (still with the properties ${{\varepsilon}_{k}}={{\varepsilon}_{k+\pi}}$ and ${{\varepsilon}_{k}}\ne {{\varepsilon}_{k+\pi /n}}$ for generic k and n  >  1) the anomalous terms could be factorised as well. Generally in such situations the factors have a very simple time dependence, e.g. a single oscillation frequency. As a consequence, relaxation is ruled out.

In this appendix we prove the lemmas of section 5.

Lemma C.1. If $\mathcal{O}\in \mathcal{E}$, the operator norm (i.e. the maximal eigenvalue in absolute value) of $\mathcal{O}/L$ is bounded.

Proof. The proof is straightforward. Let us expand $\mathcal{O}/L$ as in (4.9):

$$\begin{eqnarray}&&\frac{\mathcal{O}}{L}=\underset{j=1}{\overset{N}{\sum}}\,\frac{1}{{{L}^{{{n}_{j}}}}}\mathcal{O}_{1}^{(j)}\cdots \mathcal{O}_{{{n}_{j}}}^{(j)}\end{eqnarray} \tag{ C.1 }$$

where $\mathcal{O}_{m}^{(j)}$ have local densities, that is to say, they can be written as follows

$$\begin{eqnarray}&&\mathcal{O}_{m}^{(j)}=\underset{\ell}{\sum}\,o_{m;\ell}^{(j)},\end{eqnarray} \tag{ C.2 }$$

with $o_{m;\ell}^{(j)}$ local operators. We immediately find the chain of inequalities

$$\begin{eqnarray}&&\parallel \frac{\mathcal{O}}{L}\parallel \leqslant \underset{j=1}{\overset{N}{\sum}}\,\frac{1}{{{L}^{{{n}_{j}}}}}\parallel \mathcal{O}_{1}^{(j)}\cdots \mathcal{O}_{{{n}_{j}}}^{(j)}\parallel \leqslant \underset{j=1}{\overset{N}{\sum}}\,\parallel \frac{\mathcal{O}_{1}^{(j)}}{L}\parallel \cdots \parallel \frac{\mathcal{O}_{{{n}_{j}}}^{(j)}}{L}\parallel \leqslant \underset{j=1}{\overset{N}{\sum}}\,{{\left(\underset{m,\ell}{{\max}}\,\parallel o_{m;\ell}^{(j)}\parallel \right)}^{{{n}_{j}}}}.\end{eqnarray} \tag{ C.3 }$$

The right hand side is clearly O(L⁰) because N and n_j are finite by definition and $o_{m;\ell}^{(j)}$ are local.

Lemma C.2. If $\mathcal{O},\tilde{\mathcal{O}}\in \mathcal{E}$, then $\left[\mathcal{O},\tilde{\mathcal{O}}\right]\in \mathcal{E}$ as well.

Proof. Without loss of generality, we can restrict to two single terms of the expansions (4.9) of $\mathcal{O}$ and $\tilde{\mathcal{O}}$. The commutator of the two terms reads as

$$\begin{array}{l} [\frac{1}{L^{n_i-1}}\mathcal O_1^{(i)}\cdots \mathcal O_{n_i}^{(i)},\frac{1}{L^{n_j-1}}\tilde{\mathcal O}_1^{(j)}\cdots \tilde{\mathcal O}_{n_j}^{(j)}]\nonumber \\ \qquad=\frac{1}{L^{n_j+n_j-2}}\sum_{k,p}\mathcal O_1^{(i)}\cdots \mathcal O_{k-1}^{(i)}\tilde{\mathcal O}_1^{(j)}\cdots \tilde{\mathcal O}_{p-1}^{(j)}[\mathcal O_k^{(i)},\tilde{\mathcal O}_p^{(j)}]\mathcal O_{k+1}^{(i)} \cdots\mathcal O_{n_i}^{(i)}\tilde{\mathcal O}_{p+1}^{(j)} \cdots\tilde{\mathcal O}_{n_j}^{(j)}. \end{array} \tag{ C.4 }$$

Since $\left[\mathcal{O}_{k}^{(i)},\tilde{\mathcal{O}}_{p}^{(j)}\right]$ have local densities (the commutator of two local operators is nonzero only if there is a region on which they both act nontrivially; in addition, its range is smaller than the sum of the ranges of the two operators), the number of extensive operators exceeds by one the exponent of 1/L. Thus, $\left[\mathcal{O},\tilde{\mathcal{O}}\right]\in \mathcal{E}$.

Lemma C.3. (viz. lemma 5.1) Let $\mathcal{O}\in \mathcal{E}$ and $\mid \Psi\rangle$ a state with cluster decomposition properties. The expectation value of $\mathcal{O}/L$ in $\mid \Psi\rangle$ can be reduced to the expectation values of the local translation invariant operators it consists of:

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\,\langle \Psi\mid \frac{{{H}_{1}}}{L}\cdots \frac{{{H}_{n}}}{L}\mid \Psi\rangle =\underset{L\to \infty}{{\lim}}\,\underset{j}{\prod}\,\frac{\langle \Psi\mid {{H}_{j}}\mid \Psi\rangle}{L}.\end{eqnarray} \tag{ C.5 }$$

Proof. Let us consider a term (5.1) of the expansion (4.9). Its expectation value (per unit of length) is given by

$$\begin{eqnarray}&&\langle \Psi\mid \frac{{{H}_{1}}}{L}\cdots \frac{{{H}_{n}}}{L}\mid \Psi\rangle =\frac{1}{{{L}^{n}}}\underset{{{\ell}_{1}},\ldots ,{{\ell}_{n}}}{\sum}\,\langle \Psi\mid {{h}_{1,{{\ell}_{1}}}}\cdots {{h}_{n,{{\ell}_{n}}}}\mid \Psi\rangle ,\end{eqnarray} \tag{ C.6 }$$

where ${{h}_{j,{{\ell}_{j}}}}$ are local operators acting nontrivially only around ${{\ell}_{j}}$ and such that

$$\begin{eqnarray}&&{{H}_{j}}=\underset{\ell}{\sum}\,{{h}_{j,{{\ell}_{j}}}}.\end{eqnarray} \tag{ C.7 }$$

By cluster decomposition we have

$$\begin{eqnarray}&&\underset{\begin{array}{c} {{\ell}_{1}},\ldots ,{{\ell}_{n}} \\ \mid {{\ell}_{j}}-{{\ell}_{{{j}^{\prime}}}}\mid >\xi \gg 1\left(\forall j\ne {{j}^{\prime}}\right) \\ \end{array}}{\sum}\,\langle \Psi\mid \frac{{{h}_{1,{{\ell}_{1}}}}}{L}\cdots \frac{{{h}_{n,{{\ell}_{n}}}}}{L}\mid \Psi\rangle =\underset{\begin{array}{c} {{\ell}_{1}},\ldots ,{{\ell}_{n}} \\ \mid {{\ell}_{j}}-{{\ell}_{{{j}^{\prime}}}}\mid >\xi \gg 1\left(\forall j\ne {{j}^{\prime}}\right) \\ \end{array}}{\sum}\,\frac{\langle \Psi\mid {{h}_{1,{{\ell}_{1}}}}\mid \Psi\rangle}{L}\cdots \frac{\langle \Psi\mid {{h}_{n,{{\ell}_{n}}}}\mid \Psi\rangle}{L}+f(\xi ,L),\end{eqnarray} \tag{ C.8 }$$

where

$$\begin{eqnarray}&&\underset{\xi \to \infty}{{\lim}}\,\underset{L\to \infty}{{\lim}}\,f(\xi ,L)=0.\end{eqnarray} \tag{ C.9 }$$

The difference between (C.6) and the left hand side of (C.8) can be bounded from above as follows

$$\begin{eqnarray}&&\mid \underset{\begin{array}{c} {{\ell}_{1}},\ldots ,{{\ell}_{n}} \\ \mid {{\ell}_{j}}-{{\ell}_{{{j}^{\prime}}}}\mid \leqslant \xi \text{}\text{}\left(\exists j\ne {{j}^{\prime}}\right) \\ \end{array}}{\sum}\,\langle |\Psi\frac{{{h}_{1,{{\ell}_{1}}}}}{L}\cdots \frac{{{h}_{n,{{\ell}_{n}}}}}{L}|\Psi\rangle \mid \leqslant \left(\begin{array}{c} n \\ 2 \\ \end{array}\right)\frac{\xi}{L}\underset{\left\{\ell \right\}}{{\max}}\,\mid \langle \Psi\mid {{h}_{1,{{\ell}_{1}}}}\cdots {{h}_{n,{{\ell}_{n}}}}\mid \Psi\rangle \mid \to 0.\end{eqnarray} \tag{ C.10 }$$

Analogously

$$\begin{eqnarray}&&\mid \underset{\begin{array}{c} {{\ell}_{1}},\ldots ,{{\ell}_{n}} \\ \mid {{\ell}_{j}}-{{\ell}_{{{j}^{\prime}}}}\mid \leqslant \xi \text{}\text{}\left(\exists j\ne {{j}^{\prime}}\right) \\ \end{array}}{\sum}\,\frac{\langle \Psi\mid {{h}_{1,{{\ell}_{1}}}}\mid \Psi\rangle}{L}\cdots \frac{\langle \Psi\mid {{h}_{n,{{\ell}_{n}}}}\mid \Psi\rangle}{L}\mid \leqslant \left(\begin{array}{c} n \\ 2 \\ \end{array}\right)\frac{\xi}{L}\underset{\left\{\ell \right\}}{{\max}}\,\mid \underset{j=1}{\overset{n}{\prod}}\,\langle \Psi\mid {{h}_{j,{{\ell}_{j}}}}\mid \Psi\rangle \mid \to 0\text{}\text{},\end{eqnarray} \tag{ C.11 }$$

so that

$$\begin{eqnarray}&&\mid \langle \Psi\mid \frac{{{H}_{1}}}{L}\cdots \frac{{{H}_{n}}}{L}\mid \Psi\rangle -\underset{j}{\prod}\,\frac{\langle \Psi\mid {{H}_{j}}\mid \Psi\rangle}{L}\mid \leqslant \mid f(\xi ,L)\mid +O(1/L)\end{eqnarray} \tag{ C.12 }$$

Being $\xi$ arbitrary, we can take the limit ${{\lim}_{\xi \to \infty}}\,\,\underset{L\to \infty}{{\lim}}\,$, obtaining (C.5).

Lemma C.4. If the state $\mid {{\Psi}_{0}}\rangle$ has cluster decomposition properties and $\mathcal{O}\in \mathcal{E}$, the mean-field Hamiltonian defined in section 5 satisfies the following identity:

$$\begin{eqnarray}&&\quad \underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\left[H{{_{\text{MF}}^{{{\Psi}_{0}}}}^{{}}}_{{}}(t),\frac{\mathcal{O}}{L}\right]\bar{U}(t)\mid {{\Psi}_{0}}\rangle =\underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\left[H,\frac{\mathcal{O}}{L}\right]\bar{U}(t)\mid {{\Psi}_{0}}\rangle ,\end{eqnarray} \tag{ C.13 }$$

where $\bar{U}$ was defined in (5.4).

Proof. Let us consider a generic term (5.1) of the expansion (4.9) of H

$$\begin{eqnarray}&&\tilde{H}=\frac{1}{{{L}^{n-1}}}{{H}_{1}}\cdots {{H}_{n}}.\end{eqnarray} \tag{ C.14 }$$

The corresponding term (5.3) of the mean-field Hamiltonian (5.3) is given by

$$\begin{eqnarray}&&\tilde{H}_{\text{MF}}^{{{\Psi}_{0}}}(t)=\underset{\ell =1}{\overset{n}{\sum}}\,\underset{j\ne \ell}{\prod}\,\frac{\langle {{\Psi}_{0}}\mid {{{\bar{U}}}^{\dagger}}(t){{H}_{j}}\bar{U}(t)\mid {{\Psi}_{0}}\rangle}{L}{{H}_{\ell}}.\end{eqnarray} \tag{ C.15 }$$

By taking the commutators with $\mathcal{O}$ we find

$$\begin{eqnarray} &&\left[\tilde H,\frac{\mathcal O}{L}\right]=\sum_{\ell=1}^n \prod_{j=1}^{\ell-1} \frac{H_j}{L}\frac{[H_{\ell},\mathcal O]}{L}\prod_{j=\ell+1}^n \frac{H_j}{L} \end{eqnarray} \tag{ C.16 }$$

$$\begin{eqnarray} &&\left[\tilde H_{\rm MF}^{\Psi_0}(t),\frac{\mathcal O}{L}\right]=\sum_{\ell=1}^n \prod_{j\neq \ell} \frac{\langle{\Psi_0\vert \bar U^\dagger(t)H_j\bar U(t)\vert \Psi_0}\rangle}{L}\frac{[H_\ell,\mathcal O]}{L}. \end{eqnarray} \tag{ C.17 }$$

Because $\mid {{\Psi}_{0}}\rangle$ has cluster decomposition properties and the mean-field Hamiltonian is local at any time, the state $\bar{U}(t)\mid {{\Psi}_{0}}\rangle$ has cluster decomposition properties as well (the only difference with respect to $\mid {{\Psi}_{0}}\rangle$ is that the function f of (C.8) is now time dependent). Finally, by lemma C.3, in the thermodynamic limit the expectation values of (C.16) and (C.17) in the state $\bar{U}(t)\mid {{\Psi}_{0}}\rangle$ are identical, that is to say (C.13).

Lemma C.5. Let $\mid {{\Psi}_{0}}\rangle$ be a translation invariant state with cluster decomposition properties and $H,\mathcal{O}\in \mathcal{E}$. The time derivatives of the expectation value of $\mathcal{O}/L$ in the state evolving with $H_{\text{MF}}^{{{\Psi}_{0}}}(t)$ fulfil

$$\begin{equation} \quad\frac{{\rm d}^{\,n}}{{\rm dt}^{\,n}}\lim_{L\rightarrow\infty}\langle{\Psi_0\vert \bar U^\dagger(t)\frac{\mathcal O}{L}\bar U(t) \vert \Psi_0}\rangle=i^n\lim_{L\rightarrow\infty}\langle{\Psi_0\vert \bar U^\dagger(t)\underbrace{[H,[H,...[H}_{n},\frac{\mathcal O}{L}]...]\bar U(t) \vert \Psi_0}\rangle \end{equation} \tag{ C.18 }$$

Proof. We proceed by induction. First of all we see that for n   =   0 the property is trivially satisfied; let then the property be true for n, we have

$$\begin{array}{l} \frac{{\rm d}^{\,n+1}}{{{\rm d}t}^{\,n+1}}\lim_{L\rightarrow\infty}\langle{\Psi_0\vert \bar U^\dagger(t)\frac{\mathcal O}{L}\bar U(t) \vert \Psi_0}\rangle=i^n\frac{{\rm d}}{{\rm d}t}\lim_{L\rightarrow\infty}\langle{\Psi_0\vert \bar U^\dagger(t)\underbrace{[H,[H,...[H,}_{n}\frac{\mathcal O}{L}]...]\bar U(t) \vert \Psi_0}\rangle\nonumber\\ \qquad\qquad\qquad\,\,\, = i^{n+1}\lim_{L\rightarrow\infty}\langle{\Psi_0\vert \bar U^\dagger(t)[H_{\rm MF}(t),\underbrace{[H,[H,...[H}_{n},\frac{\mathcal O}{L}]...]]\bar U(t) \vert \Psi_0}\rangle\nonumber\\ \qquad\qquad\qquad\,\,\, = i^{n+1}\lim_{L\rightarrow\infty}\langle{\Psi_0\vert \bar U^\dagger(t)\underbrace{[H,[H,...[H}_{n+1},\frac{\mathcal O}{L}]...]\bar U(t) \vert \Psi_0}\rangle. \end{array} \tag{ C.19 }$$

In the second step we used lemmas C.2 and C.4. This concludes the proof.

Lemma C.6. (viz. lemma 5.2) Let $\mid {{\Psi}_{0}}\rangle$ be a translation invariant state with cluster decomposition properties and $H,\mathcal{O}\in \mathcal{E}$. Let the expectation value of $\mathcal{O}$ in the state that time evolves with $H_{\text{MF}}^{{{\Psi}_{0}}}(t)$ be an analytic function of t in the strip $\mid \text{Im}\left[t\right]\mid <r$, with r a nonzero constant. In the thermodynamic limit, the time evolution with H can be replaced by the time evolution with the mean-field Hamiltonian:

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}Ht}}\frac{\mathcal{O}}{L}{{\text{e}}^{-\text{i}Ht}}\mid {{\Psi}_{0}}\rangle =\underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\frac{\mathcal{O}}{L}\bar{U}(t)\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ C.20 }$$

Proof. We define

$$\begin{eqnarray}&&f(t,s)=\underset{L\to \infty}{{\lim}}\,\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t){{\text{e}}^{\text{i}Hs}}\frac{\mathcal{O}}{L}{{\text{e}}^{-\text{i}Hs}}\bar{U}(t)\mid {{\Psi}_{0}}\rangle .\end{eqnarray} \tag{ C.21 }$$

By lemma C.5 we have

$$\begin{eqnarray}&&\frac{{{\partial}^{n}}}{\partial {{t}^{n}}}f(t,0)=\frac{{{\partial}^{n}}}{\partial {{s}^{n}}}{{\mid}_{s=0}}f(t,s),\end{eqnarray} \tag{ C.22 }$$

indeed

$$\begin{equation} \left.\frac{\partial^n}{\partial t^n}\right\vert _{t=0} {\rm e}^{{\rm i} H t}\frac{\mathcal O}{L} {\rm e}^{-{\rm i} H t}=i^n \underbrace{[H,[H,...[H}_{n},\frac{\mathcal O}{L}]...]. \end{equation} \tag{ C.23 }$$

By assumption, f (t, 0) (which corresponds to the time evolution with the mean-filed Hamiltonian) is analytic in the strip $\mid \text{Im}\left[t\right]\mid <r$, so the convergence radius of the Taylor expansion at t   =   0 is larger than or equal to r. Thus we have

$$\begin{eqnarray}&&f(\tau ,0)={\sum}_{n}\frac{{{\tau}^{n}}}{n!}\frac{{{\partial}^{n}}}{\partial {{t}^{n}}}{{\mid}_{t=0}}f(t,0)={\sum}_{n}\frac{{{\tau}^{n}}}{n!}\frac{{{\partial}^{n}}}{\partial {{t}^{n}}}{{\mid}_{t=0}}f(0,t)=f(0,\tau ),\quad \mid \tau \mid <r.\end{eqnarray} \tag{ C.24 }$$

Let us call ${{t}_{*}}$ a time such that f (t, 0)   =    f (0, t) for any $0\leqslant t<{{t}_{*}}$. As before, the function $f(t+\tau ,0)$ is analytic in the strip $\mid \text{Im}\left[\tau \right]\mid <r$, so we have

$$\begin{eqnarray}&&f(t+\tau ,0)=\underset{n}{\sum}\,\frac{{{\tau}^{n}}}{n!}\frac{{{\partial}^{n}}}{\partial {{t}^{n}}}f(t,0)=\underset{n}{\sum}\,\frac{{{\tau}^{n}}}{n!}\frac{{{\partial}^{n}}}{\partial {{t}^{n}}}f(0,t)=f(0,t+\tau ),\quad \mid \tau \mid <r.\end{eqnarray} \tag{ C.25 }$$

That is to say

$$\begin{equation}\qquad\qquad f(t,0)=f(0,t) \quad \forall t< t_\ast\quad\Longrightarrow\quad f(t,0)=f(0,t) \quad \forall t< t_\ast+\tau. \end{equation} \tag{ C.26 }$$

Since $\tau$ is finite and (C.24) holds, we conclude

$$\begin{eqnarray}&&f(t,0)=f(0,t)\quad \forall t,\end{eqnarray} \tag{ C.27 }$$

which is exactly (C.20).

Corollary C.7. (viz. corollary 5.3) Lemma C.6 holds true in particular for local operators.

Proof. By translation invariance, the expectation value of any local operator $\mathcal{O}$ is equal to the expectation value per unit of length of the operator ${{\mathcal{O}}_{*}}\in \mathcal{E}$, obtained by shifting $\mathcal{O}$ along the chain and summing all the (L) terms.

Corollary C.8. (viz. corollary 5.4) Let $\mid {{\Psi}_{0}}\rangle$ a translation invariant state with cluster decomposition properties and $H\in \mathcal{E}$. In the thermodynamic limit, the time evolution of the reduced density matrix (RDM) of some spin block S is equal to the RDM in the state that time evolves with the mean-field Hamiltonian:

$$\begin{eqnarray}&&{{\rho}_{S}}(t)=\text{t}{{\text{r}}_{{\bar{S}}}}\left[{{\text{e}}^{-\text{i}Ht}}\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid {{\text{e}}^{\text{i}Ht}}\right]=\text{t}{{\text{r}}_{{\bar{S}}}}\left[\bar{U}(t)\mid {{\Psi}_{0}}\rangle \langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\right].\end{eqnarray} \tag{ C.28 }$$

Proof. This is a direct consequence of corollary C.7.

Corollary C.9. (viz. corollary 5.5) Let $H\in \mathcal{E}$ and $\mid \Psi\rangle$ a state with cluster decomposition properties. If $\mid \Psi\rangle$ is an excited state of the corresponding mean-field Hamiltonian $H_{\text{MF}}^{\Psi}$

$$\begin{eqnarray}&&H_{\text{MF}}^{\Psi}\mid \Psi\rangle ={{E}_{\Psi}}\mid \Psi\rangle ,\end{eqnarray} \tag{ C.29 }$$

the expectation value of local observables in ${{\text{e}}^{-\text{i}Ht}}\mid \Psi\rangle$ is independent of time. Therefore, $\mid \Psi\rangle$ behaves locally as an excited state of H.

The reverse is also true. If an excited state of H is locally equivalent to a state with cluster decomposition properties, then the latter is an excited state of the corresponding mean-field Hamiltonian.

Proof. Clearly the mean-field Hamiltonian $H_{\text{MF}}^{\Psi}$ is the solution of (5.3). Being ${{\text{e}}^{-\text{i}H_{\text{MF}}^{\Psi}t}}\mid \Psi\rangle \propto \mid \Psi\rangle$, by corollary C.7 the expectation value of local observables is independent of time. The reverse holds true for analogous reasons.

Here we show that neglecting the anomalous term $L\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ (see (B.6)) in the time averaged Hamiltonian is a self-consistent approximation. To this aim, we consider the time evolution of a Slater determinant $\mid {{\Psi}_{0}}\rangle$ under the Hamiltonian $\tilde{H}$, obtained from H by removing $L\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$. In section 5 and appendix C we proved that, as long as $\mathcal{O}$ is a local operator (but the class of allowed operators is in fact larger), ${{\text{e}}^{-\text{i}\tilde{H}t}}$ can be replaced by the mean-field time evolution operator $\bar{U}(t)$ (5.4). Here we show that inserting $L\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ back at time t does not change the expectation value of local observables. In other words we are going to prove

$$\begin{eqnarray}&&\underset{L\to \infty}{{\lim}}\,L\langle {{\Psi}_{0}}\mid {{\bar{U}}^{\dagger}}(t)\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}},\mathcal{O}\right]\bar{U}(t)\mid {{\Psi}_{0}}\rangle =0\quad \forall t,\end{eqnarray} \tag{ D.1 }$$

where $\mathcal{O}$ is a generic local operator.

Using the notations of appendix B, any local operator can be written as a linear combination of operators of the form

$$\begin{eqnarray}&&\quad \mathcal{O}=\mathbf{c}_{{{\ell}_{1}}}^{{{\gamma}_{1}}}\cdots \mathbf{c}_{{{\ell}_{n}}}^{{{\gamma}_{n}}}=\frac{1}{{{L}^{n/2}}}\underset{\left\{{{p}_{i}}\right\}}{\sum}\,{{\mathcal{F}}_{n}}_{\left\{{{\sigma}_{i}}\right\}}^{\left\{{{\gamma}_{i}}\right\}}\left(\left\{{{\ell}_{i}}\right\}\mid \left\{{{p}_{i}}\right\}\right){{\mathbf{b}}_{{{\sigma}_{1}}}}\left({{p}_{1}}\right)\cdots {{\mathbf{b}}_{{{\sigma}_{n}}}}\left({{p}_{n}}\right),\end{eqnarray} \tag{ D.2 }$$

where

$$\begin{eqnarray}&&{{\mathcal{F}}_{n}}_{\left\{{{\beta}_{i}}\right\}}^{\left\{{{\alpha}_{i}}\right\}}\left(\left\{{{j}_{i}}\right\}\mid \left\{{{p}_{i}}\right\}\right)\equiv \underset{i=1}{\overset{n}{\prod}}\,\left({{\text{e}}^{\text{i}{{\alpha}_{i}}{{j}_{i}}{{p}_{i}}}}U\left({{p}_{i}}\right)_{{{\beta}_{i}}}^{{{\alpha}_{i}}}\right).\end{eqnarray} \tag{ D.3 }$$

If n is odd then (D.1) is trivially satisfied because $\bar{U}(t)\mid {{\Psi}_{0}}\rangle$ is a Slater determinant by assumption and hence the expectation value of an odd number of fermions vanishes. We therefore focus on the case n   =   2m. The anomalous term $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ of (B.6) can be written as follows

$$\begin{eqnarray}&&\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}=\frac{1}{{{L}^{2}}}\underset{{{k}_{1}},{{k}_{2}}}{\sum}\,\underset{{{{\bar{k}}}_{3}},{{{\bar{k}}}_{4}}}{\sum}\,{{\mathcal{F}}_{2}}_{\left\{{{\beta}_{i}}\right\}}^{\left\{{{\alpha}_{i}}\right\}}\left(\left\{{{j}_{i}}\right\}\mid \left\{{{k}_{i}}\right\}\right){{\mathcal{F}}_{2}}_{\left\{{{\beta}_{i}}\right\}}^{\left\{{{\alpha}_{i}}\right\}}\left(\left\{{{j}_{i}}\right\}\mid \left\{{{\bar{k}}_{i}}\right\}\right){{\mathbf{b}}_{{{\beta}_{1}}}}\left({{k}_{1}}\right){{\mathbf{b}}_{{{\beta}_{2}}}}\left({{k}_{2}}\right){{\mathbf{b}}_{{{\beta}_{3}}}}\left({{\bar{k}}_{3}}\right){{\mathbf{b}}_{{{\beta}_{4}}}}\left({{\bar{k}}_{4}}\right)\end{eqnarray} \tag{ D.4 }$$

where ${{\bar{k}}_{3}}$ and ${{\bar{k}}_{4}}$ are the anomalous solutions of system (B.12), i.e. they are implicit functions of k₁ and k₂ defined by the system (B.12) and in addition fulfilling

$$\begin{eqnarray}&&{{k}_{1,2}}\pm {{\bar{k}}_{3,4}}\ne 0,\quad {{\bar{k}}_{3}}\pm {{\bar{k}}_{4}}\ne 0,\end{eqnarray} \tag{ D.5 }$$

almost everywhere. Since $\vert{\tilde{\Psi}_t}\rangle=\bar U(t)\vert{\Psi_0}\rangle$ is a Slater determinant, we can use the Wick theorem to compute expectation values. We then have

$$\begin{eqnarray}&&L\langle {{{\tilde \Psi}}_{t}}\mid \mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\mid {{{\tilde \Psi}}_{t}}\rangle =L\langle {{{\tilde \Psi}}_{t}}\mid \mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mid {{{\tilde \Psi}}_{t}}\rangle \langle {{{\tilde \Psi}}_{t}}\mid \mathcal{O}\mid {{{\tilde \Psi}}_{t}}\rangle +{{\mathcal{C}}_{2}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}+{{\mathcal{C}}_{4}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}\end{eqnarray} \tag{ D.6 }$$

where ${{\mathcal{C}}_{2}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}$ contains terms in which two of the $\mathbf{b}$'s in $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ are contracted together and the other two are contracted with two $\mathbf{b}$'s in $\mathcal{O}$; ${{\mathcal{C}}_{4}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}$ contains all the terms in which any $\mathbf{b}$ in $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$ is contracted with a $\mathbf{b}$ in $\mathcal{O}$.

According to the definition of $\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}$, any Wick contraction among $\mathbf{b}$'s in it gives zero (because of (D.5)), hence the only non zero contribution to (D.6) arises from ${{\mathcal{C}}_{4}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}$. To conclude the proof we will show that the terms in ${{\mathcal{C}}_{4}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}$ scale as O(L⁻¹) and in the thermodynamic limit their contribution can thus be neglected. To this end it is sufficient to consider a typical element of ${{\mathcal{C}}_{4}}{{\left[\mathcal{A}_{\left\{n\right\}}^{\left\{\alpha \right\}}\mathcal{O}\right]}_{t}}$

where sums over the momenta {p_i}, {k_i} and the indices are understood, and we defined

$$\begin{eqnarray}&&G_{\mathbf{j},\mathbf{k},\boldsymbol{\ell },\mathbf{p}}^{\boldsymbol{\alpha },\boldsymbol{\beta },\boldsymbol{\gamma },\boldsymbol{\sigma }}\equiv {{\mathcal{F}}_{2}}_{\left\{{{\beta}_{i}}\right\}}^{\left\{{{\alpha}_{i}}\right\}}\left(\left\{{{j}_{i}}\right\}\mid \left\{{{k}_{i}}\right\}\right){{\mathcal{F}}_{2}}_{\left\{{{\beta}_{i}}\right\}}^{\left\{{{\alpha}_{i}}\right\}}\left(\left\{{{j}_{i}}\right\}\mid \left\{{{\bar{k}}_{i}}\right\}\right){{\mathcal{F}}_{n}}_{\left\{{{\sigma}_{i}}\right\}}^{\left\{{{\gamma}_{i}}\right\}}\left(\left\{{{\ell}_{i}}\right\}\mid \left\{{{p}_{i}}\right\}\right).\end{eqnarray} \tag{ D.9 }$$

The 2m+2 sums over the momenta are reduced to m by the Kronecker deltas arising from the Wick contractions. Because the number of factors L⁻¹ exceeds by one the number of sums, the term turns out to be O(L⁻¹). The validity of equation (D.1) is then established.

E.1. Ground state and maximal energy state

Here we show that the TFIC ground state $\mid {{\psi}_{{{g}_{0}}}}\rangle$, with g₀ satisfying (7.16), is the Slater determinant that minimises the energy (7.11) of the Hamiltonian (7.1).

The symbol of the correlation matrix of the most general reflection symmetric one-site shift invariant Slater determinant is

$$\begin{eqnarray}&&\Gamma(k)={{n}_{x}}(k){{\sigma}^{x}}+{{n}_{y}}(k){{\sigma}^{y}}+{{n}_{z}}(k){{\sigma}^{z}},\end{eqnarray} \tag{ E.1 }$$

with ${{n}_{\alpha}}(k)$ real functions that satisfy

$$\begin{eqnarray}&&n_{x}^{2}+n_{y}^{2}+n_{z}^{2}=1\end{eqnarray} \tag{ E.2 }$$

$$\begin{eqnarray}&&{{n}_{x,z}}(k)=-{{n}_{x,z}}(-k)\end{eqnarray} \tag{ E.3 }$$

$$\begin{eqnarray}&&{{n}_{y}}(k)={{n}_{y}}(-k).\end{eqnarray} \tag{ E.4 }$$

Equation (E.2) manifests the fact that the initial state is pure, which indeed implies that the eigenvalues of $\Gamma(k)$ are ±1. The other conditions simply mean that the correlation matrix is a purely imaginary skew-symmetric matrix. The absence of a term proportional to the identity in (E.1) is a consequence of reflection symmetry.

The energy (7.11) can be written as follows:

$$\begin{eqnarray}&&\quad \varepsilon =\lambda {{\left({\int}_{0}^{\pi}\frac{\text{d}k}{\pi}{{n}_{y}}(k)\right)}^{2}}-\tilde{g}{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}{{n}_{y}}(k)+{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\left(\cos k{{n}_{y}}(k)-\sin k{{n}_{x}}(k)\right).\end{eqnarray} \tag{ E.5 }$$

The minimisation can be worked out by zeroing the variation of the functional

$$\begin{eqnarray}&&\Phi\left[{{n}_{x}},{{n}_{y}},{{n}_{z}},\mu \right]=\varepsilon -{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\mu (k)\left(n_{x}^{2}(k)+n_{y}^{2}(k)+n_{z}^{2}(k)\right)\end{eqnarray} \tag{ E.6 }$$

with respect to its arguments. We immediately see that the variation with respect to n_z results in

$$\begin{eqnarray}&&\mu (k){{n}_{z}}(k)=0.\end{eqnarray} \tag{ E.7 }$$

If we assume $\mu (k)\ne 0$, then n_z(k)   =   0 and we can enforce (E.2) by setting ${{n}_{x}}=-\sin\theta$ and ${{n}_{y}}=\cos\theta$. Instead of working with the functional $\Phi$ we can therefore express the energy as

$$\begin{eqnarray}&&\varepsilon \left[\theta \right]=\lambda {{\left({\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\cos{{\theta}_{k}}\right)}^{2}}-\tilde{g}{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\cos{{\theta}_{k}}+{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\left(\cos k\cos{{\theta}_{k}}+\sin k \sin{{\theta}_{k}}\right)\end{eqnarray} \tag{ E.8 }$$

and consider its variation with respect to ${{\theta}_{k}}$. We obtain

$$\begin{eqnarray}&&\left[\tilde{g}-2\lambda \left({\int}_{0}^{\pi}\frac{\text{d}p}{\pi}\cos{{\theta}_{p}}\right)-\cos k\right]\sin{{\theta}_{k}}+\sin k \cos{{\theta}_{k}}=0\end{eqnarray} \tag{ E.10 }$$

which is solved by

$$\begin{eqnarray}&&{{\text{e}}^{\text{i}{{\theta}_{k}}}}=\frac{{{g}_{0}}-{{\text{e}}^{\text{i}k}}}{\sqrt{1+g_{0}^{2}-2{{g}_{0}}\cos k}},\end{eqnarray} \tag{ E.10 }$$

where g₀ satisfies (7.16). One of the solutions of (7.16) minimises the energy. The energy is instead maximal for

$$\begin{eqnarray}&&{{\text{e}}^{\text{i}{{\theta}_{k}}}}=-\frac{{{g}_{1}}-{{\text{e}}^{\text{i}k}}}{\sqrt{1+g_{1}^{2}-2{{g}_{1}}\cos k}},\end{eqnarray} \tag{ E.11 }$$

where g₁ is a solution of

$$\begin{eqnarray}&&\tilde{g}={{g}_{1}}-2\lambda \left({\int}_{0}^{\pi}\frac{\text{d}p}{\pi}\frac{{{g}_{1}}-\cos p}{\sqrt{1+g_{1}^{2}-2{{g}_{1}}\cos p}}\right).\end{eqnarray} \tag{ E.12 }$$

We emphasise that numerical data obtained by exact diagonalization of (7.1) for small chains are compatible with the two states of (E.10) and (E.11) being the actual ground state and the maximal energy state, respectively, of (7.1) in the thermodynamic limit.

E.2. Excited states

Following corollary 5.5, we construct an infinite number of states for which the expectation value of local observables that time evolve under the Hamiltonian (7.1) is independent of time. To all intents and purposes, those states can be considered eigenstates of (7.1).

A generic stationary solution of (7.7) should have ${{\tilde{y}}^{\prime}}=0$ and

$$\begin{eqnarray}&&4(h-\cos k){{\phi}_{k}}+16(h\cos k-1){{y}_{k}}=0,\end{eqnarray} \tag{ E.13 }$$

namely

$$\begin{eqnarray}&&-16(h-\cos k)\text{tr}\left[\Gamma(k){{\sigma}^{y}}{{\text{e}}^{\text{i}k{{\sigma}^{z}}}}\right]+16(h\cos k-1)\text{tr}\left[\Gamma(k){{\sigma}^{y}}\right]=0.\end{eqnarray} \tag{ E.14 }$$

This is solved by

$$\begin{eqnarray}&&{{\Gamma}_{m}}(k)=-m(k)\frac{\left({{h}_{m}}-\cos k\right){{\sigma}^{y}}+\sin k{{\sigma}^{x}}}{\sqrt{1+h_{m}^{2}-2{{h}_{m}}\cos k}},\end{eqnarray} \tag{ E.15 }$$

where m(k)   =   m(−k) and we used ${{\tilde{y}}^{\prime}}=0$ to remove any term proportional to ${{\sigma}^{z}}$. The parameter h_m is the solution of the equation

$$\begin{eqnarray}&&{{h}_{m}}=\tilde{g}+2\lambda {\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}m(k)\frac{{{h}_{m}}-\cos k}{\sqrt{1+h_{m}^{2}-2{{h}_{m}}\cos k}}.\end{eqnarray} \tag{ E.16 }$$

We recognise ${{\Gamma}_{m}}(k)$ as the symbol of the correlation matrix of a reflection symmetric (Slater determinant) excited state of the TFIC model with magnetic field h_m [67, 82, 87, 88]. Within this interpretation, if h_m was independent of m, $m\left({{k}_{j}}\right)=\pm 1$ (for finite chains) would have provided an orthonormal basis of reflection symmetric excited states. We remind the reader that in the thermodynamic limit infinite excited states become locally equivalent and, in sufficiently regular cases, can be characterised by a function $-1\leqslant m(k)\leqslant 1$. However, in the latter characterisation, which exactly describes the local properties of the excited states, the manifest orthonormality properties of the basis are lost. This problem is even more pronounced in our case, where h_m depends on the excitation, and a naive finite volume regularisation does not produce an orthonormal basis also for $m\left({{k}_{j}}\right)=\pm1$ (this is not unexpected, being our description valid only in the thermodynamic limit). Nevertheless, the noninteracting states with correlation matrices (E.15) are locally equivalent to the excited states of (7.1). In addition, it seems that they span the entire space of reflection symmetric states.

We notice that in the thermodynamic limit generally there are infinite m that give rise to the same h_m. Indeed we can always add to m(k) an even function $\delta m(k)$ such that $-1\leqslant m+\delta m\leqslant 1$ and

$$\begin{eqnarray}&&{\int}_{0}^{\pi}\frac{\text{d}k}{\pi}\delta m(k)\frac{{{h}_{m}}-\cos k}{\sqrt{1+h_{m}^{2}-2{{h}_{m}}\cos k}}=0,\end{eqnarray} \tag{ E.17 }$$

which physically means the the excitation must have the same magnetisation per unit of length in the z direction.

For the sake of completeness we report the expression of the energies of the excited states:

$$\begin{eqnarray}&&{{\varepsilon}_{m}}=\frac{h_{m}^{2}-{{{\tilde{g}}}^{2}}}{4\lambda}-{\int}_{-\pi}^{\pi}\frac{\text{d}k}{2\pi}m(k)\frac{{{h}_{m}}\cos k-1}{\sqrt{1+h_{m}^{2}-2{{h}_{m}}\cos k}}.\end{eqnarray} \tag{ E.18 }$$

We conclude this appendix with a comparison between finite chain data and expectations based on the mean-field correspondence, which strongly relies on the thermodynamic limit. In particular, we tried to apply corollary C.9 to a finite chain to estimate the lowest and highest excitations. Figures E1 and E2 show that the mean-field predictions are in fairly good agreement with numerical data even in small chains (notice that there can be $\mathcal{O}(1/L)$ corrections in both directions).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size