# Power-law decay of correlations after a global quench in the massive XXZ
chain

Flávia B. Ramos  Physics Department and Research Center OPTIMAS, Technische
Universität Kaiserslautern, 67663 Kaiserslautern, Germany  Andrew Urichuk
Department of Physics and Astronomy and Manitoba Quantum Institute, University
of Manitoba, Winnipeg R3T 2N2, Canada  Fakultät für Mathematik und
Naturwissenschaften, Bergische Universität Wuppertal, 42097 Wuppertal, Germany
Imke Schneider  Physics Department and Research Center OPTIMAS, Technische
Universität Kaiserslautern, 67663 Kaiserslautern, Germany  Jesko Sirker
Department of Physics and Astronomy and Manitoba Quantum Institute, University
of Manitoba, Winnipeg R3T 2N2, Canada

###### Abstract

We investigate the relaxation dynamics of equal-time correlations in the
antiferromagnetic phase of the XXZ spin-1/2 chain following a global quantum
quench of the anisotropy parameter. We focus, in particular, on the relaxation
dynamics starting from an initial Néel state. Using state-of-the-art density-
matrix renormalization group simulations, the exact solution of an effective
free-fermion model, and the quench-action approach within the thermodynamic
Bethe ansatz, we show that the late-time relaxation is characterized by a
power-law decay $\sim t^{-3/2}$ independent of anisotropy. This is in contrast
to the previously studied exponential decay of the antiferromagnetic order
parameter. Remarkably, the effective model describes the numerical data
extremely well even on a quantitative level if higher order corrections to the
leading asymptotic behavior are taken into account.

##  I Introduction

The high control available in experiments on ultracold atomic gases has opened
the door to investigate non-equilibrium dynamics in many fundamental quantum
models of Statistical Physics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Probing
relaxation and thermalization becomes possible in these systems because
decoherence and dissipation remain small for sufficiently long time scales.
The typical out-of-equilibrium scenario in cold atomic gases is the so-called
quantum quench where a control parameter is suddenly changed.

From a theoretical point of view, investigating how quantum many-body systems
equilibrate is an enormous challenge. Except for Gaussian problems, which can
either be solved analytically or at least simulated numerically for long
times, the exponential growth of the Hilbert space with the system size limits
exact diagonalizations to small systems. Since recurrence times are relatively
short in small systems, it is often very difficult to draw any conclusions
from such calculations about the thermodynamic limit. On the other hand,
methods based on a truncation of the Hilbert space such as the density-matrix
renormalization group (DMRG) [11, 12] or, more generally speaking, approaches
based on matrix-product states or operators, suffer from accuracy loss due to
the growth of the entanglement entropy (EE) with time. For global quenches in
clean one-dimensional systems, in particular, the EE grows linearly in time
which can be understood in terms of quasiparticles spreading in a light-cone-
like fashion [13].

The stationary state of a system is determined by its (quasi-)local conserved
charges. If the system is non-integrable then its subsystems are generally
expected to relax to a density matrix described by a Gibbs ensemble with an
effective temperature set by the energy density [14, 15, 16, 17]. In contrast,
the infinite number of local conserved charges in an integrable model prevent
thermalization. In this case, a relaxation to stationary values given by
generalized Gibbs ensembles (GGE) is expected [18, 19, 20]. Quench dynamics
has been studied in the context of effective free models, Luttinger liquids,
conformal field theories, and interacting systems, see for example Refs. [21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46].

While the relaxation dynamics after a quench is by now quite well understood
in Gaussian models, much less is known for interacting systems. The XXZ chain,
in particular, is well suited to serve as a simple model for the investigation
of interaction effects. Due to its integrability, methods based on the Bethe
ansatz can be used [32, 47, 48, 49, 50, 47]. In this context, one important
method is the quench-action approach (QA): If the form factors of the initial
state with the Bethe eigenstates are known, this method can be used to
calculate the stationary values of local observables at long times after the
quench. So far, however, the relaxation dynamics leading to equilibration has
remained out of reach of Bethe ansatz based approaches.

In this work, we investigate the time-dependent behavior of equal-time
correlators in the XXZ chain after an interaction quench within the gapped
antiferromagnetic phase. More specifically, we consider the quench from the
classical Néel state (the ground state of the XXZ chain for infinite
anisotropy) to a finite anisotropy. Using the state-of-art light-cone
renormalization group (LCRG) [12], we obtain numerical results in which the
long-time regime of local correlation functions can be studied. On the
analytical side, inspired by Barmettler _et al._ [25, 26], we consider the
large-anisotropy limit, in which the time evolution of the system can be
related to the dynamics governed by the XZ chain. In this effective framework,
we obtain exact formulas for the long-time behavior of short-distance
correlators. In contrast to the non-oscillatory exponential decay of the order
parameter [25, 26], we show based on a saddle-point analysis that the
correlations decay asymptotically in a power-law fashion with interaction-
independent exponents. In addition, we use the QA method to compute the
stationary values of the time-dependent correlations in the XXZ model.
Finally, we define an asymptotic formula, based on the analytical results for
the effective model, to fit the LCRG data. A remarkable, quantitative
agreement between the XZ results with renormalized parameters and the
numerical results for the XXZ chain is reported.

This paper is organized as follows: In Sec. II, we introduce the model, quench
protocol, and methods we use to investigate the problem of interest. We devote
Sec. III to the discussion of the large-anisotropy limit. Here, we map the
interacting problem to a spinless free-fermion model and derive its exact
solution as well as the asymptotic behavior of short-distance correlators. In
Sec. IV, we present our LCRG results for three different correlators and
discuss them in terms of the effective theory. Finally, we provide concluding
remarks in Sec. V.

##  II Model and methods

We consider the isolated XXZ chain described by the Hamiltonian

$H(\Delta)=J\sum_{j}\left(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}+\Delta\sigma_{j}^{z}\sigma^{z}_{j+1}\right),$
(1)

