Quantum quenches in the anisotropic spin-1/2 Heisenberg chain: different approaches to many-body dynamics far from equilibrium

Recent experimental achievements in controlling ultracold gases in optical lattices open a new perspective on quantum many-body physics. In these experimental setups it is possible to study coherent time evolution of isolated quantum systems. These dynamics reveal new physics beyond the low-energy properties usually relevant in solid-state many-body systems. In this paper we study the time evolution of antiferromagnetic order in the Heisenberg chain after a sudden change of the anisotropy parameter, using various numerical and analytical methods. As a generic result we find that the order parameter, which can show oscillatory or non-oscillatory dynamics, decays exponentially except for the effectively non-interacting case of the XX limit. For weakly ordered initial states we also find evidence for an algebraic correction to the exponential law. The study is based on numerical simulations using a numerical matrix product method for infinite system sizes (iMPS), for which we provide a detailed description and an error analysis. Additionally, we investigate in detail the exactly solvable XX limit. These results are compared to approximative analytical approaches including an effective description by the XZ-model as well as by mean-field, Luttinger-liquid and sine-Gordon theories. This reveals which aspects of non-equilibrium dynamics can as in equilibrium be described by low-energy theories and which are the novel phenomena specific to quantum quench dynamics. The relevance of the energetically high part of the spectrum is illustrated by means of a full numerical diagonalization of the Hamiltonian.


4
in this paper we intend to support ongoing efforts at improving the control of ultracold atomic gases.
We will study the emerging dynamics of the order parameter of an XXZ Heisenberg chain prepared in the classical (uncorrelated) Néel state, which can be realized in experiment, but in order to get a deeper insight into the problem, general antiferromagnetic initial states are also considered. Our special interest concerns the effect of the quantum phase transition that can be triggered by tuning the magnetic anisotropy parameter.
Exact results based on numerical calculations are presented. Furthermore, alternative approximative approaches are applied. The applicability of the analytical tools, which have been very successful in the description of equilibrium phenomena, turns out to be strongly restricted for the non-equilibrium problem under consideration. We identify the apparent problems in the standard approximations and point out in which direction these approaches should be extended in order to capture the main features of the quantum quench dynamics.

Brief review of quantum quenches in extended systems
In relation to transport phenomena (e.g. [36]) and impurity problems such as quantum dots [37] and spins in a dissipative environment [38], non-equilibrium dynamics have been subjected to intensive theoretical investigations over many years. However, non-equilibrium transport can be seen as a result of perturbations (voltage biases) at the edges of the system and quantum dots are zero-dimensional systems. This is fundamentally different from quench dynamics in translation-invariant systems considered here, where the parameter change is global and the energy scales involved in the dynamics scale with the system size. More closely related to a quantum quench are highly excited electronic states in solids, generated in femtosecond pump-probe spectroscopy [39]- [41]. Nevertheless, in these systems, decoherence times are short and the dissipative processes strongly contribute to the emerging dynamics. Consequently, concepts developed for transport phenomena and dynamics in condensed matter systems are not necessarily appropriate for quenches in ultracold atomic systems. Except for pioneering works on quench dynamics in the 1970s [42]- [45], specific theoretical research on quench dynamics started only recently, stimulated by the experimental developments in ultracold atomic physics. In these works, which shall be briefly summarized in this section, two main lines have been followed. A first aspect is the study of the nature of the quasi-stationary states in the long-time limit. As demonstrated by an experiment of Kinoshita et al [16], these non-equilibrium states can exhibit striking properties for specific types of interactions. Another approach explicitly focuses on the characteristics of the time evolution after the quench-experimental examples are the oscillations [14] or the dephasing [17] of the superfluid phase. It turns out to be an ambitious challenge to establish relations between dynamical phenomena and the details of the microscopic model, such as integrability and dimensionality. Although numerous remarkable theoretical efforts revealed a number of interesting phenomena, many aspects of relaxation dynamics and equilibration, which shall be discussed in detail in this work, remain unclear.
The effective description of many-body systems by means of low-energy theories, captured within the renormalization group framework [46], has proven sufficient for the theoretical understanding of a broad range of equilibrium phenomena. Therefore the application of renormalization group ideas to non-equilibrium dynamics seems to be a promising approach. Along this way, diagrammatic techniques [47]- [50] and the solutions of the dynamics of field-theoretical models at the renormalization group fixed point [51]- [59] were developed. 5 The flow-equation method as a unitary perturbative approach has been applied to Fermi gases [60,61] and the sine-Gordon model [59]. Providing a generic view on the quench problem for critical theories, the work of Calabrese and Cardy [54,55] based on conformal field theory has to be highlighted. Although, for continuum systems, field-theoretical models were successfully applied to generic quantum quenches [17,56], it has to be clarified under what conditions they all provide an accurate description of lattice systems. The range of applicability of semiclassical theories is also unclear [62]- [64]. For a restricted class of problems the time evolution can be calculated exactly, e.g. for the Jordan-Wigner diagonalizable XY chains [65,66], where the relaxation of order parameters [43,67] and correlation functions [44], [68]- [70] has been investigated for various quench scenarios with ferromagnetic interactions. Also the quench in the 1 r -Hubbard-chain [71] allows for an exact solution. A major drawback of these exactly solvable models is that the possibility of their representation in terms of non-interacting particles apparently leads to very specific relaxation phenomena, which are not generic, not only for non-integrable but also for more complicated integrable models. For instance, it is questionable whether the generalized Gibbs ensemble, which has been proposed for the description of quasi-stationary states of integrable models [69], is a useful concept beyond the simple Jordan-Wigner diagonalizable cases [70,72]. For the more general Bethe-ansatz solvable models, it has not yet been possible to extract dynamics, except for the Richardson [73] and the Lieb-Liniger models [74].
In view of the high complexity of the quench dynamics, efficient unbiased numerical approaches are crucial to gain deeper insight. Using exact diagonalization [75]- [77] it is possible to calculate the dynamics of small systems over exceedingly long times. For larger (but one-dimensional (1D)) systems the density matrix renormalization group (DMRG) [78]- [81] can be applied. Although only for finite times, the dynamics of spin-chain (respectively spinless fermions) [82]- [88] and bosonic [76], [89]- [92] and fermionic [93] lattice models have been evaluated. Recently, the dynamical mean field theory has been applied to fermionic models in the limit of infinite dimensions [94]- [96].