where $\sigma^{x,y,z}_{j}$ are the Pauli matrices acting at site $j$. $J>0$ is
the exchange coupling and $\Delta$ is the anisotropy parameter. This model
exhibits a U(1) symmetry corresponding to invariance under rotation by an
arbitrary angle around the $z-$axis. Furthermore, the XXZ chain is exactly
solvable by Bethe ansatz and its ground-state phase diagram is characterized
by a gapless phase with quasi-long-range order for $|\Delta|\leq 1$ and a
gapped antiferromagnetic (ferromagnetic) phase for $\Delta>1$ $(\Delta<-1)$
[51]. The limits $\Delta\rightarrow\infty$ and $\Delta\rightarrow-\infty$
correspond to the classical Ising antiferromagnet and ferromagnet,
respectively.

Here, we will focus on the relaxation dynamics after a global quantum quench
of the anisotropy in the antiferromagnetic phase of the XXZ chain. The quench
protocol is as follows: The system is prepared in the Néel state defined as

$|\psi_{N}\rangle=|\uparrow\downarrow\uparrow\downarrow\cdots\uparrow\downarrow\rangle,$
(2)

which is a ground state of $H(\infty)$. At time $t=0$, the anisotropy
parameter is suddenly quenched to a finite $\Delta$, so that
$|\psi_{N}\rangle$ is no longer an eigenstate of $H(\Delta)$. The unitary time
evolution of the system is governed by the time-dependent Schrödinger equation

$|\psi(t)\rangle=e^{-iH(\Delta)t}|\psi_{N}\rangle.$ (3)

We focus, in particular, on the equal-time correlations along the
$z$-direction at time $t$ after the quench which are given by

$C(\ell,t)=\langle\psi_{N}|\sigma_{0}^{z}(t)\sigma_{\ell}^{z}(t)|\psi_{N}\rangle.$
(4)

The number of excitations (spin flips) on top of the Néel state depends on the
anisotropy parameter $\Delta$ during the time evolution. For $\Delta\gg 1$,
the Néel state is close to the ground state of $H(\Delta)$ so that the density
of excitations during the time evolution is expected to be low. If only few
spin flips (kinks) are present, the time evolution can be approximated by that
of an effective XZ chain [52, 53, 54]. This effective model can be mapped onto
a Gaussian fermionic model and the long-time behavior of $C(\ell,t)$ can be
extracted from a saddle-point analysis (see Sec. III). In addition, the
stationary values for the full XXZ chain (1) can be calculated using the QA
approach. Some stationary values of short distance correlators at large
$\Delta$ are already available from Ref. [48].

We will numerically calculate the time dependence of local operators in the
thermodynamic limit, e.g. $C(\ell,t)$, using the LCRG method. In contrast to
the infinite time-evolving block decimation, LCRG does not require
translational invariance. The method relies on a Suzuki-Trotter (ST)
decomposition of the time evolution operator and a systematic truncation of
the Hilbert space that is iteratively carried out using a reduced density-
matrix [12]. In comparison to the time-dependent DMRG, the major advantage of
LCRG lies in the computational speed-up achieved due to the light-cone
structure of transfer matrices, which allows a reduction of bond operators in
the ST decomposition.

As already mentioned in the introduction, the entanglement imposes limitations
on the simulation times achievable by DMRG-based methods. Keeping up to 65,000
states to represent the truncated Hilbert space and using a second-order ST
decomposition, we can reach times $Jt\lesssim 5$ while keeping the truncation
error below $10^{-9}$. It is worth to keep in mind that if one wants to
compare this simulation time with previous studies that we are using Pauli
matrices in the definition of the Hamiltonian instead of spin-$1/2$ operators.
If we rewrite the Hamiltonian (1) in terms of spin operators, then our
simulation times correspond to $Jt\lesssim 20$, i.e. in this case
$t\rightarrow 4t$. The numerical results are more sensitive to the finite
Trotter step $\delta t$ the deeper one goes into the massive phase. This is
because the amplitude of oscillations in the correlators becomes smaller with
increasing $\Delta$ while the oscillation frequency increases. By analyzing
the convergence of our numerical results for different values of $\delta t$,
we find that accurate results are obtained for $\delta t=0.005$ $(\delta
t=0.001)$ if $\Delta\leq 12$ $(\Delta>12)$. The Trotter error for the
decomposition used is always of order $(\delta t)^{2}$. Additional numerical
details are reported in Sec. IV.

Let us now briefly comment on how stationary values are obtained from the QA
approach. The overlap between the initial state, here the Néel state, and
Bethe states can be determined [33, 55] and inserted into the quench action.
Once these overlaps are known, the minimization of the quench action yields
integral relations—known as the overlap thermodynamic Bethe ansatz
equations—which determine the saddle-point string densities [48, 56] as

$0=\frac{\delta
S_{\text{QA}}[\vec{\rho}]}{\delta\rho_{n}}\bigg{|}_{\vec{\rho}=\vec{\rho}^{\text{sp}}}\quad\text{for}\quad
n\geq 1.$ (5)

These saddle-point solutions are combined with quantum transfer matrix
formulas to determine the long-time limit as

$\lim_{t\to\infty}C(\ell,t)=\langle\rho^{\text{sp}}|\sigma_{0}^{z}\sigma_{\ell}^{z}|\rho^{\text{sp}}\rangle.$
(6)

In order to have a guide to analyze our numerical results, we start with the
analytical investigation of the strong-anisotropy limit.

##  III Strong-anisotropy limit

Using the raising and lowering operators
$\sigma^{\pm}=(\sigma^{x}\pm\text{i}\sigma^{y})/2$ we can write
$\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}=2(\sigma^{+}_{j}\sigma^{-}_{j+1}+\sigma^{-}_{j}\sigma^{+}_{j+1})$
while
$2\sigma^{x}_{j}\sigma^{x}_{j+1}=2(\sigma^{+}_{j}\sigma^{-}_{j+1}+\sigma^{-}_{j}\sigma^{+}_{j+1})+2(\sigma^{+}_{j}\sigma^{+}_{j+1}+\sigma^{-}_{j}\sigma^{-}_{j+1})$.
For $\Delta\gg 1$, states which have two adjacent parallel spins will
contribute little to the quench dynamics when starting from the Néel state so
that replacing $\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}$
by $2\sigma^{x}_{j}\sigma^{x}_{j+1}$ is a reasonable approximation. Thereby,
the XXZ chain effectively becomes the XZ model, whose Hamiltonian is

$H_{{XZ}}=J\sum_{j}\left(2\sigma^{x}_{j}\sigma^{x}_{j+1}+\Delta\sigma^{z}_{j}\sigma^{z}_{j+1}\right).$
(7)

We note that this approximation leads to an explicit U(1)-symmetry breaking
due to the addition of the terms $\sigma^{+}_{j}\sigma^{+}_{j+1}$ and
$\sigma^{-}_{j}\sigma^{-}_{j+1}$ which do not conserve the total
magnetization.

The XZ chain can be mapped to the exactly solvable XY chain by a 90∘-degree
rotation around the $x$-axis. After the rotation, the Hamiltonian in Eq. (7)
can be expressed in terms of the ladder operators as

$\displaystyle
H_{{XZ}}=J\sum_{j}\left[(2+\Delta)(\sigma^{+}_{j}\sigma^{-}_{j+1}+\sigma^{-}_{j}\sigma^{+}_{j+1})\right.$
$\displaystyle\left.+(2-\Delta)(\sigma^{+}_{j}\sigma^{+}_{j+1}+\sigma^{-}_{j}\sigma^{-}_{j+1})\right].$
(8)

The system described by Eq. (8) corresponds to a spinless free-fermion model
and its exact solution is obtained by means of the Jordan-Wigner
transformation

$\sigma^{+}_{j}\rightarrow
c^{\dagger}_{j}e^{i\pi\phi_{j}},\quad\sigma^{-}_{j}\rightarrow
c_{j}e^{-i\pi\phi_{j}},$ (9)

where $c^{\dagger}_{j}$ and $c_{j}$ are the fermionic creation and
annihilation operators acting at site $j$ and
$\phi_{j}=\sum_{\ell=0}^{j-1}c^{\dagger}_{\ell}c_{\ell}$ is the Jordan-Wigner
string. In Fourier space, the fermionic Hamiltonian can then be diagonalized
by performing a Bogoliubov transformation that maps the operators
$c_{k}=\frac{1}{\sqrt{L}}\sum_{j}e^{-i\pi j}c_{j}$, where $L$ is the system
size, to a new set of fermionic operators $a_{k}$. The corresponding
transformation $U_{\Delta}$ is defined by

$\displaystyle\left(\begin{array}[]{c}c_{-k}\\\
c_{k}^{\dagger}\end{array}\right)=U_{\Delta}\left(\begin{array}[]{c}a_{-k}\\\
a_{k}^{\dagger}\end{array}\right),$ (14) $\displaystyle
U_{\Delta}=\left(\begin{array}[]{cc}\cos\theta_{k}&-i\sin\theta_{k}\\\
-i\sin\theta_{k}&\cos\theta_{k}\end{array}\right),$ (17)

with $\tan{\left(2\theta_{k}\right)=\frac{2-\Delta}{2+\Delta}}\tan(k)$. The
Hamiltonian in the new operators then reads

$H_{XZ}=\frac{1}{2}\sum_{k}\left(\begin{array}[]{cc}a_{k}^{\dagger}&a_{-k}\end{array}\right)\left(\begin{array}[]{cc}\varepsilon_{k}&0\\\
0&-\varepsilon_{k}\end{array}\right)\left(\begin{array}[]{c}a_{k}\\\
a_{-k}^{\dagger}\end{array}\right),$ (18)

where $\varepsilon_{k}=4J\sqrt{1+\Delta^{2}/4+\Delta\cos(2k)}$ is the
dispersion relation. Hence, the ground-state phase diagram of the system is
characterized by a critical point at $\Delta=2$ separating two
antiferromagnetic phases with alignment along the $x-$ and $z-$direction for
$0<\Delta<2$ and $\Delta\geq 2$, respectively. Note that the classical Néel
state, $|\psi_{N}\rangle$, is the ground state of $H_{XZ}$ for
$\Delta=\infty$.

Let us now turn our attention to the asymptotic behavior of $C(\ell,t)$. Using
the transformation in Eq. (9), $C(\ell,t)$ is written as

$C(\ell,t)=-\langle(c^{\dagger}_{0}-c_{0})(c^{\dagger}_{\ell}e^{i\pi\phi_{\ell}}-e^{-i\pi\phi_{\ell}}c_{\ell})\rangle.$
(19)

Defining $A_{j}=(c^{\dagger}_{j}+c_{j})$ and $B_{j}=(c^{\dagger}_{j}-c_{j})$,
and observing that $e^{i\pi c^{\dagger}_{j}c_{j}}=-B_{j}A_{j}$, we have

$C(\ell,t)=(-1)^{\ell}\langle A_{0}B_{1}A_{1}\cdots
A_{\ell-1}B_{\ell}\rangle.$ (20)

The $\ell$-point correlation can then be decomposed into a combination of
pairwise contractions via Wick’s theorem [57]. Thereby, the problem is reduced
to calculating the contractions $\langle A_{m}A_{n}\rangle$, $\langle
A_{m}B_{n}\rangle$, $\langle B_{m}A_{n}\rangle$, and $\langle
B_{m}B_{n}\rangle$. To do so, it is convenient to write the time-dependent
fermionic operators in the basis which diagonalizes $H_{XZ}$ for
$\Delta=\infty$

$\left(\begin{array}[]{c}c_{-k}(t)\\\
c_{k}^{\dagger}(t)\end{array}\right)=U_{\Delta}\left(\begin{array}[]{cc}e^{-i\varepsilon_{k}t}&0\\\
0&e^{i\varepsilon_{k}t}\end{array}\right)U_{\Delta}^{-1}U_{\infty}\left(\begin{array}[]{c}a_{-k}^{\infty}\\\
a_{k}^{\infty\dagger}\end{array}\right),$ (21)

so that $a^{\infty}_{k}|\psi_{N}\rangle=0$. For $m\neq n$, the aforementioned
contractions can then be written as