Basic setup and general discussion
The Heisenberg model is a paradigm in the theory of magnetism and strongly correlated systems in general. In appendix A, we derive how the model can be realized with ultracold two-level atoms in various geometries of optical lattices. For instance, it is possible to generate a one-dimensional XXZ Heisenberg model, where the sign and the strength of the exchange coupling J and can be tuned dynamically. The XXZ model is integrable and its eigenstates can be constructed by the Bethe ansatz. In the case of antiferromagnetic couplings J > 0, the anisotropy parameter triggers a quantum phase transition from a gapless 'Luttinger liquid' phase (0 < 1) to a gapped, Ising-ordered antiferromagnetic phase ( > 1). The main features of the model at equilibrium and its fieldtheoretical formulation are given in appendix B.
The non-equilibrium dynamics (1) shall be investigated in the following quantum quench: at time t < 0, the system is prepared in a ground state |ψ 0 with long-range antiferromagnetic order. The corresponding anisotropy parameter is denoted by 0 , 0 > 1. Among the 6 antiferromagnetic equilibrium states the Néel state, which corresponds to the limit 0 → ∞, has already been realized in experiment [19] and will attract our special attention. At t = 0, the system is pushed out of equilibrium by changing the strength of the interaction, < 0 , and the dynamics emerging at t > 0 are studied.
In the context of optical lattices, where the system is well isolated and no phonons are present, dissipation can be neglected in a first approximation. Also, being interested in quantum effects, we set T = 0. Finite temperature may become relevant for weak magnetic exchange interactions in the ultracold atomic setup, but how to investigate efficiently the non-equilibrium problem at T > 0 is still an unsolved problem. Under these assumptions, the dynamics is formally described by the solution of the Schrödinger equation, We seth = 1 throughout this paper. Involving a priori all the energy scales of the manybody Hamiltonian, the calculation of the time evolution of the wave function (3) is highly complex. When approaching the problem analytically, one is forced to introduce an appropriate approximation-the advantages and drawbacks of various approaches will be investigated in this work. When using numerics, the dynamics (3) can be solved by fully diagonalizing the Hamiltonian H . In section 8, we apply the full numerical diagonalization approach. Highly efficient routines have been developed for this purpose [97], which can nevertheless be used only for small system sizes (up to 20 lattice sites). More efficient and applicable directly in the thermodynamic limit are matrix product states (MPS), which will be used for the simulation of the general quench dynamics of the XXZ model. For a detailed description of the MPS method, see appendix C.
To describe the dynamics of the state |ψ(t) , we mainly focus on the antiferromagnetic order parameter, Since the state |ψ 0 is invariant under translation and subsequent spin inversion, ±m s (t) corresponds to the local magnetization at any site of the lattice. It will also be useful to look at the frequency distribution f m s ( ), which resolves the contributions to the dynamics in energy space, The staggered magnetization is not only the natural observable characterizing the ordering of antiferromagnetic states, but also reflects the properties of the local density matrix of a single site. For describing non-local properties, we choose the equal-time connected spin-spin correlation function, Before going into the study of the many-body dynamics of the Hamiltonian (1), it is worthwhile considering the case of only two spins. A corresponding experiment has been carried out by Trotzky et al [19] by loading 87 Rb atoms in the hyperfine states |↓ = |F = 1, 7 m F = −1 , |↑ = |F = 1, m F = 1 , into an array of double-well potentials. The initial Néel state was generated using a magnetic field gradient transferring the effective spins in each double well from a triplet-bond state into an antiferromagnetic one, |ψ 0 = |↑↓ . The dynamics are in this special case independent of and can be described as Rabi oscillations between |↑↓ and |↓↑ states, Hence, the antiferromagnetic order parameter describes oscillatory behavior, m s (t) = 1 2 cos(J t), where the Rabi frequency is set by the exchange coupling J , which was indeed observed in the experiment [19].
Although, as we shall see in section 2, in a many-body system such Rabi-like oscillations may survive, the dynamics become much more intricate when going to large system sizes. On a heuristic level, the initial state may be regarded as a bunch of excitations of the Hamiltonian H , whose dynamics gives rise to the propagation of correlations throughout the system. For spin models with sufficiently local interactions, Lieb and Robinson [98] have proven that this propagation takes place within a light cone-the deviation of a correlation function from its initial value becomes exponentially small for distances > 2ut, where u is the maximum velocity of excitations in the system. For an isolated but arbitrarily large system, this means that relaxation to a stationary state can only be observed for subsystems of size < 2ut. This lightcone effect has been more precisely described in the framework of boundary conformal field theory [99], which predicts an exponential decay of the correlations in the long-time limit. These short-range correlations are in contrast with the entanglement properties of the non-equilibrium problem. It has been shown [99] that the entanglement entropy of a subsystem of size grows linearly with time if 2ut < and saturates to a value proportional to if 2ut > .
It is an open question under what conditions the stationary state in the long-time limit can be described by a statistical ensemble at a finite temperature, meaning that thermalization occurs. There are several examples for which this is not the case [69,77,84,85,89], and the extended Gibbs ensemble [69], which takes into account the constraints of the non-dissipative dynamics, and the micro-canonical ensemble [72,75,92,93,100] are possible candidates for describing steady states. Whether the integrability is a necessary condition for the absence of thermalization remains unclear. It has been suggested that the absence of thermalization could be associated with non-perturbative behavior, which is not related to the integrability of the underlying Hamiltonian [77].
Here, it will be shown that in the long-time limit the antiferromagnetic order vanishes in all cases; hence, at least for this local quantity, thermalization is observed-in a one-dimensional system, no long-range order is possible at finite temperatures. This does not necessarily imply thermalization for correlation functions. Indeed, in section 3, we present results that indicate the absence of thermalization in the spin-spin correlations (6). However, the correlation functions exhibit somewhat slow relaxation dynamics and it is difficult to extract steady-state properties from the rather short accessible times that can be achieved numerically.
Nevertheless, interesting dynamical effects are present also at short times. Their characterization as a function of the initial state and the interaction parameter will be investigated. The magnetic order parameter turns out to be a good observable for the quantitative extraction of nontrivial time scales. Here, where the initial state can be characterized by the gap parameter s (more precisely, the inverse correlation length), one expects that the typical time scale of the relaxation dynamics is given by −1 s and the length scales, which depend on the 8 momentum distributions in the initial states, should be of the order of u/ s , where u is given by the velocity of quasi-particles (spin waves).
In the solution of the quench dynamics for conformally invariant theories [101] of Calabrese and Cardy [54,55], these qualitative arguments were put on a solid ground: the initial state enters into the framework of quantum field theory as a finite slab width, τ e , the extrapolation length that stands for the renormalization-group distance of the initial state from the fixed point of the gapped theory [102]. To first order, this is given by the inverse gap, here τ e ∼ −1 s . Using a conformal transformation, the slab geometry is mapped onto the semi-infinite plane, for which, by means of boundary conformal field theory [103], the properties of the correlation functions can be extracted.
The results of Calabrese and Cardy [55] do apply to the quench in the XXZ model if the discussion is restricted to the low-energy modes in the gapless regime | | 1, here captured by the Luttinger model (see appendix B, equation (B.9)). For the staggered magnetization as a local observable, the outcome is where τ e ∼ −1 s . However, several remarks concerning the applicability of the conformal field theory results to the quench in the XXZ chain are in place. Firstly, the initial state is treated on a perturbative level in terms of a renormalization-group distance from the fixed point and simply characterized by the gap parameter. It is questionable whether in this framework it is possible to correctly take into account the physics of the antiferromagnetic states, especially those close to the critical point (i.e. far from the antiferromagnetic fixed point), described by the sine-Gordon model. Secondly, within the field theory it is impossible to treat lattice effects, which are expected to emerge if the energy of the quasi-particles forming the initial state is of the order of the bandwidth -a situation that is realized for instance by the Néel state (2). As a simple example of a lattice effect we presented the Rabi oscillations in the two-spin system (7), with the frequency set by the magnetic exchange J . Macroscopic order parameter oscillations following a quantum quench have been predicted to appear in a variety of systems [62,67], [104]- [106]. In this work, we will characterize Rabi-like oscillations and investigate origins of dephasing in the presence of many-body correlations. A particular property of the quench in the XXZ chain illustrates the novel aspect of the non-equilibrium dynamics in many-body lattice models: the time evolution of m s (t) is invariant under the change of sign → − . Ferro-and antiferromagnetic Hamiltonians exhibit identical dynamics despite their completely different elementary excitations. As a third point restricting the applicability of the conformal field theory result, we mention that a conformal theory does not capture the case of a parameter quench into the gapped phase, > 1. Here this regime will be addressed using a sine-Gordon description of the XXZ model.

Summary of the results
The further content of the paper is organized as follows: the non-equilibrium dynamics in the XX limit of the Heisenberg model, which can be solved in a simple way by means of the Jordan-Wigner transformation, is analyzed in section 2. Numerical results for the general case are given in section 3 and approximative approaches in sections 4-6. In section 8, an exact diagonalization analysis of the spectrum of the XXZ model is carried out before presenting Section 2: Exact analytical calculation in the XX limit Néel the conclusions. In appendix A, we describe the experimental realization of quantum magnetic systems in optical lattices. The well-established properties of antiferromagnetic states and equilibrium phase transitions in the context of the Heisenberg model in one dimension are reviewed in appendix B. The description and an error analysis of the matrix product algorithm is provided in appendix C.
Our results for the non-equilibrium dynamics of the staggered magnetization are summarized in table 1. We find essentially two types of relaxation dynamics: non-oscillatory dynamics, characterized by a relaxation time τ 1 , and oscillatory dynamics with a frequency ω and an associated relaxation time τ 2 . An important result is that for nonzero , we find a fundamentally new mode of many-body dynamics, which always leads to exponential decay of the staggered moment regardless of whether the short-time dynamics is oscillatory or not. In contrast with the oscillation frequency, which is set by the exchange interaction, the relaxation time is an emergent scale generated by the highly correlated dynamics and, hence, cannot be simply related to the microscopic parameters. We find divergent relaxation times, τ 1 → ∞ in the limit → 0 and τ 2 → ∞ if → ∞. For the particular case of the Néel state, we find that the relaxation times essentially vanish in the vicinity of the critical point, 1. Table 1 also shows to what extent approximative methods, which take into consideration only a particular aspect of the Hamiltonian, are applicable to the non-equilibrium problem.
The mean-field approximation, for example, leads to contradictions with the unbiased numerical results-an algebraic decay for 1 and a non-vanishing asymptotic value of the staggered moment for > 1 [104]. In the case of the initial Néel state, comparing the low-energy result of conformal field theory with the numerics, the immediate relaxation τ 1 ≈ 0 is, in principle, in agreement with s → ∞ in equation (8). However, the oscillations dominate the long-time dynamics and are, as expounded before, not captured by the field theory. If the initial state is close to the critical point, an exponential relaxation similar to equation (8) is found; however, an additional algebraic prefactor appears to be present. In our treatment of the Luttinger model this effect is also not seen, but the results from conformal field theory (8) are reproduced.

The XX model, = 0
It is particularly illustrative to study the exactly solvable case of zero anisotropy ( = 0), where the Heisenberg Hamiltonian (1) can be represented in terms of free spinless fermions with a cosine dispersion relation (B.4). For free fermions the non-equilibrium dynamics can be solved analytically [107]. We study two cases: firstly, the Néel state as the initial condition and, secondly, the case of the initial spin-density-wave (SDW) state.