$\displaystyle\langle A_{m}A_{n}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{-ik(m-n)}\left[1-\sin{(2\varepsilon_{k}t)}\sin{(2\tilde{\theta}_{k})}\right],$
$\displaystyle\langle A_{m}B_{n}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{-i[k(m-n)+2\theta_{k}]}\left[\right.\cos(2\tilde{\theta}_{k})$
(22)
$\displaystyle\qquad+\text{i}\cos{(2\varepsilon_{k}t)\sin(2\tilde{\theta}_{k})}\left.\right],$

where $\tilde{\theta}_{k}=\theta_{k}-\theta_{k}^{\infty}$ with
$\theta_{k}^{\infty}=-k/2$. Besides, we have $\langle
B_{m}B_{n}\rangle=\langle A_{m}A_{n}\rangle$ and $\langle
A_{m}B_{n}\rangle=-\langle B_{n}A_{m}\rangle$. More details about the
calculations discussed in this section can be found in App. A. For simplicity,
we will set $J=1$ as the energy scale henceforth.

Now, we focus on the distance-1 correlator,

$\displaystyle
C(1,t)=-\frac{1}{2\pi}\int_{-\pi}^{\pi}dk\,e^{i(2\theta_{k}+k)}\left[\cos(2\theta_{k}+k)\right.$
$\displaystyle-\text{i}\left.\cos(2\varepsilon_{k}t)\sin(2\theta_{k}+k)\right].$
(23)

We use the saddle-point method to determine the asymptotics of $C(1,t)$. By
deforming the complex contour in the vicinity of the roots of
$\varepsilon_{k}^{\prime}=0$, this method provides the leading contributions
to the long-time regime. In particular, by expanding $\varepsilon_{k}$ up to
fourth order around the inflection points, we obtain

$\displaystyle C(1,t)$ $\displaystyle\approx$
$\displaystyle-\left(1-\frac{2}{\Delta^{2}}\right)-\frac{\left(1+\frac{2}{\Delta}\right)^{1/2}}{32\sqrt{\pi}(1+\Delta^{2}/4+\Delta)}\left\\{\left(1+\frac{2}{\Delta}\right)\frac{\sin\left[\tilde{\omega}_{1}(\Delta)t-\frac{\pi}{4}\right]}{t^{3/2}}+\frac{3(1+\Delta^{2}/4-\Delta/2)}{16\Delta^{2}}\right.$
(24) $\displaystyle\times$
$\displaystyle\left.\frac{\cos\left[\tilde{\omega}_{1}(\Delta)t-\frac{\pi}{4}\right]}{t^{5/2}}\right\\}-\frac{\left(1-\frac{2}{\Delta}\right)^{1/2}}{32\sqrt{\pi}(1+\Delta^{2}/4-\Delta)}\left\\{\left(1-\frac{2}{\Delta}\right)\frac{\sin\left[\tilde{\omega}_{2}(\Delta)t-\frac{3\pi}{4}\right]}{t^{3/2}}\right.$
$\displaystyle+$
$\displaystyle\frac{3(1+\Delta^{2}/4+\Delta/2)}{16\Delta^{2}}\left.\frac{\cos\left[\tilde{\omega}_{2}(\Delta)t-\frac{3\pi}{4}\right]}{t^{5/2}}\right\\}+O(t^{-7/2}),$

where $\tilde{\omega}_{1}(\Delta)=8(\Delta/2+1)$ and
$\tilde{\omega}_{2}(\Delta)=8(\Delta/2-1)$. The first term in Eq. (24) is the
stationary value of $C(1,t)$. For $\Delta\rightarrow\infty$, this result is
consistent with the energy density of the classical Néel state. Moreover, we
observe that the time-dependent part contains two oscillation frequencies
proportional to the anisotropy and two exponents, independent of anisotropy,
that define the dominant power-law decay of the correlation function.

Following the same rationale, one can obtain distance $\ell=2$ and $\ell=3$
correlators. It turns out that $C(2,t)$ is particularly simple because only
$\langle A_{0}B_{1}\rangle\langle A_{1}B_{2}\rangle$ does not vanish after
applying Wick’s theorem. In this case, $C(2,t)=\left(C(1,t)\right)^{2}$. On
the other hand, for $C(3,t)$, two non-vanishing terms remain

$C(3,t)=-\langle A_{0}B_{1}\rangle^{3}+\langle A_{0}B_{1}\rangle\langle
A_{2}B_{1}\rangle\langle A_{0}B_{3}\rangle.$ (25)

Up to fourth order in the saddle-point analysis, each contraction in Eq. (25)
contributes with one time-independent term, two terms with distinct
frequencies that decay with power-law exponent $3/2$ and two other with
exponent $5/2$ (see App. A). In particular, the saturation value $C(3,\infty)$
is given by

$C(3,\infty)\approx-\left(1-\frac{2}{\Delta^{2}}\right)^{3}-\frac{\left(\Delta^{2}-2\right)(\Delta^{2}-4)}{\Delta^{6}}.$
(26)

It is worth mentioning that in contrast to short-distance correlators, the
relaxation of the order parameter is difficult to calculate because it is
nonlocal in the fermionic basis. However, a numerical evaluation for large
system sizes is possible [25, 26] and Barmettler _et al._ find that the long-
time behavior of the staggered magnetization for $\Delta\geq 2$ shows a non-
oscillatory exponential decay with a relaxation time scale proportional to
$\Delta^{2}$. It is also worth noting that the time-dependent behavior of
correlators for a quench within the ferromagnetic phase of the transverse-
field Ising chain was studied in Ref. [31].

##  IV Numerical results

Figure 1: Log-log scale plot of the slope $\alpha$ extracted from the linear
growth of the EE as a function of $\Delta$. The inset shows typical results of
the time dependence of $S(t)$ for $\Delta=4$ and $12$.

Figure 2: Correlator $C(1,t)$ as a function of time for different Trotter
steps and anisotropies: (a) $\Delta=5$ and (b) $\Delta=16$. The dashed lines
indicate the stationary values obtained from the QA approach.

In this section, we present numerical results for the quench dynamics of the
XXZ chain obtained by the LCRG. As discussed in the previous section, in the
massive regime the quench dynamics is approximately governed by the effective
XZ model for which the time-dependent correlations can be exactly calculated
via a fermionization. Furthermore, the saturation values of $C(\ell,t)$ can be
directly computed for the XXZ model using the QA method. We are therefore in a
position where we can provide a systematic analysis of the LCRG data for the
quench dynamics.