Initial Néel state, 0 = ∞
In the fermionic picture, the Néel state reads as a charge density wave, The fermionic operators are easily represented in the Heisenberg picture, Hence, the dynamics of the XX chain prepared in the Néel state, in analogy with the twosite model (7), takes the form of Rabi oscillations between SDW with different sublattice magnetizations, The relaxation of the staggered magnetization can be seen as a dephasing process, driven by inhomogeneous Rabi frequencies in k-space, In the thermodynamic limit, where J 0 denotes the zeroth Bessel function of the first kind. The underlying frequency distribution (5) ranges over a band of width 4J , θ( ) being the Heaviside function. High-energy modes with a vanishing velocity at the band edge, | k | = J , dominate the long-time limit of (13) and give rise to the oscillations with a frequency set by the bandwidth, The exponent of the t −1/2 decay is a consequence of the quadratic dispersion at k = 0. In the XX limit it is also possible to express the correlation function, G zz c ( , t), in terms of Bessel functions, This results in slowly decaying, spatially oscillating correlations, Figure 1 shows how the correlations evolve within the light cone 2t. The magnitude of the wave front decays as a power law in time. The negative sign reflects spinon characteristics [108] of the propagating correlations.
Although it is possible to carry out the analysis of the XX model without any approximation, it is useful to investigate the result of restriction to a particular part of the spectrum. This provides information on the range of applicability of low-energy theories, which are candidates for treating the more complicated case of interacting systems.

12
In the case of the linearized theory (appendix B, equation (B.8)), the dynamics of the magnetization is characterized by oscillations with a 1/t decay and cutoff-dependent period, The cutoff gives the correct periodic behavior if it is equal to the bare bandwidth ( = 2J ). The oscillatory behavior, a consequence of the presence of the lattice, is indeed not captured in the continuum limit /J → ∞, where the oscillations disappear. The power-law decay appears in the linear approximation being independent of the cutoff, but the exponent is overestimated by a factor of two compared to the case of the full dispersion. The energy distribution corresponding to the magnetization (17) is simply flat, A seemingly (in the context of equilibrium theories) unconventional approach is the development of the modes in the vicinity of the band edges, In the present case of non-equilibrium dynamics, we find, however, that the corresponding energy distribution, provides the correct long-time limit if the cutoff is sufficiently large, We now clearly understand the mechanism behind the dephasing process in the freefermion models: Rabi oscillations are present if there is a sharp step at the edge of the band. The dephasing of the oscillations is algebraic, t −α , α = 1 if the frequency distribution is homogeneous and α = 1 2 in the case of the quadratic dispersion at the band edge. For the longtime behavior it is sufficient to stick to the modes at the edge of the band; the low-frequency part is effective only at short times t ∼ J −1 . The reason for such behavior is best illustrated in the analysis of the correlation functions for the linear spectrum. The result, as shown in figure 1, is a single coherent spinon mode traveling in the light cone |2t − | = 0. For the staggered magnetization as a local observable, this means that it relaxes as soon as the spinon mode moves over more than one lattice distance 2t > 1. In contrast to the case of the full dispersion, there are no oscillations within the light cone. We note that this immediate decay is in agreement with the result of conformal field theory (8), which predicts zero relaxation time for the Néel state due to its vanishing correlation length (inverse gap).

Initial spin-density wave
As an introduction to our discussion of quenches from correlated antiferromagnetic states (i.e. quenches with 1 < 0 < ∞), we consider the time evolution of weakly antiferromagnetic SDW states under the XX Hamiltonian (see appendix B, equations (B.4) and (B.20)). This section will provide a benchmark for the numerical results in section 3 and also discusses the applicability of effective low-energy theories to this quench. 13 The time evolution of the staggered magnetization m s (t) in the XX model starting from an SDW state at t = 0 ( where we have taken the thermodynamic limit in the last equation. With the coefficients obeying , the dephasing process in energy representation reads For a weak SDW state ( s 1) there are two main contributions to the integral in equation (23). The first comes from the Fermi points = 0, whereas the second originates in the square root singularities at = ±J . We write these two contributions separately, In comparison with the case of the initial Néel state, in addition to identical algebraically decaying oscillations (13) a non-oscillatory decay stemming from the low-energy part of the spectrum is obtained. This exponential behavior with an algebraic prefactor is characterized by the relaxation time τ = (2 s ) −1 . Hence, for t > −1 s ln(J/ s ) the oscillations on top of the non-oscillatory decay dominate the order-parameter dynamics. Nevertheless, unlike the case of the initial Néel state, the low-energy modes contribute to the non-equilibrium dynamics over significant periods of time.

Interaction quench in the XXZ model-numerical study
In this section, we first study the quench in the XXZ model starting from the Néel state. Subsequently, ground states of the XXZ models at finite = 0 will be considered.
Unlike for = 0, the problem is no longer analytically treatable and we have to resort to numerical techniques. In the iMPS algorithm (appendix C) we use 2000 states and a secondorder Suzuki-Trotter decomposition with a time step δ ∼ 10 −3 J −1 for large and up to 7000 states with δ ∼ 10 −2 J −1 for small . An intermediate time regime J t 16 can be reached, which exceeds in general greatly the short transient time.

Initial Néel state, 0 = ∞
An overview of the results for the initial Néel state is presented in figure 2. For small anisotropies, we find oscillations of the order parameter similar to those in the XX limit, but with the decay time decreasing upon approaching the isotropic point = 1. In the easy-axis regime > 1 of the XXZ model, the relaxation slows down again for increasing and we observe non-oscillatory behavior for 1. Figure 3 focuses on easy-plane anisotropy 0 < < 1. The results for 0 < 0.4 are well described, for accessible time scales, by exponentially decaying oscillations   The oscillation frequency is almost independent of the anisotropy, while the relaxation time τ 2 increases with decreasing . Logarithmic divergence of the relaxation time in the limit → 0 is suggested by the fit shown in figure 4. The picture is less clear closer to the isotropic point. For the range 0.5 < 1 there appears to be an additional time scale after which the oscillations start to decay even faster than exponentially; simultaneously the period of the oscillations is reduced. Therefore, the relaxation times plotted in figure 4 are only valid within an intermediate time window, whose width shrinks upon approaching the critical point.
For intermediate easy-axis anisotropies 1 3, the magnetization does not reach a stable regime within the numerically accessible time window (figure 5(a)). The complicated behavior of m s (t) in this parameter range can be ascribed to the interplay of processes at all energy scales. Nevertheless, the numerical data suggest that the relaxation is fastest close to the isotropic point, in the range between = 1 and = 1.6. A simple generic type of behavior is recovered for large anisotropies 3. The numerical data in figure 5(b) indicate exponential relaxation of the staggered magnetization The relaxation time scales roughly quadratically with (figure 4). Oscillations do persist on top of the exponential decay, but they fade out quickly. We briefly describe the relaxation of the spin-spin correlation functions (6) as presented in figure 6. A more detailed study of these has been carried out by Manmana et al [87]. For weak interactions (e.g. = 0.6) the dynamics of correlation functions is still dominated by the spinon mode moving according to the light cone [98,99] set by the spin-wave velocity u (see appendix B, equation (B.11)), as is the case at = 0 (equation (16), figure 1). For larger , this mode is smeared out; instead, antiferromagnetic correlations build up. The strength of the short-range antiferromagnetic correlations increases as the anisotropy is augmented. With the numerical method, however, we are unable to reach sufficiently long times to calculate the  25) and (26)). (b) Comparison of the XXZ chain (symbols) and the XZ chain (dashed lines) for strong anisotropies; solid lines correspond to an exponential fit. The dynamics of the staggered magnetization of the XXZ and XZ chains converge towards each other in the large-limit.
quasi-stationary correlation length. It becomes, nevertheless, clear that the correlations cannot be described in terms of a thermal ensemble. We evaluated the equilibrium correlation functions at a temperature corresponding to the energy of the system by means of quantum Monte Carlo simulations 7 . The resulting correlation functions depicted in figure 6 decay considerably faster than the non-equilibrium ones.

Initial antiferromagnet,
The Néel state is an entirely classical state with no quantum correlations. In order to generalize our results, we first study the case of small but finite correlations starting from the ground state for 0 = 4.0. We find that the picture gained from the initial Néel state remains qualitatively valid-the dynamics of m s (t) is very similar to that in the case of the initial Néel state (figure 2). The corresponding relaxation times and periods are plotted in figure 7. For close to zero, the behavior of τ 2 is again close to a logarithmic law and the divergence of the relaxation time for T T denotes quantum Monte Carlo results for the XXZ model at equilibrium at a temperature fixed by the energy of the non-equilibrium system.

T T T T T T T T T T T T T T T T T T
We expect qualitatively different behavior for a weakly ordered (more strongly correlated) initial state. In section 2, we have seen that for an initial SDW state and the XX Hamiltonian, in addition to the algebraically decaying oscillations, an exponential relaxation exists, whose relaxation rate is proportional to the gap of the initial state. In figure 8, where we show the results for the quench from an initial state with 0 = 1.5, oscillations are found on top of non-oscillatory relaxation. At = 0, for sufficiently large t, the dynamics is similar to the SDW result (24), to very high accuracy, despite the fact that the SDW is a different wave function from the ground state of the XXZ chain. The relaxation time τ 1 ≈ 5.1J is slightly smaller than the one predicted by the SDW calculations, (2 s ) −1 ≈ 5.8J . The difference may be explained by the importance of short-range effects, which are supposed to contribute to the non-equilibrium dynamics. As illustrated in figure B.3, the correlations decay much faster at shorter distances than in the large distance asymptotics.
For 0, in correspondence with the result for the initial Néel state, we find that the oscillations are exponentially damped, while the non-oscillatory part remains qualitatively the same as in the XX limit, Figure 9. Relaxation times and oscillation period T = 2π/ω as a function of anisotropy in the XXZ model for the system prepared in the ground state at 0 = 1.5. Solid lines are a guide to the eye. τ 1 is comparable to (J/K s ).
In figure 9, we plot the fitting parameters for small (0 < 0.6) where formula (28) is well obeyed. τ 1 behaves similarly to (J/K s )-a law that is the natural extension of the non-interacting SDW result (24) to finite anisotropies using the same scaling as derived for the Luttinger model (see section 6, equation (48)). The logarithmic behavior of τ 2 , apparent for 0 1, is not observed here. Oscillatory and non-oscillatory terms are superimposed. In the non-oscillatory term of (28), absence of the algebraic prefactor, as suggested by the field theoretical result (8), can be clearly excluded on the basis of the numerical results. Pure exponential law (26) is, however, found for 1. The intermediate regime 0.6 1 cannot be described by either of the laws (26) and (28).

Mean field
Time-dependent mean-field theory is one possible way to treat the dynamics of the XXZ model approximately. The mean-field approximation of the Hamiltonian (B.2) at an instant of time t is defined by expanding the interaction term to linear order in fluctuations δn j around the mean density, n j = n j + δn j , and by setting n j = 1/2 + (−1) j m s , where the mean-field staggered magnetization m s (t) has to be determined self-consistently. For developing an intuition it is worthwhile to imagine the dynamics of pseudo-spins in k-space by defining pseudo-spin operators σ z Note that these momentum-space pseudo-spins are different from the original spins on the chain. In pseudo-spin representation the staggered magnetization is given by the average x-projection per pseudo-spin, m s = 1 N π/2 k=−π/2 σ x k , and the mean-field Hamiltonian can be written as The Néel state as an initial condition corresponds to all pseudo-spins pointing in the x-direction at t = 0. Then they start to precess due to a Zeeman field that depends on the instantaneous average orientation of the x-projection of the spins. In these terms it is easy to understand the evolution of the staggered magnetization m s (t) for = 0. We simply have a collection of independent pseudo-spins subject to constant Zeeman fields J cos k in the z-direction. Because the field magnitude varies from spin to spin over a bandwidth, they precess at frequency 2J .
Since the band of precession frequencies is continuous, the spins gradually dephase, leading to the 1/ √ t decay (15) of the oscillation envelope of m s (t). The situation is more complicated in the case of = 0. Now there is also a field in the x-direction, which is the same for all spins but changes in time according to the instantaneous orientation of the spins. To lowest order in , i.e. setting m s (t) = m 0 (t) = J 0 (2J t)/2 in the Hamiltonian (30), the additional Zeeman field in the x-direction tilts the precession axis, giving rise to a smaller average x-projection of the spins and thus leading to a faster decay of m(t). The numerical results for the time evolution of the staggered magnetization according to (30) are shown in figure 10 for different values of . Finite leads to accelerated dephasing of the oscillations very much like in the unbiased calculations (section 3). However, the asymptotic law as extracted from the numerical solution by Hastings and Levitov [104] for 0 < | | 1 exhibits algebraically decaying oscillations with a t −2/3 envelope, This algebraic decay, as well as the two frequencies, which lead to a revival phenomenon [104], is in contradiction with the MPS calculations for the full Hamiltonian (1). For > 1, the staggered magnetization saturates to a nonzero value for t → ∞, which is presumably also an artifact of the mean-field approach not corroborated in the unbiased treatment. We conclude that the approach provides only a very rough picture of the order-parameter dynamics, which confirms the importance of collective effects, apparently not captured by the effective non-interacting mean-field Hamiltonian (29).

The XZ model-effective description for 1
In this section, we study the time evolution of the staggered magnetization m s (t) following a quench from the Néel state in the analytically treatable XZ model. This serves as a complementary analytical approach to the numerical investigation of the quench dynamics in the XXZ model in the regime of large anisotropies 1 and allows for a discussion of the long-time asymptotic behavior of m s (t). The XZ model is defined by the Hamiltonian At equilibrium the XZ model exhibits a quantum phase transition at = c = 2 which separates two gapped phases with antiferromagnetically ordered ground states in the z-direction for > c and in the x-direction for < c . It differs from the XXZ model (1) by terms violating the conservation of S z tot = j S z j , but has the advantage of being analytically diagonalizable. In the following, we will prove that the staggered magnetization in this model vanishes for all finite > c in the long-time limit after a quench from the Néel state and calculate the exact time evolution of m s (t) semi-analytically up to times J t ≈ 100, thus going beyond the time window accessible by the MPS calculation for the XXZ chain.
Using the Jordan-Wigner transformation for S x j and S z j and going over to momentum representation, the Hamiltonian of the XZ model (32) takes the form with a † k and a k denoting, respectively, creation and annihilation operators of spinless Jordan-Wigner fermions with quasi-momentum k. This Hamiltonian can be diagonalized by the Bogoliubov transformation, which maps (33) to a model of free fermions, with a dispersion ε k = J 1 + 2 /4 + cos 2k.
Since the initial Néel state is the ground state of (33) with = 0 → ∞, it is convenient to express the time-dependent (Heisenberg) Jordan-Wigner fermion operators a k (t) in terms of the Bogoliubov quasi-particle operators α 0 k which diagonalize the Hamiltonian (33) for the initial value = 0 , This reduces the computation of correlation functions at arbitrary time to the evaluation of ground-state expectation values. To calculate the time evolution of the staggered magnetization m s (t) following a quench in the XZ model, we define the two-spin correlation function, from which the square of the staggered magnetization is obtained by taking the infinite-range limit, In the fermionized picture of the XZ model, the two-spin correlator takes the form with A j = a † j + a j and B j = a † j − a j being the Majorana operators at lattice site j [66]. Using Wick's theorem, this correlation function can be expressed as a Pfaffian of pairwise contractions [44]. For the quench problem studied here, the explicit form of these contractions follows from (36) and is given by with φ k = θ k − θ 0 k (see also [67], where identical expressions have been derived for the transverse-field Ising model). We have taken the thermodynamic limit and converted the sums into integrals in the expressions above. In the limit t → ∞ for > c the evaluation of (38) reduces to the computation of a Toeplitz determinant, since the contractions (40) of the A j 's and B j 's among themselves vanish. Szegö's theorem can then be used to calculate the asymptotics of the Toeplitz determinant, yielding the result lim t→∞ C( , t) Thus, after a quench from the Néel state in the XZ model, the staggered magnetization vanishes for all finite > c at large times. At finite times, when the contractions (40) do not vanish, the Pfaffian representing the two-spin correlator (39) can be evaluated numerically at arbitrary times for a given distance. Due to the so-called light-cone effect [55,98], two spins at a distance are not causally connected at times smaller than ut < /2, since the correlation length of the initial Néel state is zero. Here u denotes the maximum (classical) speed of quasi-particles, which in the XZ model is given by u = max k (∂ k ε k ) = 2J . Exploiting this light-cone effect, the staggered magnetization can be calculated in terms of a finite-range correlation function, This method significantly reduces the computational effort at short times. We remark, however, that the light cone is not completely sharp in quantum-mechanical systems [55]. Nevertheless, for practical finite-precision calculations, the infinite-range limit of the two-spin correlator is reached for distances just a few lattice sites beyond the light cone. The results for the time evolution of the staggered magnetization following a quench from the Néel state in the XZ model are displayed in figure 11. As is the case for the XXZ chain, an explicit analytical expression for m s (t) in the XZ model can be derived for a quench to = 0, which is given by m s (t) = 0.5 cos 2 (J t). For < c , the numerical data for m s (t) at large times fits very well exponentially decaying oscillations of the form m s (t) ∝ e −t/τ 2 (cos 2 (ωt) − const). (44) In this regime, the behavior of m s (t) in the XZ model is qualitatively different from that in the XXZ model, as can be seen from the period of the magnetization oscillations. In the XZ model the period diverges at the critical point (see figure 12), whereas it becomes smaller upon approaching the isotropic point in the XXZ model (see figure 4). Furthermore, the critical point exactly marks the crossover between oscillatory and non-oscillatory behavior of m s (t) in the XZ model. For c , the staggered magnetization decays exponentially in the XZ model and shows no oscillations at large times. Interestingly, the numerical results for m s (t) in the XXZ and XZ models are almost indistinguishable at large anisotropies 1, as can be seen from figure 11. We have extracted the relaxation times from exponential fits to the numerical data, obtaining a clearly pronounced minimum right at the isotropic point (see figure 12). The relaxation time scales as τ 2 ∝ −1 for c and as τ 1 ∝ 2 for c .