We start discussing the two main sources of errors in our numerical results:
the truncation error and the Trotterization of the time evolution operator. To
keep the largest truncation error of order $10^{-9}$ or below at late times,
we include up to $65,000$ states in the truncated Hilbert space. We observe
that longer simulation times can be reached for larger anisotropy $\Delta$.
This trend can be understood by considering the bipartite entanglement entropy
$S(t)$ which grows linearly in time, $S=\alpha t$. As can be seen in Fig. 1,
the slope $\alpha$ is decreasing approximately in a power-law fashion. The
typical time dependence of $S(t)$ is shown in the inset of Fig. 1 for two
distinct values of $\Delta$. On the other hand, the oscillation frequency
increases with $\Delta$ while the amplitude of the oscillation decreases.
Thus, the time step (Trotterization) plays an important role for the accuracy
of the simulations at strong anisotropy. In Fig. 2, we show how our numerical
results depend on $\delta t$. Stationary values obtained from the QA are also
indicated in the figure. We use these exact values to determine a time step
$\delta t$ which leads to errors that are small compared to the amplitude of
the oscillations.

Figure 3: Equal-time correlations for distances $\ell=1,2$ and (a) $\Delta=4$,
(b) $\Delta=15$. The symbols are LCRG results and the red solid lines are fits
to our data using Eq. (27). To display the results for both distances in the
same plot, we added an offset to $C(2,t)$ as indicated. (c) Oscillatory part
of the distance-1 correlator on a log-log scale for $\Delta=4$ and 15. The
dashed red lines are power laws with exponent $-3/2$. Figure 4: Stationary
values of $C(1,t)$ obtained by fitting the LCRG data, using QA, and the exact
solution of the XZ chain. The inset shows estimates of the oscillation
frequencies $\omega_{\alpha}$ as a function of $\Delta$ for $C(1,t)$. The
analytical results $\omega_{1,3}=8(\Delta/2+1)$ and
$\omega_{2,4}=8(\Delta/2-1)$ for the strong-anisotropy regime are also shown
in the plot.

Let us now discuss the temporal behavior of $C(\ell,t)$. In Fig. 3, we show
the time-dependent correlations for $\Delta=4$ and $\Delta=15$. Using the
analytical results for $C_{\infty}=\lim_{t\to\infty}C(\ell,t)$ from the QA and
the asymptotics for the XZ model discussed in Sec. III, we define a function
to fit distance-1 and distance-2 correlators

$\displaystyle|C(\ell,t)|$ $\displaystyle=$
$\displaystyle\left\\{C_{\infty}+\frac{A_{1}\sin(\omega_{1}t+\phi_{1})+A_{2}\sin(\omega_{2}t+\phi_{2})}{t^{3/2}}\right.$
(27) $\displaystyle+$
$\displaystyle\left.\frac{A_{3}\cos(\omega_{3}t+\phi_{3})+A_{4}\cos(\omega_{4}t+\phi_{4})}{t^{5/2}}\right\\}^{\ell}.$

We keep the exponents of the leading and subleading terms fixed to reduce the
number of free parameters in the fitting equation and because we do not expect
them to change as compared to the XZ model. We notice that close estimates of
these exponents are obtained if they are keptf2 as free parameters. Note that
Eq. (27) perfectly fits the time dependence of the LCRG results. Similar high-
quality fits can be obtained for all $\Delta\gtrsim 3$. A clear power-law
decay is observed in Fig. 3(c). The red dashed lines in the figure indicate
that the leading time-dependent contribution decays with exponent 3/2. In Fig.
4, we summarize the estimates for $\omega_{1-4}$ and $C_{\infty}$ obtained
from such fits as a function of $\Delta$ for the distance-1 correlator. A
comparison with the exact solution of the XZ chain and the QA values for
$C_{\infty}$ is also shown. Overall, we observe an excellent agreement between
the fitted oscillation frequencies for the XXZ model and those for the exactly
solvable XZ chain. The fitted stationary values $C_{\infty}$ are slightly
shifted as compared to those for the XZ model but agree very well with the
exact QA results. Very similar results are obtained for the case of $\ell=2$.
It is worth stressing that the remarkable agreement for all shown $\Delta$
values is surprising given that the effective theory might only be expected to
hold at very large $\Delta$. To summarize, we observe an excellent agreement
with the predictions of the XZ model with a small, $\Delta$-dependent shift in
the stationary value—fully covered and understood by the QA approach—as well
as a renormalization of the oscillation amplitudes. The oscillation
frequencies, on the other hand, remain essentially unchanged.

The analysis of the regime $1<\Delta<3$ is challenging. First, in the regime
$\Delta>2$, where the XZ model remains ordered along the $z-$direction, higher
order corrections to the asymptotics derived in Sec. III become more important
and longer times are needed before they can be neglected. This can be seen by
comparing with a numerical integration of Eq. (23). Reaching much longer time
scales is, however, not possible with current algorithms and computing
resources. For $1<\Delta<2$, there might also be additional contributions to
the asymptotics which are not captured by the XZ model.

Finally, let us also discuss the distance-3 correlator. In this case, due to
the combination of different Wick contractions, many free parameters would be
required to define a fit function and achieving convergence of the fit becomes
more difficult. Instead, we concentrate on large $\Delta$ where we can
directly evaluate Eq. (25) by a numerical integration. In Fig. 5, we show a
typical result for such a correlator for $\Delta=15$ and compare it directly
to the XZ model and the stationary value obtained by the QA. In the inset of
Fig. 5, we show that the stationary value obtained for the XZ model, Eq. (26),
agrees well with the QA result.

Figure 5: Equal-time correlation $C(3,t)$ for $\Delta=15$. The filled circles
are LCRG results and the solid red line is the result for the XZ model, Eq.
(25) [see also Eqs. (38)-(41) in App. A]. The black dashed line is the exact
stationary value obtained from QA. The inset compares the QA stationary values
for the XXZ chain and the stationary value for the XZ chain, Eq. (26).