Gapless theory-the Luttinger model
In the analysis of the XX limit (section 2), it became evident that, if the initial gap is sufficiently small, the non-oscillatory relaxation of the order-parameter dynamics is determined by low-energy modes, which motivates the application of the Luttinger model to a quench to the 0. gapless phase < 1 of the XXZ model. We note that our treatment is essentially equivalent to the one of Iucci and Cazalilla [58]. The inverse scenario, a quench starting from a gapless phase, has been studied in [53].
In appendix A, the Luttinger model, has been introduced as a low-energy effective theory for the XXZ chain in the easy-plane regime. The bosonized form of the staggered magnetization is given by m s ∼ cos(2φ) x=0 , where we have made use of the translational invariance. The remaining problem amounts to computing the time evolution of cos(2φ) , starting from a state where the field φ is initially pinned at 0 or π/2. We remark that this problem is essentially the dual of the dephasing problem studied in [112], and thus we expect an exponential decay of m s with a characteristic time scale τ ∼ J/(K s ). A convenient technique for solving this problem is the truncated Wigner method [63], which is exact for quadratic Hamiltonians such as (B.9). Using this approach, the time-dependent expectation value of the staggered magnetization can be written as a functional integral over the Wigner transform W (φ 0 ,φ 0 ) of the initial density matrix: Here, the functional δ-distribution ensures that one integrates only over solutions of the classical equations of motion and φ cl (x, t) denotes the classical solution of the 1D wave equation corresponding to the initial conditions φ 0 (x) andφ 0 (x). We have also used the fact that the operator cos(φ) is diagonal in the φ-representation. The solution φ cl (t) can be explicitly constructed using d'Alembert's formula. After switching to dual-field representation using K u∂ x θ =φ, we obtain Since in the initial state φ is pinned at φ 0 = 0, we factor out the φ-dependent part of the integral, obtaining where the brackets with the index 0 denote the expectation value taken with respect to the initial state. The rhs of equation (47) can be estimated within a semiclassical analysis, where the ground state of the Luttinger Hamiltonian (B.9) with an additional mass term ∼ s φ 2 is used as the initial state. This finally leads to where s again denotes the gap of the initial state. In contrast with the empirical rule (28) for the XXZ model, the Luttinger model, being a continuum theory, does not reproduce oscillations. The non-oscillatory relaxation in (48) is characterized by a relaxation time inversely proportional to the gap and to the Luttinger parameter, τ = 2 π K s , behavior identical to the conformal field theory result (8) and similarly observed in the numerical calculation for the quench in the XXZ model. However, the algebraic prefactor present in the case of the XXZ model (28) and the SDW initial state under the XX Hamiltonian (24) is not present in this treatment of the Luttinger model. Since the Luttinger model includes the XX limit at K = 1, we conclude that the missing algebraic prefactor is a shortcoming of the initial state, which has been approximated as the ground state of the Klein-Gordon model (B.16). More accurate results could provide a treatment using the sine-Gordon Hamiltonian, which as we shall see in the next section strongly complicates the problem.

Gapped theory-the sine-Gordon model
In this section, we analyze the quench in the sine-Gordon model (see appendix B), as a possible continuum approach to the quantum quench in the XXZ chain for > 1. In what follows, we use the boundary-state formalism as a convenient tool for describing the nonequilibrium problem [57]. In this formalism, the initial state, which is not the eigenstate of the Hamiltonian, can be thought of as a special superposition of pairs of eigenmodes of the quantum Hamiltonian with opposite momenta, which sums up into a squeezed state of eigenmodes [113]. Unlike for the gapless Luttinger-liquid theory, we cannot present a full solution of the dynamics. Possible directions to be followed in future are pointed out.
Since the sine-Gordon model has relativistic (Lorentz) invariance, we can exchange the time and space directions x ↔ t and consider the following boundary-in-time Hamiltonian (using more conventional notation x again for the imaginary time direction) where θ(x) is a theta-function and δ(x) takes care of the initial condition. In order to implement the Néel state as an initial condition we send 0 → ∞, which corresponds to the Dirichlet boundary (initial) condition. This boundary (initial) condition formulation can be reformulated in a boundary-state formalism of the boundary sine-Gordon model (bSG). The initial condition is expressed as a squeezed state of bulk degrees of freedom. We note that for non-interacting particles or Luttinger liquid this correspondence can be seen directly. Since for K < 1/2 there are only solitons and antisolitons in the spectrum (repulsive regime of the sine-Gordon model), we obtain the boundary state in the following form: Here A † a,b (θ ) is an operator of creation of the soliton (a) or antisoliton (b) and K ab D (θ) is a reflection matrix of the soliton-antisoliton pair corresponding to the Dirichlet boundary condition. The rapidity θ is related to the momentum P = M s sinh θ and energy E = M s cosh θ , where the soliton mass M s is given by [114] where we define ξ = β 2 /(1 − β 2 ). The evolution is trivial in the soliton basis, because the bulk Hamiltonian is diagonal in soliton-antisoliton operators, K ab D (θ, t) = K ab D (θ ) exp(2it M s cosh(θ)). In the boundary state formulation the evolution of the magnetization is equivalent to the computation of the following quantity: In general, the squeezed state represented by the boundary state |B(t) should be expanded as a series in powers of the reflection matrices. This produces multiple dynamical processes, which include solitons and antisolitons. Multi-particle expectation values of the operators, like cos(2φ), are called form factors. To compute the correlation functions in the massive theories at equilibrium, only a small number of lowest form-factor contributions is necessary. However, our evaluation of the lowest order contributions in our case provided results contradictory to the numerical simulations. The reason for this will be found in the spectral analysis of the next section, which hints that not only soliton-antisoliton form factors are important (which is the case for the spectral function of the sine-Gordon model for small energies), but also multiple processes, which include energies well above the spectral gap (soliton mass), are necessary to be considered. Technically, the problem of inclusion of multi-soliton form factors is rather difficult. The θ -integrals corresponding to evaluation of different multi-particle contributions become even more complicated because of the reflection matrices K ab D (θ). A possible alternative approach to this form-factor evaluation could be a resummation of the leading divergencies of the scattering processes in the presence of a boundary state. Since the ultraviolet (UV) energies make an important contribution, one can try to proceed by considering the logarithm of the one-or (two-) point function and to sum the leading contributions as proposed previously [115]- [118]. However, the complexity of the boundary reflection matrix does not allow us to realize this program. We hope to return to this problem in future.
In view of the high complexity of the boundary-state formalism, it may be worthwhile to establish phenomenological analogies between non-equilibrium dynamics and equilibrium dynamics of the sine-Gordon model. This would be useful since for calculating dynamical structure factors a powerful machinery has been developed over the last decades. Arguing that the initial state can be described by a thermal ensemble (with some effective temperature considered as a fitting parameter) instead of the boundary state, we can relate the dynamics of the magnetic order parameter to the two-point function, where the average is taken over some thermal ensemble characterized by temperature T eff . The operation of cos(2φ(0)) onto the thermal state is one possible way to introduce some magnetic order or, stated otherwise, to establish an analogy to equation (54); cos(2φ(t)) acting on the thermal state is a way to mimic a boundary-in-time state. We note that such an approach has been successfully applied for studying the dynamics of a non-local observable in quench in the quantum Ising chain [119]. The dynamics of the two-point function (55) is separated into two regimes: large temperatures T eff M s and low temperatures T eff M. It is known that for large energies (UV) massive models, like the sine-Gordon model, have the conformal field theory asymptotics. Therefore, in the large-temperature regime the behavior of the correlation functions should be the same as in the high-temperature limit of the corresponding conformal field theory. Hence, for T eff M s , the large-time asymptotics of the correlation function is given by an exponential decay This conformal field theory behavior is universal also for the gapless phase, where, at least in some regimes of weakly magnetized initial states, setting T eff = s this behavior is a good first approximation of the dynamics of the order parameter in the quench problem (see table 1). However, in the gapped phase, we cannot find a reasonable way to define T eff . For example, the temperatures corresponding to the Boltzmann ensembles used in the following section do not reproduce the numerical findings at all. In the other regime, T eff M s , the structure of the massive theory is important. In this case, the leading order behavior comes from the zero-momentum exchange processes and depends on the structure of scattering matrix S(0) in this limit. Resummation of the kinematical singularities leads again to the exponential decay for the two-point correlation function [120], in agreement with a quasi-classical formula from [121]. Implementing the results of [120] to our situation, we obtain where the proportionality coefficient depends on the power of M s . Such behavior, however, is in disagreement with our numerical findings, where in the limit of large M s (large ) we find a decay rate proportional to −2 . We conclude that although the sine-Gordon is a valuable candidate for describing the dynamics following a quantum quench in the XXZ model, the evaluation of the corresponding form factors is difficult and demands further efforts. A relation of the coherent dynamics of the order parameter to dynamical structure factors, circumventing this problem, is not straightforward to be established.

Spectral analysis
For a deeper understanding of the relaxation dynamics, it is useful to consider the problem in energy space. The idea is to associate the properties of the spectrum of the Hamiltonian to the dynamical phenomena observed in the simulation of the time evolution and to clarify the possibility of separating energy scales-a question that is especially important for improving analytical descriptions of the non-equilibrium dynamics.
Using the Lehmann representation, the time evolution of an operator O takes the form of a Fourier transform over the eigenlevels of the Hamiltonian, For a more convenient continuum description, we introduce the quenched probability distribution, which determines the properties of the stationary state at t → ∞ [84], [122]- [124] (the frequencies are shifted by E 0 -the ground-state energy of H). It can be compared to the thermal (Boltzmann) distribution of the grand canonical ensemble, where the temperature is set by the energy of the initial state, d ρ B ( ) = ψ 0 |H |ψ 0 . In general, it is known that the thermal distribution can deviate strongly from the quenched distribution [69,84], [122]- [124], which leads to the phenomenon of absence of thermalization.
Here the thermal distributions are used only as a reference and questions in connection with thermalization phenomena will not be investigated. While the quenched probability distribution, ρ ψ 0 , captures the effect of the initial state, the distribution of the expectation value, reflects the specific spectral properties of the given observable. ρ ψ 0 and O( , ) provide the contributions to the weighted expectation value, which, via the distribution function represents the dynamics of the observable in frequency space, The spectral properties of the XXZ chain prepared in the Néel state, |ψ 0 = |↑↓↑ · ·· ↓↑ , with the staggered magnetization as the observable, O = m s , are calculated by means of exact diagonalizations for small system sizes. Figure 13 displays the results for a chain of length N = 14. The small system size results in strongly peaked distributions, but we have made sure that the qualitative features we extract in the following analysis are stable against variations of the system size, both towards larger, N = 18, and smaller, N = 10, values. Quantitative information cannot be extracted from this simple analysis, but this might be possible when going to larger system sizes by means of more involved techniques, such as the Lanczos method [77].
In the non-interacting limit ( = 0), where the Hamiltonian has a free-fermion representation as discussed in section 2, all distributions are centrosymmetric about = E 0 . The quenched distribution ρ ψ 0 ( ) exhibits peaks at − E 0 = ±J , which are not present in the thermal distribution. From the discussion in section 2 and the finite size study (figures 13(a)-(c)), we can separate two contributions to W m s ( , ):  The isolated peak (v) at zero energy is irrelevant for the dynamics. It would correspond to a finite asymptotic value at t → ∞, which apparently vanishes in the thermodynamic limit.
In summary, on the basis of the analysis of the frequency distributions, we can now draw a qualitative crossover picture from the oscillatory to the non-oscillatory behavior as varies. Approaching the isotropic point from small values of , the band edges (i) are smeared out, leading to a decreasing relaxation time τ 2 . Starting from large values of , the low-frequency peak merges into the homogeneous distribution upon approaching the isotropic point. Hence, the characteristics of (i) and (iii) contributions, dominating at small (respectively large) values of , are both lost at intermediate values of , where the interplay of all energy scales apparently leads to a non-generic dynamical relaxation of the order parameter. We also note that, even in the regime 1, where the initial state is rather close to the ground state of the Hamiltonian, the relevant part of the spectrum located above the gap is a multi-particle continuum, difficult to treat analytically.

Conclusions and outlook
We have analyzed the dynamics of the staggered magnetization in quantum spin chains following a quantum quench, considering various antiferromagnetically ordered initial states by using a number of complementary numerical and analytical approaches. In the numerical MPS study we have essentially found three types of relaxation dynamics for the order parameter: (i) for highly ordered initial states ( 0 1) and sufficiently small anisotropy parameters of the Hamiltonian at t > 0, there are Rabi oscillations, which dephase exponentially in time away from the XX limit; (ii) for strong anisotropies (  1) there is an exponential decay and the relaxation time scales as τ ∝ 2 ; (iii) for initial states close to the phase transition, we found evidence for algebraic corrections to the exponential decay. There is a crossover phenomenon between oscillatory (small ) and non-oscillatory dynamics (large ), but no clear point of transition can be identified. Either both types of dynamics are superimposed, as is the case for 0 1, or both vanish in an extended transition regime, in which non-generic behavior is found (the case of 0 1).
We have shown that a precise description of the Rabi oscillations (i) is possible only when the full spectrum is taken into account. Therefore mean-field as well as low-energy approaches lead to incomplete results-an analytical treatment of this type of dynamics is feasible only by novel approaches. It has become clear that quasi-particles relevant for Rabi oscillations are located at the band edges. In contrast to equilibrium properties, where many-body effects can be incorporated into the relevant linearizable Fermi-level excitations within the bosonization formalism, it is not clear how to treat interactions in combination with the quadratic dispersion relation at the edges of the band. A treatment along the lines of previously developed concepts dealing with non-linear effects in dynamical phenomena [52,125] might be a possible solution to this problem.
While the exponential decay at large anisotropies (ii) appears to be rather generic behavior and is also reproduced in the exactly solvable XZ model, the algebraic prefactor in front of the exponential law (iii) is a more intricate phenomenon. In the XX limit, where we approximated the initial state by an SDW, such order-parameter dynamics are reproduced. However, standard field-theoretical approaches, such as conformal field theory or the description by Luttinger model adopted here, capture roughly the scaling of the corresponding relaxation time, but miss the prefactor. Our results suggest that this is a consequence of the phenomenological description of the initial state in terms of a simple massive theory (Klein-Gordon). A more elaborate treatment of the initial state using the sine-Gordon model could resolve this deficiency of the field-theoretical descriptions.
The sine-Gordon and the original XXZ model are integrable. Analyzing the quench by using the integrability of these models amounts to evaluation of form factors of particles with nontrivial statistics. While for two or three particles this problem can be solved [126,127], we have found that non-equilibrium dynamics requires the evaluation of higher order form factors-a yet unsolved and highly complex problem. A promising approach is to use the structure of the Bethe-ansatz solution in combination with a numerical algorithm [73].
Similarly to studies of other models [90,96] we find as a generic feature that relaxation times become small in the vicinity of the quantum critical point. However, the behavior is far too rich to be attributed to a generic dynamical phase transition [96]. In our analysis we have shown that effective descriptions by low-energy theories belonging to the universality class of the model cannot capture all relevant processes. A sort of dynamical phase transition occurs, however, in the mean-field description (a finite magnetization is found in the long-time limit for > 1), which treats interaction terms as on an infinite-dimensional lattice. This suggests that the existence of a sort of dynamical critical behavior may, very much like for equilibrium phase transitions, depend on dimensionality. A first step towards understanding the role of dimensionality is the study of non-equilibrium dynamics in infinite dimensions, for example by using dynamical mean-field theory [94,96]. How to treat coherent dynamics in a 2D or 3D system is still an open question.
Finally, we have to mention that the Heisenberg chain is a simplified model, well suited for numerical and analytical investigations, but not necessarily appropriate for a full description of experimental systems. Although in experiments with two-level atoms in optical lattices, behavior similar to our theoretical prediction is observed [19,21], the model has to be adjusted to provide an accurate description of the experimental results. For example, the effect of density fluctuations beyond the purely magnetic model needs to be investigated. Although the matrix product algorithm is an efficient method to study the relaxation dynamics, the numerically reachable times are fundamentally restricted by growing entanglement. The runaway time may become very small in models more involved than spin-1 2 chains, in particular when particle fluctuations need to be taken into account [86]. Recently, schemes have been proposed that allow us to go beyond what is possible within the conventional matrix product algorithms [128]- [131]. The main idea in these approaches is to calculate directly the dynamics of an observable rather than explicitly follow the evolution of the state. Another important aspect neglected here, but relevant in experiments, is temperature. How to efficiently include the effects of finite temperature within a time-dependent matrix product algorithm is still an unsolved problem [129,132]. degree of freedom σ represents typically a hyperfine state and can be identified with (pseudo) spin- 1 2 , σ =↑, ↓. The hopping integral t i jσ and the interaction parameters U σ σ depend on the geometry and depth of the lattice and can be expressed in terms of overlaps of the Wannier orbitals. If, for concreteness, we consider a periodic, spin-dependent lattice potential with an isotropic spacing a, with K = 2π a , the recoil energy needs to be much smaller than the lattice depth, E r =¯h V µσ , so that the atoms are in the lowest harmonic level and single-band description (A.1) is valid. In reality, the Gaussian shape of the laser beams introduces inhomogeneity in the lattice depth in addition to the harmonic trapping potential µ i . The superposition of polarized laser beams generates spin-dependent potentials [31,133]. In general, the orbitals extend only over short distances and the hopping can be restricted to nearest neighbors, where [134] µ↓ ) 2 is the spin-average potential in each direction, V µσ σ = V µσ , and a sσ σ is the s-wave scattering length between atoms of spin σ and σ .
In the case of strong on-site repulsion, t µσ U σ σ , and integer filling, the system is in the Mott phase, where occupation number fluctuations are essentially suppressed. In this case, the effective basis contains locally only singly occupied spin-up |↑ or -down |↓ states (for the sake of simplicity we choose n i↑ + n i↓ = 1, although higher occupation numbers are also possible). In this subspace, neglecting terms of order t 4 i j /U 3 σ σ , the Hubbard model (A.1) can be mapped onto a spin-1/2 Heisenberg model (XXZ model) [135]- [138], The superexchange interaction constants J i j z and J i j ⊥ are given by An analogous treatment can be carried out for the fermionic Hubbard model. In the resulting magnetic Hamiltonian (A.5), J ⊥ has the opposite sign compared to equation (A.7) and in the expression for J z the last two terms are absent since double occupancy is forbidden by the Fermi statistics.
For appropriately chosen lattice and interaction parameters, the anisotropy of the spin exchange, i j = J i j z /J i j ⊥ , is tunable to a large extent. For example, in the bosonic case with symmetric on-site repulsions, a ferromagnet with possible easy-axis anisotropy, i j = 1 2 ( t i j↑ t i j↓ + t i j↓ t i j↑ ) 1, is realized. In addition, as demonstrated recently [19], double-well potentials can be 35 used to change the sign of the exchange interactions from ferro-(J i j z < 0) to antiferromagnetic (J i j z > 0) [86]. Although the Heisenberg Hamiltonian (A.5) is a good first approximation to strongly interacting two-component Bose gases in optical lattices, we note that the measurements of Trotzky et al [19] clearly show the limitations of the purely magnetic picture. The strong repulsion leads to superexchange interaction, which reaches the order of currently realistic temperatures, J/k B ∼ 10 −9 K, and in non-equilibrium experiments the dynamics slow down, so that the effects of inhomogeneous laser beams become strong. For larger tunnelings, density fluctuations are important and introduce an additional higher frequency; excitations to higher Bloch bands may also become possible. We conclude that, although the experimental progress looks promising, further improvements in experimental setups are still needed in order to produce a clean realization of a quantum magnet.