##  V Conclusions

We have studied the relaxation dynamics after a global quantum quench in the
massive antiferromagnetic phase of the XXZ chain using the Néel state as
initial state. By adding terms that explicitly violate the conservation of the
total magnetization in the XXZ model, we can obtain a XZ model which is
exactly solvable via a Jordan-Wigner transformation. We have argued that at
large anisotropy $\Delta$, the terms violating the conservation of the total
magnetization have little effect on the dynamics. The XZ model is thus a
useful effective model for which the long-time asymptotics of correlation
functions can be obtained in analytic form via a saddle-point integration. In
addition, we have used the quench action approach to obtain exact results for
the stationary values of correlation functions in the XXZ chain at long times
after the quench.

These analytical results then formed the basis for our analysis of numerical
data for the quench dynamics in the XXZ chain obtained by the light-cone
renormalization group, a variant of the density-matrix renormalization group.
Quite remarkably, we find that the XZ asymptotics describes the numerical data
extremely well for $\Delta>3$ if a shift in the stationary value—exactly known
from the QA—and a small renormalization of the amplitudes is taken into
account. We find, in particular, that the power-law exponents of the decay
towards the stationary value as well as the oscillation frequencies remain
unchanged.

In summary, we find that the behavior of correlators
$\langle\Psi_{N}|\sigma^{z}_{0}(t)\sigma^{z}_{\ell}(t)|\Psi_{N}\rangle$ at
long times $t$ after the quench is described by an oscillatory behavior with
oscillation amplitudes which decay as a power law with interaction-independent
exponents and a leading power $t^{-3/2}$. We note that in contrast, the
stationary values, oscillation frequencies, and amplitudes of the oscillations
all do depend on the anisotropy $\Delta$ of the quench Hamiltonian.

For the future, it would be interesting to extend this study to quenches where
the anisotropy of the quench Hamiltonian is within the critical phase of the
XXZ chain. If we start from the Néel state this will be, however, never a
small quench and there seems to be no effective model which could potentially
aid in the analysis of the numerical data. Furthermore, the entanglement
entropy in the numerical simulations does grow much more rapidly in this case,
further limiting the obtainable simulation times.

###### Acknowledgements.

This work was supported by the DFG via the SFB/Transregio 185 (A4) and the
research unit FOR 2316. The authors thank Compute Canada and the high
performance cluster “Elwetritsch” for providing computational resources. J.S.
acknowledges support by the Natural Sciences and Engineering Research Council
(NSERC) of Canada and thanks F. Essler and J.-S. Caux for discussions.

## References

  * Kinoshita  _et al._ [2004] T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 1125 (2004).
  * Paredes  _et al._ [2004] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, Nature 429, 277 (2004).
  * Bloch [2005] I. Bloch, Nature Physics 1, 23 (2005).
  * Kinoshita  _et al._ [2006] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
  * Sadler  _et al._ [2006] L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore, and D. M. Stamper-Kurn, Nature 443, 312 (2006).
  * Hofferberth  _et al._ [2007] S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm, and J. Schmiedmayer, Nature 449, 324 (2007).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Trotzky  _et al._ [2012] S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Nature Physics 8, 325 (2012).
  * Cheneau  _et al._ [2012] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, Nature 481, 484 (2012).
  * Jepsen  _et al._ [2020] P. N. Jepsen, J. Amato-Grill, I. Dimitrova, W. W. Ho, E. Demler, and W. Ketterle, Spin transport in a tunable heisenberg model realized with ultracold atoms, Nature 588, 403 (2020).
  * White and Feiguin [2004] S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
  * Enss and Sirker [2012] T. Enss and J. Sirker, New Journal of Physics 14, 023008 (2012).
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Journal of Statistical Mechanics: Theory and Experiment 2005, P04010 (2005).
  * Deutsch [1991] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Sirker  _et al._ [2014] J. Sirker, N. P. Konstantinidis, F. Andraschko, and N. Sedlmayr, Phys. Rev. A 89, 042104 (2014).
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
  * Rigol  _et al._ [2006] M. Rigol, A. Muramatsu, and M. Olshanii, Phys. Rev. A 74, 053616 (2006).
  * Langen  _et al._ [2015] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Science 348, 207 (2015).
  * Barthel and Schollwöck [2008] T. Barthel and U. Schollwöck, Phys. Rev. Lett. 100, 100601 (2008).
  * Manmana  _et al._ [2009] S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, Phys. Rev. B 79, 155104 (2009).
  * Iucci and Cazalilla [2009] A. Iucci and M. A. Cazalilla, Phys. Rev. A 80, 063619 (2009).
  * Cazalilla [2006] M. A. Cazalilla, Phys. Rev. Lett. 97, 156403 (2006).
  * Barmettler  _et al._ [2009] P. Barmettler, M. Punk, V. Gritsev, E. Demler, and E. Altman, Phys. Rev. Lett. 102, 130603 (2009).
  * Barmettler  _et al._ [2010] P. Barmettler, M. Punk, V. Gritsev, E. Demler, and E. Altman, New Journal of Physics 12, 055017 (2010).
  * Cramer and Eisert [2010] M. Cramer and J. Eisert, New Journal of Physics 12, 055020 (2010).
  * Cassidy  _et al._ [2011] A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. Lett. 106, 140405 (2011).
  * Calabrese  _et al._ [2011] P. Calabrese, F. H. L. Essler, and M. Fagotti, Phys. Rev. Lett. 106, 227203 (2011).
  * Calabrese  _et al._ [2012a] P. Calabrese, F. H. L. Essler, and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2012, P07016 (2012a).
  * Calabrese  _et al._ [2012b] P. Calabrese, F. H. L. Essler, and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2012, P07022 (2012b).
  * Caux and Essler [2013] J.-S. Caux and F. H. L. Essler, Phys. Rev. Lett. 110, 257203 (2013).
  * Wouters  _et al._ [2014] B. Wouters, J. De Nardis, M. Brockmann, D. Fioretto, M. Rigol, and J.-S. Caux, Phys. Rev. Lett. 113, 117202 (2014).
  * Pereira  _et al._ [2014] R. G. Pereira, V. Pasquier, J. Sirker, and I. Affleck, Journal of Statistical Mechanics: Theory and Experiment 2014, P09037 (2014).
  * Bucciantini  _et al._ [2014] L. Bucciantini, M. Kormos, and P. Calabrese, Journal of Physics A: Mathematical and Theoretical 47, 175002 (2014).
  * Karrasch  _et al._ [2015] C. Karrasch, R. G. Pereira, and J. Sirker, New Journal of Physics 17, 103003 (2015).
  * Essler  _et al._ [2014] F. H. L. Essler, S. Kehrein, S. R. Manmana, and N. J. Robinson, Phys. Rev. B 89, 165104 (2014).
  * Chudzinski and Schuricht [2016] P. Chudzinski and D. Schuricht, Phys. Rev. B 94, 075129 (2016).
  * Calabrese and Cardy [2016] P. Calabrese and J. Cardy, Journal of Statistical Mechanics: Theory and Experiment 2016, 064003 (2016).
  * Sirker [2020] J. Sirker, SciPost Phys. Lect. Notes , 17 (2020).
  * Collura and Essler [2020] M. Collura and F. H. L. Essler, Phys. Rev. B 101, 041110 (2020).
  * Babenko  _et al._ [2021] C. Babenko, F. Göhmann, K. K. Kozlowski, J. Sirker, and J. Suzuki, Phys. Rev. Lett. 126, 210602 (2021).
  * Robinson  _et al._ [2021] N. J. Robinson, A. J. J. M. de Klerk, and J.-S. Caux, SciPost Phys. 11, 104 (2021).
  * Ruggiero  _et al._ [2021] P. Ruggiero, P. Calabrese, L. Foini, and T. Giamarchi, SciPost Phys. 11, 055 (2021).
  * Robertson  _et al._ [2022] N. F. Robertson, J. Surace, and L. Tagliacozzo, Phys. Rev. B 105, 195103 (2022).
  * Göhmann  _et al._ [2022] F. Göhmann, K. K. Kozlowski, J. Sirker, and J. Suzuki, SciPost Phys. 12, 158 (2022).
  * Pozsgay  _et al._ [2014a] B. Pozsgay, M. Mestyán, M. A. Werner, M. Kormos, G. Zaránd, and G. Takács, Phys. Rev. Lett. 113, 117203 (2014a).
  * Brockmann  _et al._ [2014] M. Brockmann, B. Wouters, D. Fioretto, J. D. Nardis, R. Vlijm, and J.-S. Caux, Journal of Statistical Mechanics: Theory and Experiment 2014, P12009 (2014).
  * Caux [2016] J.-S. Caux, Journal of Statistical Mechanics: Theory and Experiment 2016, 064006 (2016).
  * Alba and Calabrese [2016] V. Alba and P. Calabrese, Journal of Statistical Mechanics: Theory and Experiment 2016, 043105 (2016).
  * Bethe and Heitler [1934] H. Bethe and W. Heitler, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 146, 83 (1934).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
  * Barouch and McCoy [1971] E. Barouch and B. M. McCoy, Phys. Rev. A 3, 786 (1971).
  * Sengupta  _et al._ [2004] K. Sengupta, S. Powell, and S. Sachdev, Phys. Rev. A 69, 053616 (2004).
  * Pozsgay  _et al._ [2014b] B. Pozsgay, M. Mestyán, M. A. Werner, M. Kormos, G. Zaránd, and G. Takács, Phys. Rev. Lett. 113, 117203 (2014b).
  * Mestyán  _et al._ [2015] M. Mestyán, B. Pozsgay, G. Takács, and M. A. Werner,  2015, P04001 (2015).
  * Caianiello [1958] E. R. Caianiello, Il Nuovo Cimento 9, 569 (1958).

##  Appendix A Asymptotic formulas for large-$\Delta$ limit

As mentioned in the main text, to obtain the asymptotics of distance-$\ell$
correlators we must calculate the expectation value of an operator product
[see Eq.(20)] that can be written as a Pfaffian of pairwise contractions. The
problem is then reduced to calculating two-point correlations. In particular,
for values of $\ell$ up to 3, the contractions involved are $\langle
A_{0}B_{1}\rangle$, $\langle A_{2}B_{1}\rangle$, and $\langle
A_{0}B_{3}\rangle$. In this appendix, we provide some details about the
calculation of such contractions.

The integral form of $\langle A_{m}B_{n}\rangle$ is shown in Eq. (22). Let us
first focus on the stationary values given by

$\langle
A_{m}B_{n}\rangle_{\infty}=\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{-i[k(m-n)+2\theta_{k}]}\cos(2\theta_{k}+k).$
(28)

Using Euler’s formula, trigonometric relations and observing that
$\cos(2\theta_{k})=2(2+\Delta)\cos(k)/\varepsilon_{k}$, one can write Eq. (28)
as

$\displaystyle\langle A_{m}B_{n}\rangle_{\infty}$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\frac{2dk}{\pi\varepsilon_{k}^{2}}\left\\{\left\\{2\cos[k(r+1)]+\Delta\cos[k(r-1)]\right\\}\right.$
(29) $\displaystyle\times$
$\displaystyle\left.\left[2\cos(2k)+\Delta\right]\right\\},$

where $r=n-m$. Thus, the stationary values of the contractions $\langle
A_{0}B_{1}\rangle$, $\langle A_{2}B_{1}\rangle$, and $\langle
A_{0}B_{3}\rangle$ are

$\displaystyle\langle A_{0}B_{1}\rangle_{\infty}$ $\displaystyle=$
$\displaystyle\frac{\Delta^{2}-2}{\Delta^{2}},$ (30) $\displaystyle\langle
A_{2}B_{1}\rangle_{\infty}$ $\displaystyle=$ $\displaystyle\frac{1}{\Delta},$
(31) $\displaystyle\langle A_{0}B_{3}\rangle_{\infty}$ $\displaystyle=$
$\displaystyle\frac{4-\Delta^{2}}{\Delta^{3}}.$ (32)

Now, we turn to the time-dependent contributions given by