Appendix B. Equilibrium properties of the XXZ model in one dimension
In a spatially anisotropic optical lattice, the Heisenberg chain, a paradigm in the theory of magnetism and strongly correlated systems in general, can be realized experimentally as proposed in the appendix A. Here we give an overview of the equilibrium phases of the Heisenberg chain, focusing on antiferromagnetic exchange interactions. At the same time, the important concepts and notations to be used in the ensuing discussion of the non-equilibrium problem are introduced.
The 1D spin-1/2 Heisenberg chain, is integrable-the eigenstates and an infinite number of conserved operators can be obtained using the Bethe ansatz [139]- [143]. A number of equilibrium properties can be exactly computed for the Bethe wave function-examples are the energy and momentum of lowlying states, or local observables such as the staggered magnetization [144,145]. For some specific cases, non-local properties can also be calculated analytically [146,147] or by means of a combination of the Bethe ansatz with numerical algorithms [73,148,149]. A simplified insight into the physics of the Heisenberg chain can be gained from a continuum description via the bosonization technique [108,150]. Here, results from both approaches, Bethe ansatz and bosonization, will be presented. The ground-state phase diagram of the XXZ model is presented in figure B.1. Without loss of generality the coupling J can be considered to be positive and the phases are simply characterized by . The long-range ordered antiferromagnetic phase for Ising-like anisotropies > 1 exhibits a spectral gap. In the easy-plane regime | | 1, a critical gapless phase is found. The phase for < −1 is ferromagnetically ordered.
A useful equivalent representation of (1) is a model of interacting spinless fermions, obtained from (1) by Jordan-Wigner transformation from spin operators to spinless fermion operators [65], In the case of a 1D Hamiltonian with nearest-neighbor interactions, particle statistics is irrelevant, and alternatively the fermions can also be replaced by hardcore bosons [108]. The fermionic picture is especially useful in the non-interacting case ( = 0, also known as the XX limit), where (B.2) is diagonal in Fourier space, In the case of zero magnetization, which is of interest here, the ground state is described by the half-filled Fermi sea, where |0 is the fermionic vacuum, c k |0 = 0. The (longitudinal) spin-spin correlation function, which characterizes magnetic ordering, can be calculated exactly in the XX limit [66], [151]- [156]. The result is a superposition of quasi-long-range ferromagnetic and antiferromagnetic correlations, decaying by a power law, For finite , the extraction of correlation functions from the Bethe-ansatz solution is highly nontrivial and only possible for some special cases (e.g. [147]). In order to obtain a continuum description of the Heisenberg chain, the spectrum of the non-interacting model is linearized at the Fermi points and the modes are separated into left and right movers, The cutoff is of the order of the bandwidth. Starting from (B.8), interactions can be included using the bosonization formalism [108]. At the renormalization-group fixed point, which captures the long-distance properties, the Luttinger model, provides an effective description for | | < 1. (x) and φ(x) are conjugate bosonic fields, . We note that the excitations described by the Luttinger model correspond to linearly dispersed spin waves with velocity u. The values of both, u and the Luttinger-liquid parameter K , can be derived from the Bethe ansatz [157], where α ∼ 1 . Here, the lattice spacing is set to one, so that the original sites are located at x = i, i = 1, . . . , N (N being the number of lattice sites).
For the quadratic Luttinger Hamiltonian (B.9), the correlation functions can be evaluated [108], The constants C 1 and C 2 have been calculated in [158]. Hence, in the whole planar phase (| | < 1), the correlations exhibit critical behavior and fall off algebraically.
A different situation has to be faced for 1. In the renormalization-group treatment, backscattering terms become important. For 1, the sine-Gordon Hamiltonian, is the effective model. At the isotropic point, = 1, K = 1/2, the cosine term is marginally relevant and leads to logarithmic corrections to the correlation function (B.12). For Isinglike anisotropies, > 1, 0 < K < 1/2, the cosine term is relevant-a spectral gap, s , opens and the phase φ becomes pinned at 0 or π/2. Hence, = 1 marks a phase transition to an antiferromagnetically ordered phase with a finite asymptotic value of the spin-spin correlations, (B.14) The two degenerate ground states, corresponding to φ = 0 or π/2, exhibit staggered magnetization, m s , of opposite signs, The spectral gap as well as the staggered magnetization are continuous in all derivatives in -the phase transition is of Berezinskii-Kosterlitz-Thouless type [159,160].  [161] or appendix C for the description of this method). However, the relation (B.19) is only valid for sufficiently large gaps.
In order to avoid dealing with the complicated structure of the antiferromagnetic states in the XXZ model, we introduce the SDW state, The coefficients of the wave function are related to the gap parameter s by The concept of MPS [162]- [165] as a generalization of valence-bond states [150,166,167] has been developed parallel in time with the DMRG algorithm [168,169]. DMRG established quickly as one of the most powerful numerical approaches for solving (quasi) 1D correlated many-body problems at equilibrium. Although DMRG was originally introduced as a realspace renormalization group, it can be understood as a variational optimization procedure in the space of matrix product states [171]. This identification of DMRG with MPS is especially useful for the implementation of the ideas of DMRG in the thermodynamic limit [161,171] and for time-dependent problems [78]- [80], [172]. In the following, we present a formulation of a DMRG-like algorithm, which is most suitable for both time-dependent and infinite-size calculations. The procedure is identical to the infinite-size time-evolving block decimation algorithm iTEBD [161], except that different matrices, introduced in the context of static DMRG [171], are used in order to improve the stability of the algorithm. Since neither the density matrix nor the renormalization group idea appears explicitly in this formulation, we refer to the algorithm as MPS or iMPS if the infinite-size limit is to be emphasized. Error analysis will be given for a specific case of a non-equilibrium problem in the thermodynamic limit, where we find that the behavior of the error can be considered identical to the case of the time-dependent DMRG for finite lattices [82].