$\displaystyle\langle A_{m}B_{n}\rangle_{t}=i\int_{-\pi}^{\pi}\frac{dk}{2\pi}$
$\displaystyle\left\\{e^{-i[k(m-n)+2\theta_{k}]}\cos(2\varepsilon_{k}t)\right.$
(33) $\displaystyle\times$ $\displaystyle\left.\sin(2\theta_{k}+k)\right\\}.$

It is convenient to write Eq. (33) in the following form

$\displaystyle\langle A_{m}B_{n}\rangle_{t}$ $\displaystyle=$
$\displaystyle\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{4dk}{\pi\varepsilon_{k}^{2}}\left\\{e^{i2\varepsilon_{k}t}\\{2\sin[k(r+1)]\right.$
(34) $\displaystyle+$
$\displaystyle\left.\Delta\sin[k(r-1)]\\}\sin(2k)+\text{h.c.}\right\\}.$

To obtain the asymptotics of $\langle A_{m}B_{n}\rangle_{t}$, we carry out a
saddle-point analysis. The leading terms in the long-time regime are
determined by the stationary points of $\varepsilon_{k}$, which are $k_{0}=0$
and $\pm\pi/2$. Expanding the dispersion relation in the vicinity of $k_{0}=0$
and $\pm\pi/2$ up to fourth order, we respectively have

$\displaystyle\varepsilon_{k}$ $\displaystyle\approx$ $\displaystyle
4\left(\frac{\Delta}{2}+1\right)-\frac{8\Delta}{\Delta+2}k^{2}+\frac{8\Delta(\Delta^{2}-2\Delta+4)}{3(\Delta+2)^{3}}k^{4},$
$\displaystyle\varepsilon_{k}$ $\displaystyle\approx$ $\displaystyle
4\left(\frac{\Delta}{2}-1\right)+\frac{8\Delta}{\Delta-2}\left(k\pm\frac{\pi}{2}\right)^{2}$
(35) $\displaystyle-$
$\displaystyle\frac{8\Delta(\Delta^{2}+2\Delta+4)}{3(\Delta-2)^{3}}\left(k\pm\frac{\pi}{2}\right)^{4}.$

Assuming $t\gg 1$ and denoting $x^{2}=k^{2}t$, we can write $\langle
A_{m}B_{n}\rangle_{t}$ for $k_{0}=0$ as

$\displaystyle\langle A_{m}B_{n}\rangle_{t,0}$ $\displaystyle=$
$\displaystyle\frac{32e^{8i(1+\Delta/2)t}\xi(r,\Delta)}{\sqrt{\pi}t^{3/2}}\left(1+\frac{2}{\Delta}\right)^{-1/2}\int_{-\infty}^{\infty}dxe^{-16i\Delta
x^{2}/(\Delta+2)}\left\\{2x^{2}-\frac{1}{t}\left[\frac{2(r+1)^{3}+\Delta(r-1)^{3}}{6(r+1)+3\Delta(r-1)}\right.\right.$
(36) $\displaystyle+$
$\displaystyle\left.\left.\frac{4(4-8\Delta+\Delta^{2})}{3(2+\Delta)^{2}}\right]x^{4}+\frac{32i\Delta(4-2\Delta+\Delta^{2})}{3(2+\Delta)^{3}}x^{6}\right\\}+\text{h.c.}+O(t^{-7/2}),$

with $\xi(r,\Delta)$

$\xi(r,\Delta)=\frac{2(r+1)+\Delta(r-1)}{128\sqrt{\pi}(1+\Delta^{2}/4+\Delta)}\left(1+\frac{2}{\Delta}\right)^{1/2}.$
(37)

By a deformation of the complex contour, the integral in Eq. (36) can be
solved by transforming it into an exponentially decaying one. Doing so, we
have

$\displaystyle\langle A_{m}B_{n}\rangle_{t,0}$ $\displaystyle\approx$
$\displaystyle\xi(r,\Delta)\left\\{\left(1+\frac{2}{\Delta}\right)\frac{\sin\left[8\left(\frac{\Delta}{2}+1\right)t-\frac{\pi}{4}\right]}{t^{3/2}}\right.$
(38) $\displaystyle+$
$\displaystyle\left.\frac{\chi(r,\Delta)}{16\Delta^{2}}\frac{\cos\left[8\left(\frac{\Delta}{2}+1\right)t-\frac{\pi}{4}\right]}{t^{5/2}}\right\\},$

where the prefactor $\chi(\Delta)$ is given by

$\displaystyle\chi(r,\Delta)$ $\displaystyle=$
$\displaystyle\left(1+\frac{\Delta^{2}}{4}\right)\left[\frac{2(r+1)^{3}+\Delta(r-1)^{3}}{2(r+1)+\Delta(r-1)}-1\right]$
(39) $\displaystyle+$
$\displaystyle\Delta\left[\frac{2(r+1)^{3}+\Delta(r-1)^{3}}{2(r+1)+\Delta(r-1)}-\frac{11}{2}\right].$

The saddle points at $k_{0}=\pm\pi/2$ give contributions to the asymptotic
limit which can be obtained in an analogous manner to those for $k_{0}=0$. For
odd values of $|r|$, we find

$\displaystyle\langle
A_{m}B_{n}\rangle_{t,\pi/2}\approx\frac{\xi(r,-\Delta)}{2}\left\\{\left(\frac{2}{\Delta}-1\right)\frac{\sin\left[8\left(\frac{\Delta}{2}-1\right)t-\frac{3\pi}{4}\right]}{t^{3/2}}\right.$
(40) $\displaystyle+$
$\displaystyle\frac{\chi(r,-\Delta)}{16\Delta^{2}}\left.\frac{\cos\left[8\left(\frac{\Delta}{2}-1\right)t-\frac{3\pi}{4}\right]}{t^{5/2}}\right\\}(-1)^{(r+1)/2}.$

Finally, the complete expression for the contraction $\langle
A_{m}B_{n}\rangle$ is given by

$\langle A_{m}B_{n}\rangle=\langle A_{m}B_{n}\rangle_{\infty}+\langle
A_{m}B_{n}\rangle_{t,0}+2\langle A_{m}B_{n}\rangle_{t,\pi/2}.$ (41)