C.1. MPSs
In order to construct an MPS, we consider a 1D lattice model where the Hilbert space can be separated into left and right subspaces L i and R i+1 (L i including i as the rightmost site, i + 1 being the leftmost site of R i+1 ). Generally, a wave function can be written as ) are orthonormal basis vectors of the space L i (R i+1 ). i is called the bond center matrix of bond i and constructs the density matrix of the left and right subsystems, respectively. If each site is described by a set of local basis vectors |s i of dimension d i (s i = 0, . . . , d i − 1), a state of the subspace can be expanded in terms of the local basis and the remaining subspace, The orthonormality of the basis imposes on A s i the left orthonormalization constraint, Equivalently, the state of the right subspace can be expanded by means of right orthonormalized matrices, An iterative expansion of an arbitrary state |ψ on a lattice of size N is possible, providing a matrix-product expression of the state, Similarly, introducing an iterative procedure, correlation functions can be calculated [171].

C.2. Schmidt decomposition
The preceding introduction of MPS is completely general and, if infinite-dimensional matrices are allowed, any state can be formally expressed in terms of a matrix product. A class of valence-bond states [150], [162]- [167] is indeed naturally formulated in terms of MPS. Also, product states are trivially represented as MPS. MPSs are, however, especially powerful in combination with an approximative numerical algorithm, providing the optimal reduced basis set for replacing a large or possibly infinite Hilbert space. The Schmidt decomposition, as described in the following, is the procedure that allows us to select the most relevant basis states.
If only a finite number of states m is supposed to be retained (in order to keep the dimension of the Hilbert space manageable for the computer), it can be shown [169] that a state |ψ approximates best the targeted state |ψ in the form (C.1), if it is defined as the Schmidt decomposition of rank m (site indices are omitted), The Schmidt coefficients, λ α , are the dominating eigenvalues of the singular value decomposition, satisfying α λ 2 α = 1. The discarded weight, w = α>m λ 2 α , (C.14) corresponds to the mismatch, ||ψ − |ψ | = w, introduced by this truncation procedure. The new basis is given in terms of the Schmidt states, In practice, it is useful to set only an upper bound for m (rather than fixing a definite value) and instead define a threshold such that only states for which λ 2 α are retained. The applicability of this truncation procedure to a physical state depends on the characteristics of the Schmidt values or, equivalently, the spectrum of the density matrix. The more slowly the values λ α decay, the larger must be the number of retained states. A generic expression for the spectrum of the density matrix has been obtained for a critical theory [173]; for practical purposes, it is however sufficient to consider the entanglement properties of the system to get the order of the necessary number of retained states. For instance, one can consider the entanglement entropy, The fact that in 1D equilibrium states the entanglement entropy exhibits logarithmic dependence on the typical length scale ξ of the state [174] (ξ corresponds to the correlation length or, at criticality, to the size of the system) guarantees an accurate description of a large class of wave functions using a finite number of states m ∝ ξ . Away from equilibrium, however, the entanglement generally grows linearly in time [99] and a potentially exponential growth of m with time restricts the applicability of an MPS to short times. For a wave function represented at bond i by the two-site center matrix ss i , the Schmidt decomposition reads as follows: replacing in equation (C.1) the contracted two-site center matrix with the bond center matrix, the Schmidt decomposition can be carried out as presented above. The matrices are updated retaining m Schmidt states (C.15), For the evaluation of correlation functions, additionally A i+1 or B i is needed. Although the effect of loss of orthogonality is spurious when applying the direct inverse (C.10) as in the original iTEBD algorithm [161], especially in the case of real-time evolution, a procedure for recovering both left and right orthonormalized representations is needed for stabilizing the algorithm. The left (right) rotation of the matrices does the job. Starting from a single-site center matrix, s i = i−1 B s i , the left orthonormalized matrix and the rotated center matrix can be extracted from the singular value decomposition (C.13) of the re-indexed matrix, An iterative application of this procedure moves the center matrix through the lattice and brings all matrices into left orthonormal form. An analogous left-moving iteration brings the matrices into the right orthonormalized form. In the periodic iMPS a problem arises when the rightmoving center matrix reaches the edge of the unit cell. R i does not in general coincide with the former i+1 and repeating the iterations through the unit cell further changes the MPS. There exist however schemes that solve this problem by introducing an additional transformation, after which the transfer operator R i −1 i+1 becomes equal to identity (see [171,175] for detailed descriptions).

C.3. The Suzuki-Trotter decomposition
In order to calculate the time evolution of an MPS, |ψ(t) = e −iH t |ψ 0 , it is suitable to approximate the evolution operator, e −iH t , by a Suzuki-Trotter decomposition. This is possible if the global operator H contains only nearest-neighbor bond terms H = i H i,i+1 (e.g. the Heisenberg chain with H i,i+1 = S i S i+1 ). H can then be decomposed into even and odd parts H = H 1 + H 2 , (C.20) The Suzuki-Trotter decomposition can be regarded as the first-order expansion of the evolution operator using the Baker-Hausdorff formula [176], e −iH t = (e −iH 2 δ e −iH 1 δ ) n + O(δ 2 n), nδ = t.

(C.21)
This approximation is improved in a second-order expansion, e −iH t = (e −iH 1 δ/2 e −iH 2 δ e −iH 1 δ/2 ) n + O(δ 3 n), (C. 22) or, if higher accuracy is desired, using third-or higher-order expansions [177]. where |ψ 0 is some random initial state. In order to get reliable results from this procedure, the Trotter slicing has to be reduced carefully during the imaginary-time evolution [161].

C.4. Update of an iMPS
We consider now the application of a single factor of (C.24) onto an MPS. For example, for the odd bond operator, U = e −iH 1,2 δ , we havẽ Since the inverse of 2 is required, a finite threshold is necessary to guarantee the stability of this operation.

C.5. Error analysis
As an application of the iMPS method to a non-equilibrium problem, we study the quench problem in the XXZ chain, |ψ(t) = e −iH t |ψ 0 , where |ψ 0 is the ground state of the XXZ Hamiltonian at a given value = 0 , and H is characterized by an anisotropy parameter . In this case, the local basis consists of a spin-up and a spin-down state {|↓ , |↑ }. This quench problem is analyzed in detail in section 3.
First, we study the case where |ψ 0 is the Néel state, which has a trivial iMPS representation with m = 1, A s 1 = B s 1 = δ s↑ , A s 2 = B s 2 = δ s↓ . Since the z-projection of the total spin (S z tot ) is conserved, the MPS can be resolved by this quantum number [171,178]. The resulting speedup is about an order of magnitude in comparison with a simulation that exploits no symmetry. is plotted for different values of the number of retained states m, the threshold and the Trotter slicing δ. The evolution of the error can be clearly divided into two regimes by introducing the runaway time t runaway : for t < t runaway , there is a small error that does not depend on the value of m. In this case, the error is dominated by the Trotter error, which grows at most linearly in time. However, for t > t runaway , the error starts growing nearly exponentially. The approximately logarithmic dependence of t runaway on m (figure C.1(a), inset) is in agreement with the linear growth of the entanglement entropy in the non-equilibrium problem-t runaway can be understood as the point where the chosen finite number of retained states is no longer sufficient to represent the entanglement in the state. We note, however, that a strict relation between entanglement entropy and t runaway cannot be rigorously established [82]. In order to reduce the Trotter error dominating at t < t runaway , one may choose smaller values of δ. The threshold has to be decreased as well. Otherwise, due to the increased number of updates, errors associated with the discarded weight at each step may accumulate. In figure C.1(b), we plot two cases with = 10 −15 and 10 −20 . In general, it is sufficient to reduce the threshold proportional to the Trotter slicing ∝ δt.
If the Trotter slicing is chosen so that the resulting error is of the order of the accuracy goal, t runaway , then sets the time window for the validity of the numerical results (in figure C.1, the accuracy goal in the absolute error is about 10 −6 ). We note that such behavior of the error in this time-dependent infinite-size MPS algorithm is identical to that of the finite-size DMRG algorithm [79]. If the exact solution is not known, t runaway can nevertheless be determined by comparing curves from calculations with slightly different m. t runaway is the point where the difference between them starts to grow significantly. Figure  where m s (t) m=1400 is the result for 1400 retained states, we find behavior identical to the exactly solvable case of the XX chain-the curves for m < 1400 overlap completely with the one for m = 1400 up to t < t runaway and a difference can only be seen for t > t runaway . For m = 1400, t runaway is estimated in figure C.2 by extrapolation of the values for m = 600, 800, 1000. Again, the accuracy of the results for t < t runaway is dominated by the Suzuki-Trotter error, which has to be estimated separately (here, it is of the order of 10 −7 ). In practice, it is not mandatory to abort the calculation at the runaway time-from the rough behavior of δm s (t) one can estimate that even for t 10, the absolute error of the curve for m = 1400 is still of the order of 10 −6 . The presented error analysis has been carried out for a local parameter in a specific setup. As far as non-equilibrium dynamics is concerned, this behavior is completely generic, although the runaway time and the Suzuki-Trotter error have to be determined for each case. Also, the error may depend on the observable under consideration-long-distance correlation functions may exhibit shorter runaway times than local observables. We would like to emphasize that error control, which imposes criteria on the wave function [179], is in general too strict, and the presented observable-based approach can considerably extend the accessible time window.