Equilibrium dynamics of the XX chain

TL;DR

This paper re-examines the equilibrium dynamics of the spin-1/2 XX chain using a new formalism, providing high-accuracy calculations of correlation functions and comparing them with existing theoretical and numerical results.

Contribution

It introduces a formalism based on quantum transfer matrix and thermal form factors to evaluate real-time space correlation functions in the XX chain.

Findings

Reproduces exact results in limiting cases

Matches asymptotic formulas from Riemann-Hilbert approach

Shows consistency with numerical Pfaffian evaluations and Luttinger liquid theory

Abstract

The equilibrium dynamics of the spin-1/2 XX chain is re-examined within a recently developed formalism based on the quantum transfer matrix and a thermal form factor expansion. The transversal correlation function is evaluated in real time and space. The high-accuracy calculation reproduces several exact results in limiting cases as well as the well-known asymptotic formulas obtained by the matrix Riemann-Hilbert approach. Furthermore, comparisons to numerical data based on a direct evaluation of the Pfaffian as well as to asymptotic formulas obtained within non-linear Luttinger liquid theory are presented.

Figures

Tables

Table 1. Table 1: The static correlation function ⟨ σ 1 − ​ ( 0 ) ​ σ 2 + ​ ( 0 ) ⟩ delimited-⟨⟩ subscript superscript 𝜎 1 0 subscript superscript 𝜎 2 0 \langle\sigma^{-}_{1}(0)\sigma^{+}_{2}(0)\rangle for h / J = 0.5 ℎ 𝐽 0.5 h/J=0.5 and various temperatures. The Fredholm determinant is evaluated by means of an approximation by a 256 × \times 256 matrix.

T	$\delta =0.01$	$\delta =0.05$	$\delta =0.1$	exact
0.1	-0.34556050	-0.31549848	-0.31547691	-0.31547744
0.5	-0.32068971	-0.30676642	-0.30675704	-0.30675704
1	-0.28351570	-0.27640356	-0.27639925	-0.27639925

Table 2. Table 2: Results for ⟨ σ 1 − ​ ( 0 ) ​ σ m + 1 + ​ ( 0 ) ⟩ delimited-⟨⟩ subscript superscript 𝜎 1 0 subscript superscript 𝜎 𝑚 1 0 \langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(0)\rangle at h = 0.1 ​ J ℎ 0.1 𝐽 h=0.1J ( < h c absent subscript ℎ 𝑐 <h_{c} ). Upper line for each T 𝑇 T shows the results from our new formula, while lower line shows the exact values. The Fredholm determinant is estimated by means of a 512 × 512 512 512 512\times 512 matrix.

T	m=0	m=1	m=2	m=3
0.1J	0.4920331867	-0.3178821170	0.2022244472	-0.1710719802
	0.4920331808	-0.3178821255	0.2022244549	-0.1710719748
0.5 J	0.4917867108	-0.3094178153	0.1915980505	-0.1508065219
	0.4917866789	-0.3094178186	0.1915980506	-0.1508065200
J	0.4910914106	-0.2793621277	0.1561766021	-0.1020132729
	0.4910913226	-0.2793621160	0.1561766025	-0.1020132730
5 J	0.4953703567	-0.0961933456	0.0185079888	-0.0036065564
	0.4953703745	-0.0961933456	0.0185079891	-0.0036065569

Table 3. Table 3: Results for ⟨ σ 1 − ​ ( 0 ) ​ σ m + 1 + ​ ( 0 ) ⟩ delimited-⟨⟩ subscript superscript 𝜎 1 0 subscript superscript 𝜎 𝑚 1 0 \langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(0)\rangle at h = 4.1 ​ J ℎ 4.1 𝐽 h=4.1J ( > h c absent subscript ℎ 𝑐 >h_{c} ). Upper line for each T 𝑇 T shows the results from our new formula, while lower line shows the exact value. The Fredholm determinant is estimated by means of a 512 × 512 512 512 512\times 512 matrix.

T	m=0	m=1	m=2	m=3
0.1J	0.0186123689	-0.0183512823	0.0176076636	-0.0164562085
	0.0186123688	-0.0183512823	0.0176076634	-0.0164562090
0.5 J	0.0765154951	-0.0704244629	0.0563287894	-0.0402407959
	0.0765154951	-0.0704244629	0.0563287892	-0.0402407964
J	0.1186618941	-0.0983455062	0.0624189803	-0.0338002579
	0.1186618940	-0.0983455062	0.0624189800	-0.0338002583
5 J	0.3181814080	-0.0830610673	0.0160102135	-0.0029751725
	0.3181814073	-0.0830610672	0.0160102133	-0.0029751730

Table 4. Table 4: Extrapolated (to T = 0 𝑇 0 T=0 ) values of ⟨ σ 1 − ​ ( 0 ) ​ σ m + 1 + ​ ( 0 ) ⟩ delimited-⟨⟩ subscript superscript 𝜎 1 0 subscript superscript 𝜎 𝑚 1 0 \langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(0)\rangle at h / J = 0.01 ℎ 𝐽 0.01 h/J=0.01 . Exact values at h = T = 0 ℎ 𝑇 0 h=T=0 are given for comparison.

	m=5	m=6	m=7	m=8	m=9	m=10	m=20
Extrapolated	-0.132188	0.119687	-0.111465	0.103810	-0.0982175	0.0929276	0.0659939
Exact	-0.132195	0.119691	-0.111467	0.103807	-0.0982084	0.0929116	0.0657593

Equations

$${\cal H}=J\sum_{i=1}^{L}\Bigl{(}\sigma_{i}^{x}\sigma_{i+1}^{x}+\sigma_{i}^{y}\sigma_{i+1}^{y}\Bigr{)}-\frac{h}{2}\sum_{i=1}^{L}\sigma^{z}_{i},$$

$$h_c=4J.$$

$$\langle\sigma_{1}^{-}(0)\sigma_{m+1}^{+}(t)\rangle_{J,h,T}\\ =\lim_{L\rightarrow\infty}\frac{{\rm tr}_{1,\cdots,L}\bigl{(}{{\rm e}}^{-(\frac{1}{T}+it){\cal H}}\sigma^{-}_{1}{{\rm e}}^{it{\cal H}}\sigma^{+}_{m+1}\bigr{)}}{{\rm tr}_{1,\cdots,L}\bigl{(}{{\rm e}}^{-\frac{1}{T}{\cal H}}\bigr{)}}.$$

$$t_c=\frac{m}{4J}.$$

$$ϵ(p)=h-4J\cos (p).$$

$$⟨{\sigma }_1^{-}(0){\sigma }_{m+1}^{+}(t)⟩=(-1{)}^{m+1}\mathcal{A}(m)$$

$$\times \sum _{n=1}^{\infty }\frac{(-1{)}^{n-1}}{n!(n-1)!}{\int }_{\mathcal{E}}\prod _{r=1}^n\mathrm{d}{p}_r\mu (p_r){\mathrm{e}}^{i(m{p}_r-tϵ(p_r))}$$

$$\displaystyle\times\int_{\bar{\cal E}}\prod_{s=1}^{n-1}{\rm d}q_{s}\bar{\mu}(q_{s}){\rm e}^{-i(mq_{s}-t\epsilon(q_{s}))}{\cal D}_{p}^{2}\bigl{(}\{p_{r}\}_{r=1}^{n},\{q_{s}\}_{s=1}^{n-1}\bigr{)},$$

$$\mu (p)=\frac{{\mathrm{e}}^{{\sigma }_{+}(p)}}{2\pi (1-{\mathrm{e}}^{ϵ(p)/T})},\bar{\mu }(q)=\frac{{\mathrm{e}}^{-{\sigma }_{-}(q)}}{2\pi (1-{\mathrm{e}}^{-ϵ(q)/T})},$$

$$\displaystyle\sigma(p)=\int_{\cal E}\frac{dq}{2\pi i}\frac{1}{\tan\bigl{(}\frac{p-q}{2}\bigr{)}}\ln\Bigl{(}\frac{1+{\rm e}^{-\epsilon(q)/T}}{1-{\rm e}^{-\epsilon(q)/T}}\Bigr{)},$$

$$\displaystyle{\cal D}_{p}\bigl{(}\{p_{r}\}_{r=1}^{n},\{q_{s}\}_{s=1}^{n-1}\bigr{)}=$$

$$\displaystyle\phantom{ab}\frac{\underset{1\leq j<k\leq n}{\prod}\sin\bigl{(}\frac{p_{j}-p_{k}}{2}\bigr{)}\underset{1\leq j<k\leq n-1}{\prod}\sin\bigl{(}\frac{q_{j}-q_{k}}{2}\bigr{)}}{\prod_{j=1}^{n}\prod_{k=1}^{n-1}\sin\bigl{(}\frac{p_{j}-q_{k}}{2}\bigr{)}},$$

$$\tilde{V}(q_i,q_j)$$

$$V(q_i,q_j)$$

$$\tilde{v}(q_i)$$

$$v(q_i)$$

$$\Omega (m,t)$$

$$P(q_i,q_j)$$

$$\varphi(p,q)=\frac{{\rm e}^{i\frac{q-p}{2}}}{\sin\bigl{(}\frac{q-p}{2}\bigr{)}},\quad u_{m,t}(p)=\frac{i}{2}(\frac{m}{t}p-\epsilon(p)).$$

$$(\hat{V}f)(q)$$

$$(\hat{P}f)(q)$$

$$\begin{gathered} ⟨{\sigma }_1^{-}(0){\sigma }_{m+1}^{+}(t)⟩=(-1{)}^{m+1}\mathcal{A}(m)\Omega (m,t) \\ \times {\det }_{\bar{\mathcal{E}}}(1+\hat{P}-\hat{V}). \end{gathered}$$

$$\widetilde{V}(q_{i},q_{j})=\frac{{\rm e}^{i\frac{(q_{i}-q_{j})}{2}}\widetilde{v}(q_{j})-{\rm e}^{-i\frac{(q_{i}-q_{j})}{2}}\widetilde{v}(q_{i})}{\sin\bigl{(}\frac{q_{i}-q_{j}}{2}\bigr{)}}.$$

$$p_{\pm}=\frac{\pi}{2}\pm i\operatorname{arch}\bigl{(}\frac{m}{4Jt}\bigr{)}.$$

$$p_{\pm}=\frac{\pi}{2}\pm\arccos\bigl{(}\frac{m}{4Jt}\bigr{)}.$$

$$q_j=\arccos (\frac{h+ij\pi T}{4J})(j\in \mathbb{Z}),$$

$$1+{\mathrm{e}}^{-ϵ(q_{2j-1})/T}=0,1-{\mathrm{e}}^{-ϵ(q_{2j})/T}=0.$$

$$q_j^u=\{\begin{array}{ll}{q}_j& j\le 0,\\ -q_j& j\ge 1,\end{array}{q}_j^d=\{\begin{array}{ll}-q_j& j\le 0,\\ q_j& j\ge 1,\end{array}$$

$$⟨{\sigma }_1^{-}(0){\sigma }_1^{+}(0)⟩=\{\begin{array}{ll}0.0186123688\cdots & \text{exact},\\ 0.0186094\cdots & n=64,\\ 0.0186123689\cdots & n=512.\end{array}$$

$$\tilde{v}(q)={\int }_{{\mathcal{C}}_p}dp\mu (p){\mathrm{e}}^{2t{u}_{m,t}(p)}\phi (p,q)+4\pi i{\mathrm{e}}^{2t{u}_{m,t}(q)}\mu (q),$$

$$-\frac{16{\pi }^2}{\Omega (m,t)}{\mathrm{e}}^{t{u}_{m,t}(q_i)+t{u}_{m,t}(q_j)}\mu (q_i)\mu (q_j).$$

$$\Omega (m,t)$$

$$+\sqrt{\frac{\pi }{2\left|u_{m,t}^{\prime \prime }(p_{-})t\right|}}\mu (p_{-}){\mathrm{e}}^{i(m{p}_{-}-ϵ(p_{-})t)}.$$

$$\left|\mu (p_{-})\right|=\left|\frac{{\mathrm{e}}^{{\sigma }_{+}(p_{-})}}{2\pi (1-{\mathrm{e}}^{ϵ(p_{-})/T})}\right|\gg \left|\mu (p_{+})\right|$$

$$\begin{gathered} ⟨{\sigma }_0^{-}(0){\sigma }_m^{+}(t){⟩}^{(0)}\propto {\left(\frac{\pi T}{v}\right)}^{1/2} \\ \times \frac{{\text{e}}^{i\pi x/a}}{{\mathrm{sh}}^{\frac{1}{4}}\left(\frac{\pi T}{v}(x-vt-i{\eta }^{+})\right){\mathrm{sh}}^{\frac{1}{4}}\left(\frac{\pi T}{v}(x+vt+i{\eta }^{+})\right)}. \end{gathered}$$

$$\begin{gathered} ⟨{\sigma }_0^{-}(0){\sigma }_m^{+}(t){⟩}^{(1)} \\ \propto {\text{e}}^{i{\bar{p}}_{+}x}\cos (k_{F}x)⟨{\sigma }_0^{-}(0){\sigma }_m^{+}(t){⟩}^{(0)}⟨d(0,0)d^{†}(x,t)⟩. \end{gathered}$$

Full text

I Introduction

An increasing knowledge has been acquired on the correlation functions of the simplest models which are integrable in the sense of the Yang-Baxter relation. The vertex operator approach (VOA) has definitely triggered the recent advance Jimbo and Miwa (1995). It enables us to evaluate ground state correlation functions, for instance, of the XXZ chain in the antiferromagnetic regime in a vanishing magnetic field. A multiple-integral formula for the reduced density matrix has been derived naturally within this framework Jimbo et al. (1992). The evaluation of the two- and four-spinon contributions to the dynamical structure factor of the massive XXZ chain Bougourzi et al. (1998); Abada et al. (1997) and the four spinon ones in the XXX limit Caux and Hagemans (2006) were important outcomes of this method.

A complementary approach to the analysis of ground-state correlation functions of the XXZ chain starts with its algebraic Bethe ansatz solution. It utilizes the solution of the quantum inverse problem Kitanine et al. (2000a) and a special determinant formula for the scalar product of on-shell and off-shell Bethe vectors Slavnov (1989). The algebraic Bethe ansatz approach successfully confirmed Kitanine et al. (2000b) the multiple-integral formula for the density matrix obtained within the VOA. It properly takes into account the effect of a finite magnetic field in the symmetry direction of the chain. This approach can be also used to derive and analyze form factor series for dynamical correlation functions in the thermodynamic limit Kitanine et al. (2009); Kitanine et al. (2011a, b); Kitanine et al. (2012); Dugave et al. (2015); Kozlowski and Maillet (2015); Kozlowski (2018). These developments culminated in the extraction, on the basis of exact and first principle calculations, of the long-distance and large-time asymptotic behavior of two point functions in the XXZ chain at zero temperature Kozlowski (2019a). Also, it was possible to grasp on exact grounds, the full structure of the edge singularities of the spin structure factors and spectral functions in this model Kozlowski (2019b).

In a parallel development, field theoretical approaches have been extended to include non-linearities in the spectrum Imambekov and Glazman (2009); Imambekov et al. (2012). It has been argued that using non-linear Luttinger liquid theory is crucial to obtain proper universal results for threshold singularities in spectral functions and for the long time asymptotics of dynamical correlation functions. For integrable lattice models, in particular, a combination of non-linear Luttinger liquid theory and Bethe ansatz has allowed to derive parameter-free results for edge singularities and high-energy tails of spin structure factors and spectral functions Pereira et al. (2006, 2007, 2008, 2009). Of particular relevance in the present context are results for the asymptotics of the dynamical longitudinal spin-spin correlation function Karrasch et al. (2015) and the transverse spin-spin structure factor Karimi and Affleck (2011) for the XXZ chain at low but finite temperatures. These results play the same role for dynamical correlations as the well-known conformal field theory formulas do in the static case. In contrast to the latter, however, a verification of the asymptotic formulas based on a fully microscopic calculation is so far lacking in the dynamical case at low but finite temperature.

It has been observed that the simplest multiple integrals representing the reduced density matrix of short sub-segments of the XXZ chain factorize Boos and Korepin (2001); Boos et al. (2003); Sakai et al. (2003); Boos et al. (2004, 2005). The algebraic structure of the ground state expectation values of all finite-range operators was finally understood with the discovery of a ‘hidden Grassmann structure’ on the space of such operators in Boos et al. (2007a, b, 2009). While the structure turned out to be relevant for the analysis of the short-range correlation functions Sato et al. (2011); Miwa and Smirnov (2019), the form factor series seem to be more efficient in analyzing the long time, large distance behavior of the two-point functions. The microscopic verification Kitanine et al. (2011a); Kozlowski and Maillet (2015) of the predictions of conformal field theory Cardy (1986) for the XXZ chain is an important application.

Thanks to the similarity in the structure of the row-to-row transfer matrix of the six-vertex model and the quantum transfer matrix (QTM) of the XXZ chain Klümper (1993), many results for ground state correlation functions of the XXZ chain could be generalized to finite temperatures. This concerns in the first place the multiple-integral representation for the reduced density matrix of sub-chains Göhmann et al. (2004); Göhmann et al. (2005a, b). The multiple integrals factorize even at finite temperatures Boos et al. (2006, 2007c, 2008), which finds its explanation again in the hidden Grassmann structure Jimbo et al. (2009). Moreover, for the static two-point functions at finite temperature so-called thermal form-factor expansions, involving form factors of the QTM, were introduced and studied in Dugave et al. (2013, 2014, 2016). As an interesting outcome we mention the evaluation of the static two-point correlation functions in the low-temperature limit in terms of higher particle-hole excitations Dugave et al. (2016). They were conjectured to correspond to the two-, four- and six-spinon contributions in the VOA.

As compared to the ground-state case, the inclusion of the time dependence into the thermal form factors series for the two-point functions requires slightly more thought. It was realized in Sakai (2007) that the solution to the inverse problem for the QTM has to be adopted for this purpose. By combining the lattice realization of the two-point functions suggested in Sakai (2007) with the thermal form factor expansion of Dugave et al. (2013), we recently proposed Göhmann et al. (2017) a new scheme for the explicit evaluation of dynamical correlation functions of quantum spin systems in thermal equilibrium. As a first test for the validity of the new scheme, we re-derived the well-known formula Niemeijer (1967) for the longitudinal equilibrium correlation function of the spin-$\frac{1}{2}$ XX chain. We also obtained a novel thermal form factor series for the transversal two-point correlation function $\langle\sigma_{1}^{-}(0)\sigma_{m+1}^{+}(t)\rangle$ which will be analyzed in this work.

It may seem rather amazing that for a simple model as the XX chain, whose Hamiltonian can be expressed in terms of non-interacting spinless Fermions Lieb et al. (1961), the calculation of its transverse dynamical two-point function at finite temperature still poses interesting questions. It is considerably harder than the longitudinal case, because the two-point function is non-local in terms of the Jordan-Wigner Fermions. Consequently, we are facing the problem of evaluating large determinants. Analytic studies are therefore mainly restricted to the case where each matrix element has a simple structure so that the Szegö theorem McCoy et al. (1971); Brandt and Jacoby (1976) can be applied. Results at finite temperatures, that were obtained within approaches based on the use of Fermions, comprise asymptotic formulae for high temperatures Brandt and Jacoby (1976); Perk and Capel (1977) as well as the numerical evaluation for finite open systems. The latter seems most efficient if the correlation function for a finite length chain is represented in terms of a Pfaffian of two-point functions of auxiliary Fermions Derzhko and Krokhmalskii (1997); Derzhko et al. (2000). Fermion algebra and generalized Wick theorem lead to an important consequence in the XY chain and its limit, a set of nonlinear difference-differential equations. See Perk (1980); McCoy et al. (1983); Perk et al. (1984); Perk and Au-Yang (2009) and references therein.

A completely different route to the evaluation of the transverse dynamical two-point functions was taken in Colomo et al. (1993). Following the strategy of Korepin and Slavnov (1990) the authors of Colomo et al. (1993) used the coordinate Bethe ansatz to derive a Fredholm determinant representation of the correlation function in the thermodynamic limit. The integral operator in the Fredholm determinant is of integrable type Korepin et al. (1993). This formulation was then directly suitable for an asymptotic analysis at long times $t$ and large distances $m$ by means of an associated matrix Riemann-Hilbert problem amenable to a ‘non-linear steepest descent method’ Deift and Zhou (1992). In Its et al. (1993) it was shown that the long time, large distance behavior of the two-point function for a fixed ratio $m/t$ in a weak magnetic field ($h<4J$) goes like $\langle\sigma_{1}^{-}(0)\sigma_{m+1}^{+}(t)\rangle\sim At^{\nu}{\rm e}^{-m/\xi}$. The explicit expressions for $\xi$ and $\nu$ were obtained in the space-like ($m>4Jt$) and time-like ($m<4Jt$) asymptotic regimes. The result was subsequently extended to strong magnetic field ($h>4J$) Jie (1998).

This is one of three companion papers in which we revisit the problem of the evaluation of the transversal two-point function at finite temperatures. Our starting point will be the novel thermal form factor expansion for $\langle\sigma_{1}^{-}(0)\sigma_{m+1}^{+}(t)\rangle$ that was derived in a previous communication Göhmann et al. (2017). In Göhmann et al. (2019a, b) we study two different aspects of asymptotic analysis: On the one hand, the high-temperature analysis for all times and all distances and, on the other one, the evaluation of the amplitude $A$ in the space-like regime at weak magnetic field ($h<4J$). In this communication, we shall concentrate on the numerical evaluation of the series. We will remind the reader in Section II that the form factor series can be neatly re-summed in terms of a Fredholm determinant different from the one in Colomo et al. (1993). The kernel of the integral operator that defines the Fredholm determinant is a strongly oscillating function for large values of $m$ and $t$. We shall demonstrate that, nevertheless, the idea presented in Bornemann (2010) of a direct evaluation of the determinant can be applied, if combined with an appropriate choice of the integration contour in the complex plane. Then the real time evaluation can be performed in a stable manner until rather large times $t\gg\frac{m}{4J}$. Consequently, we are able to numerically check the asymptotic results of Its et al. (1993); Jie (1998). We are also able to suggest a higher-order correction for $h<4J$ from our numerical analysis. The advantage of the method of numerical evaluation proposed in this work is that it works directly in the thermodynamic limit and is free of finite size effects.

The paper is organized as follows. In Section II, we start from a brief review of the results in Göhmann et al. (2017) and then derive the novel Fredholm determinant representation of the transverse correlation function. We will discuss the analytic properties of the functions occurring in the determinant in Section III. The steepest-descent paths in the space-like and time-like regimes will be explained. They are crucial for the numerics, especially in the large-time dynamics and for the long-distance correlations. Based on these preparations, we present the result of our numerical analysis in the massive phase in Section IV. The kinematic poles will become important at a later stage, and we shall argue how they can be treated numerically. The situation becomes more complicated in the massless phase due to the presence of Fermi points on the real axis. This will be discussed in Section V. In Section VI numerical data are presented showing the difference in the oscillation amplitudes of correlation functions at odd and even distances. These differences are then explained based on asymptotic results obtained within non-linear Luttinger liquid theory. A comparison with the analytic predictions based on the Fredholm determinant representation derived in Colomo et al. (1993) is given in Section VII. Thanks to the high-precision calculation, we identify a higher order correction in the massless phase. In Section VIII we compare the new numerical scheme proposed here and a previously existing method based on the Pfaffian representation. We summarize our results and point out perspectives in Section IX. Some of the technical details as well as a comparison with the exact static values are supplemented in several appendices.

II A representation of the transverse correlation function by a Fredholm determinant

We consider the spin-$\frac{1}{2}$ XX chain,

[TABLE]

with periodic boundaries in the thermodynamic limit $L\rightarrow\infty$. We assume the strength of the spin-spin interaction $J$ and the Zeeman magnetic field $h$ to be positive. The model has two phases in its ground state phase diagram, a critical phase (low $h$) with finite magnetization that depends of the value of $h$, and a massive phase (high $h$) in which the magnetization is saturated to its maximal possible value $1/2$. The two phases are separated by the critical field

[TABLE]

The quantity of interest in this work is the transverse correlation function

[TABLE]

We shall consider it as a function of distance $m$ and time $t$ that depends parametrically on $J$ and $h$ and on the temperature $T$. The subscripts $J,h,T$ will often be suppressed below.

We divide the $m$-$t$ space-time plane into two regimes: the space-like regime with $t<t_{c}$ and the time-like regime with $t>t_{c}$, where

[TABLE]

We shall see that the qualitative behavior of the transverse correlation function falls into 4 categories: massive time-like, massive space-like, massless time-like and massless space-like. Our aim is to understand the different categories quantitatively based on the QTM method.

The application of the QTM method and the derivation of a thermal form factor expansion for the transverse correlation function (3) of the XX model was described in detail in Section 3.5 of Göhmann et al. (2017). We give a brief summary of the result obtained there in Appendix A.

II.1 A representation by a Fredholm determinant

In our previous work Göhmann et al. (2017) we expressed the integrals in the thermal form factor series in terms of rapidity variables $\lambda$ (see (29)). This is natural for integrable systems in general, but for the numerical study presented here, it turns out to be more convenient to switch to momentum variables $p$. Then the one-particle dispersion takes its familiar functional form

[TABLE]

By the common abuse of notation we use the same symbol $\epsilon$ to denote the energy as a function of the momentum or rapidity variables. The substitution $\lambda\mapsto p(\lambda)$ in (29) with $p(\lambda)$ according to (26) then leads to

[TABLE]

where ${\cal A}(m)$ is defined in (28),

[TABLE]

and $\sigma_{+}(p)$ (respectively $\sigma_{-}(p)$) is the boundary value as $p$ approaches ${\cal E}$ from above (respectively  below). We call $\{p_{r}\}_{r=1}^{n}$ the hole momenta and $\{q_{s}\}_{s=1}^{n-1}$ the particle momenta. The contour ${\cal E}$ is a straight line $[-\pi+i\delta,\pi+i\delta]$ in the massive phase, $h>h_{c}$. Similarly, $\bar{\cal E}$ is a contour just below the real axis. When $h<h_{c}$, the contours will be slightly deformed in order to avoid the Fermi points111The Fermi momentum $k_{F}$ with the proper dimension will be introduced later in Section VI. $\pm p_{F}$ (defined by $\epsilon(\pm p_{F})=0$) as shown in Figure 1. The contours ${\cal E}$ and $\bar{\cal E}$ are the images of ${\cal C}$ and $\bar{\cal C}$ (see Appendix A) under the transformation $\lambda\mapsto p(\lambda)$.

Using a technique developed in Korepin and Slavnov (1990) we can rewrite (5) as an explicit factor times a Fredholm determinant of a special form. For this purpose, let us prepare some notations in advance,

[TABLE]

where

[TABLE]

We further define integral operators,

[TABLE]

The operator $\hat{P}$ is obviously a one-dimensional projector. After these preparations, by a calculation very similar to that in Göhmann et al. (2019a), we obtain the following formula.

Theorem**.**

The transverse correlation function of the spin-$\frac{1}{2}$ XX model can be represented using a Fredholm determinant that involves the integral operator $\hat{P}-\hat{V}$,

[TABLE]

The symbol ${\rm det}_{\cal S}$ means the Fredholm determinant where the integral operator is acting on functions supported on ${\cal S}$.

Equation  (10) is a novel analytic expression representing the transverse dynamical correlation function (3) of the infinite system at finite temperatures. In recent years it has become increasingly clear Bornemann (2010) that Fredholm determinants can be very efficiently evaluated numerically by simply approximating them by determinants of finite matrices, very much like integrals are approximated by finite sums. The numerical error decreases with the number of discretization points. It has been estimated for kernels with various analytic properties in Bornemann (2010).

When we evaluate (10) numerically, we are considering the kernel as a function of the momentum variables that depends parametrically on distance $m$ and time $t$. When $m$ and $t$ become large, the kernel function becomes a rapidly oscillating function of the momentum variables on the contour $\bar{\cal E}$. This affects the numerical accuracy of the calculation in a similar way as in the case of integrals over rapidly oscillating functions. As we shall see, the cure to this problem as well is similar as in the case of integrals over rapidly oscillating functions. We have to deform the integration contour in the complex plane into a saddle point contour. Since we are dealing with meromorphic functions with infinitely many poles the deformation is an intricate issue. This will be, on a technical level, the main subject of the discussions in the following sections.

We would like to remark that $\hat{V}$ belongs to the class of integrable integral operators Korepin et al. (1993). This is easily seen if we rewrite its kernel in the following form,

[TABLE]

For integrable integral operators there exist powerful tools in order to analyze the asymptotic behavior analytically Deift and Zhou (1992). A representation of the correlation function (3) involving a Fredholm determinant with a similar but manifestly different kernel, which belongs to the family of integrable integral operators as well, was analyzed in Its et al. (1993). The leading terms of the large space and time asymptotic behavior of the transverse correlation function of the XX model have been successfully derived from that representation. Here, however, we are interested in a quantitative study on an arbitrary space scale and time scale.

III The analytic properties of integrands and the steepest

descent paths

The contours ${\cal E},\bar{\cal E}$ are optimal choices in the static limit. As time evolves, the phase factors ${\rm e}^{\pm i(mp-\epsilon(p)t)}$ bring instability and we eventually need to deform the contours for a reliable calculation. One may encounter singularities of the integrands during the deformation. We thus have to find a balance between the advantage of reducing the instability from ${\rm e}^{\pm i(mp-\epsilon(p)t)}$ and the cost of passing through many singularities of the integrands. Below we shall discuss the analytic properties of the integrands and the optimal choice of integration paths.

III.1 The saddle points and the steepest descent paths

We need to take account of the steepest descent paths for ${\rm e}^{\pm u_{m,t}(p)t}$ in the asymptotic region $m,t\gg 1$. The locations of saddle points for the system in the space-like regime are qualitatively different from those for the system in the time-like regime. On the other hand, they do not depend on whether the system is in the massive or in the massless phase.

III.1.1 The space-like regime $t<t_{c}$

The saddle points $p_{\pm}$ in the space-like regime lie at

[TABLE]

The red points in Figure 2 denote them. The dashed curve corresponds to the loci of points $p$ such that \operatorname{Im}\bigl{(}u_{m,t}(p)\big{)}=\operatorname{Im}\big{(}u_{m,t}(p_{+})\big{)}, and the shaded region satisfies $|{\rm e}^{u_{m,t}(p)}|>|{\rm e}^{u_{m,t}(p_{+})}|$. Thus, the steepest descent path for the hole momentum passes through $p_{+}$ horizontally and never enters into the shaded region. We have an upside-down figure for the path of the particle momentum which passes through $p_{-}$, as it has the conjugate phase ${\rm e}^{-u_{m,t}(p)t}$. As a result, the optimal path for the hole (particle) momentum in the space-like regime is similar to ${\cal E}$ ($\bar{\cal E}$): it stays above (below) the real axis.

III.1.2 The time-like regime $t>t_{c}$

In this regime the saddle points $p_{\pm}$ are on the real axis,

[TABLE]

They are depicted by red dots in the right panel of Figure 2. The shaded region indicates $|{\rm e}^{u_{m,t}(p)}|>1$. Thus, for the hole momentum, the steepest descent path runs through $p_{+}$ and $p_{-}$, while staying inside of the un-shaded region. Accordingly, it cannot stay above the real axis, but passes below the real axis for $p_{-}<\operatorname{Re}(p)<p_{+}$. The figure for the particle momentum is obtained by turning the figure for the hole momentum upside down.

III.2 The analyticity of $\mu$ and $\bar{\mu}$

We only have to know $\mu(p)$ ($\bar{\mu}(p)$) in (II.1) as functions of $p$ in the upper (lower) half plane in the space-like regime. As was already discussed, one needs to continue them analytically into the whole complex plane in the time-like regime. The analytic properties turn out to be quite different for the massless ($h<h_{c}$) phase and for the massive ($h>h_{c}$) phase.

III.2.1 The massive regime $h>h_{c}$

We set

[TABLE]

so that

[TABLE]

We further define “upper roots” $q^{u}_{j}$ and “lower” roots $q^{d}_{j}$ by

[TABLE]

so that $\operatorname{Im}(q_{j}^{u})>0$ and $\operatorname{Im}(q_{j}^{d})<0$.

The roots closest to the real axis are q^{u}_{0}=i\operatorname{arch}\bigl{(}\frac{h}{4J}\bigr{)} and $q^{d}_{0}=-q^{u}_{0}$.

The following analytic properties of $\mu$ are derived in B.

Lemma 1**.**

*The function $\mu(p)$ has

i) simple poles at $q^{u}_{j}$ in the upper half plane,

ii) double poles at $q^{d}_{2j+1}$ in the lower half plane,

iii) simple zeros at $q^{d}_{2j}$ in the lower half plane.

There is a strip including the real axis which is free of zeros and poles (of width 2\operatorname{arch}\bigl{(}\frac{h}{4J}\bigr{)}).*

Figure 3 (left panel) illustrates the situation. The red crosses represent single poles, while blue triangles are double poles. The circles denote the single zeros.

The zeros and poles for $\bar{\mu}(p)$ are obtained by taking the mirror image with respect to the real axis.

Note that the width of the analytic strip is independent of the temperature in the massive phase. Once a finite magnetic field is imposed, the width stays finite. Thanks to this fact, the evaluation of the correlation function is simpler in the massive regime.

III.2.2 The massless regime $h<h_{c}$

We denote by $\{q^{r}_{j}\}_{j\in\mathbb{Z}}$, the (right) “Bethe roots”, which are identical to the $q_{j}$ in eq. (14). Note that $p_{F}(=q^{r}_{0})\leq\Re{\rm e}(q^{r}_{j})<\frac{\pi}{2}$. The (left) “Bethe roots” are defined by $q^{\ell}_{j}=-q^{r}_{-j}$.

In B   we show that $\mu$ ($\bar{\mu}$) has the following analytic properties.

Lemma 2**.**

*The function $\mu(p)$ has

i) simple poles at $q^{r}_{-2j}$ and at $q^{\ell}_{-2j}$ ($j\geq 1$) in the upper half plane,

ii) a simple pole at $p=p_{F}$ and a simple zero at $p=-p_{F}$ on the real axis,

iii) double poles at $q^{r}_{2j+1}$ and at $q^{\ell}_{2j+1}$ ($j\geq 0$) in the lower half plane,

iv) simple zeros at $q^{r}_{2j}$ and at $q^{\ell}_{2j}$ ($j\geq 1$).

The zeros and poles of $\bar{\mu}(p)$ are obtained by taking the mirror image w.r.t. the origin.*

The massless case is illustrated in the right panel of Figure 3. The differences $q^{\ell}_{-2j}-q^{\ell}_{-2j-2}$ and $q^{\ell}_{2j+1}-q^{\ell}_{2j-1}$ become ${\cal O}(T)$ as $T\rightarrow 0$. Thus, singularities will approach towards the real axis as $T\rightarrow 0$ and the choice of the contours becomes difficult. This poses a technical problem in the massless regime.

IV Numerical study of the massive regime $h>h_{c}$

In this and in the next sections we present the results of our study of the transverse correlation function based on the numerical evaluation of (10). We follow the proposal in Bornemann (2010) and evaluate the Fredholm determinant by a Gauss-Legendre integration with $n$ discrete points.

We shall start from the simpler case $h>h_{c}$. In this case the Fermi points are away from the real axis. Thus, we can separately treat two technical problems, how to deal with the Fermi points and how to deal with the steepest descent path.

IV.1 The space-like regime $t<t_{c}$

In spite of having discussed the steepest descent path above, we can use simple straight integration contours in order to produce accurate results, as long as $m\lesssim 10$ in this regime. In order to demonstrate this, we first consider the static limit where many results are available. The exact static short-range correlation functions at arbitrary $T$ and $h$ were obtained, e.g., in Boos et al. (2008) (for the XXZ model in general). We compare them with the static results obtained from (10) in C. The results match with reasonable precision. We have to choose, of course, the appropriate number of discretized points $n$ to achieve agreement. For example, fixing $h=4.1J,T=0.1J$ and $m=0$ we find

[TABLE]

Thus, we can say that with $n=512$, eq. (10) practically reproduces the exact value. For larger segments, some indirect evidence is presented in D. Encouraged by this success, we adopt the straight contours in the space-like regime as they are conveniently simple for the calculation. Indeed, they produce accurate results even beyond $t_{c}$.

Figures 4 and 5 show examples for $h=4.1J$, $m=1$ and $m=8$ for various $T$. The horizontal axes indicate $Jt$. Although the space-like regime is limited to $Jt<0.25$ when $m=1$ and $Jt<2$ when $m=8$, the straight contours produce stable values even after $t_{c}$.

The real part stays almost flat in the space-like regime. This can be better seen in the case of $m=8$ in Figure 5. The amplitude is enhanced around $t=t_{c}$.

For larger values of $m$ (typically $m\sim 100$), we need to adopt contours which take into account the steepest descent paths. The saddle points in eq. (12) are located away from the real axis, and we shift the contour ${\cal E}$ ( $\bar{{\cal E}}$ ) to a straight line passing though $p_{+}$ ($p_{-}$). The integration contours cross poles of $\mu(p)$ and $\bar{\mu}(q)$ if $\frac{h}{4}<\frac{m}{4t}$, which is clearly seen in Figure 3 (left). This modifies our formula (10) slightly, following the argument in Dugave et al. (2016), which will be summarized in E for the reader’s convenience.

IV.2 The time-like regime $t>t_{c}$

The calculation with the straight integration contours gradually becomes unreliable for larger $Jt$. An example is given in Figure 6.

We thus deform the contours in the time-like regime. The most naive choices for the hole momentum ($p$) and for the particle momentum $(q)$ may be the ones depicted in Figure 7, which respect the steepest descent paths discussed in Section III.1.

The paths are determined only by the phases ${\rm e}^{tu_{m,t}(p)}$ and ${\rm e}^{-tu_{m,t}(q)}$. One, however, needs to be careful: Since $q$ and $p$ are swapped, there appears an additional term in $\tilde{v}$,

[TABLE]

due to a pole of $\varphi(p,q)$ at $p=q$. Here ${\cal C}_{p}$ denotes the red path in Figure 7. As we interpret $p$ ( respectively $q$) as the momentum of a hole ( respectively a particle), we refer to the singularity as the kinematic pole in analogy with scattering theory. Consequently, the kernel $K(q_{i},q_{j})$ of the integral operator $\hat{K}=\hat{P}-\hat{V}$ contains a term

[TABLE]

This becomes very large for $Jt\gg 1$ if $\operatorname{Im}(q)\sim{\cal O}(1)$ and makes the calculation again unstable. We thus choose the red contour in Figure 7 for the $p$ variables while we adopt a straight contour for the $q$ variable such that $\operatorname{Im}(q)\sim 0^{-}$ and $\operatorname{Im}(p)<\operatorname{Im}(q)$ for $p_{-}<\operatorname{Re}(q)<p_{+}$.

Still, there is a problem in dealing with the intersection points $p_{+},p_{-}$. We simply exclude them from the set of discretized sampling points in the integrals over the $q$ variable.

After the above modification and by increasing the number of sampling points for the Fredholm determinant ($n=100-512$, typically), a large scale calculation is possible. As an illustration, the plots for $m=2,h=4.1J$ are shown in Figure 8.

One immediately notices the very long period oscillation. This may be attributed to the oscillatory behavior of $\Omega(m,t)$ of which the period seems to depend only weakly on temperature (Figure 9).

The amplitude of $\Omega(m,t)$ decreases slowly in time and the decay of the amplitude of the correlation functions mainly comes from that of the Fredholm determinant, c.f. Figure 10, especially for higher $T$.

Unlike on temperature, the period of oscillation obviously does depend on the magnetic field, and it grows as $h$ approaches $h_{c}$ (Figure 11). It may be easier to check this by looking at the imaginary part of $\Omega(m,t)$ (cf. right panel).

The oscillatory period of $\Omega(m,t)$ is easily understood in the asymptotic region $t\gg t_{c}$ and $T\searrow 0$. A saddle point analysis yields

[TABLE]

Note that $t\gg t_{c}$, $p_{-}\sim 0$ and $p_{+}\sim\pi$. Due to the asymmetry of the integrand $q\rightarrow-q$ in (II.1) for $p=0$ or $p=\pi$, one can show that $\sigma_{+}(p_{\pm})\sim 0$ so that ${\rm e}^{\sigma_{+}(p_{\pm})}\sim 1$. On the other hand $\epsilon(p_{-})/T\sim(h-4J)/T\ll\epsilon(p_{+})/T$. Thus,

[TABLE]

and we can safely drop the first term in (IV.2). Thus, we are left with a single oscillating term and conclude that the period of the oscillation diverges as $C(h_{c}-h)^{-1}$, for some constant $C$, when $h\rightarrow h_{c}$. Note that $h=h_{c}$ is an exceptional point as the Fermi points pinch the real axis and the above argument is not valid.

V Numerical study in the massless regime $h<h_{c}$

The existence of Fermi points on the real axis makes the evaluation technically more involved than in the massive case. This is due to the fact that the Fermi-point singularity problem and the steepest-descent path problem are coupled. The situation is slightly simpler in the space-like regime from which we start our consideration.

V.1 The space-like regime $t<t_{c}$

We again transform ${\cal E}(\bar{\cal E})$ to a straight contour in the upper (lower) half plane for the hole (particle) momentum if $m$ is not too big. By this deformation, in contrast to the massive case, we inevitably pick up the contribution from the Fermi points. This modification can be treated in a similar manner as the contributions of the other poles, shortly commented on at the end of Section IV.1 on the massive phase. The details are explained in E. We choose a straight contour $[-\pi+i\delta,\pi+i\delta]$ for $p$. Although the choice of $\delta>0$ is in theory arbitrary, in practice one should not take it too small. Table 1 illustrates this with the examples of $m=1$, $t=0$ and $h=0.5J$ for various choice of $\delta$. The rightmost column gives the exactly known static values.

This can be easily understood as the contour passes near the singularities of $\mu$ (at $p_{F}$) or $\bar{\mu}$ (at $-p_{F}$). Hence, too small values of $\delta$ spoil the numerical accuracy.

With suitable choices of $\delta$ we can again adopt the straight line contours well beyond the space-like regime, see Figure 12.

The oscillation frequency in the time-like regime $t>t_{c}$ for small temperatures is given by the effective bandwidth in the Fermionic model, $4J-h$, and goes to zero as $h$ approaches $h_{c}=4J$ as demonstrated in Figure 13.

For larger values of $m$, $\delta$ is determined by the location of the saddle point (12) and one has to the take into account the pole contributions as described in E. Some numerical results will be supplemented later (cf. Section VIII).

V.2 The time-like regime $t>t_{c}$

As in the massive phase, the straight integration contours fail to give reliable results as time evolves. Figure 14 shows an example.

We thus utilize the steepest descent paths. The existence of Fermi points on the real axis makes the situation more involved than in the massive case. At the same time, one must pay attention to the singularities of $\mu(p),\bar{\mu}(q)$. The choice of the paths thus becomes a complicated technical issue, and we decided to summarize it in F.

By applying the prescription presented there, the result are now stable as shown in Figure 15, where we have chosen the same parameters as in Figure 14. The right panel represents the result up to $Jt=80$, and the real part of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{2}(t)\rangle$ decays to the order of $10^{-7}$.

Examples with slightly different parameters are given in Figure 16, showing again smooth behavior until they exhibit sufficient decay.

The fluctuation is enhanced at lower temperatures, where one has to be particularly careful. Independent calculations also suggest this. For example, the authors of Derzhko et al. (2000) calculated $\langle\sigma^{x}_{i}(0)\sigma^{x}_{i+50}(t)\rangle$ in the finite XX-chain (400 sites) using the Pfaffian representation. Their result does not converge numerically for various choices of $i$ around $Jt=30\sim 40$ for $T/J=0.04$ (in the present normalization). The authors claim that this is due to a “boundary effect”. Our results using the Pfaffian representation of the correlation function, which we present in Section VIII, do not show such instabilities. This suggests that the observed instabilities are numerical in nature (a critical step is a numerically stable direct evaluation of the Pfaffian) and not related to boundary effects. In the framework of our QTM approach we are dealing with an infinite size system from the beginning so that boundary effects are certainly absent. Nevertheless, a similar instability is observed at lower temperatures. We examined each factor in (10) and found that the instability mainly comes from the Fredholm determinant (Figure 17).

In order to achieve higher accuracy in the evaluation of the Fredholm determinant, we increase the number $n$ of discretization points. We defer the discussion of the right choice of $n$ to G. With suitable values of $n$, we are able to perform a precise evaluation of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(t)\rangle$ on a longer time scale. An example is given in Figure 18. The right panel presents a zoom of the curves for large $Jt$. This exhibits a slowly decaying oscillating pattern.

VI Even-odd effect

There are similarities and differences between the transverse correlator in the massive and in the massless regime, as we observed above. We comment on one more difference which seems to have been overlooked in past publications: the even-odd effect. By this we mean a drastic suppression of the oscillation amplitude of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(t)\rangle$ for $m$ being odd and for small $h/J$. The plots in Figure 19 present examples in the massless phase.

On the other hand, this even-odd difference is not observed in the massive case. See Figure 20.

Our formula consists of several factors (10). The behavior of $\Omega$ is the same for $m$ odd and $m$ even. We numerically find that in the massless case the oscillation phase of $\Omega(m,t)$ and that of the Fredholm determinant part almost cancel each other and that this results in the monotonous time evolution for odd $m$. By way of contrast the two contributions do not cancel each other for $m$ even or in the massive case.

The even-odd effect in the massless case can be explained straightforwardly within non-linear Luttinger liquid theory. A first step in the derivation of the long time asymptotics of the transverse spin-spin correlation function of the XXZ chain is a Jordan-Wigner transformation $\sigma^{-}_{m}=(-1)^{m}\exp(i\pi\sum_{l<m}c_{l}^{\dagger}c_{l})c_{m}$ onto spinless Fermions with annihilation operator $c_{m}$. Next, the dispersion is linearized around the Fermi points $c_{m}\sim\text{e}^{ik_{F}x}\Psi_{R}(x)+\text{e}^{-ik_{F}x}\Psi_{L}(x)$ with $x=ma$ where $a$ is the lattice constant. After using the standard Bosonization identities, this approach then allows to calculate the low-energy contributions to the transverse spin-structure factor which are located at momentum $ka\sim\pi$ and $ka\sim 2\pi M$ where $M$ is the magnetization per site. Standard Bosonization thus predicts that there are two low-energy contributions to the transverse spin-spin correlation function. At the XX point, the contribution from $ka\sim\pi$ is dominant and is given for small temperatures $T\ll 4J$ by Karimi and Affleck (2011)

[TABLE]

Here $v=4J\sin(k_{F}a)$ is the sound velocity and $\eta^{+}$ denotes a small regularization parameter. Note, in particular, that for $T=0$ and $vt\gg x$ standard Bosonization predicts a decay $\sim 1/\sqrt{vt}$ at long times. Furthermore, there is no part which oscillates in time at this level. To obtain the oscillating contribution one needs to keep—in addition to the Fermi point contributions covered by standard Bosonization—also the saddle point contributions which appear at real momenta $\bar{p}_{\pm}$ for $t>t_{c}$, see eq. (13) 222 Here the momentum has a proper dimension $\bar{p}_{\pm}={p}_{\pm}/a$. Doing so leads to non-linear Luttinger liquid theory. In the XX case, the theory is particularly simple because the saddle point and the Fermi point contributions do not interact with each other. We now use the extended ansatz $c_{m}\sim\text{e}^{ik_{F}x}\Psi_{R}(x)+\text{e}^{-ik_{F}x}\Psi_{L}(x)+\text{e}^{i\bar{p}_{+}x}d(x)$ where $d$ denotes a particle near the saddle point. It is then a fairly straightforward calculation to show that the additional saddle point contribution is Karimi and Affleck (2011)

[TABLE]

We are interested here in understanding the asymptotic behavior in the time-like regime for $vt\gg x$ and $T\ll 4J$. In this limit we have $\bar{p}_{+}\,a\to\pi$ and we can approximate the propagator for the high-energy particle by $\langle d(0,0)d^{\dagger}(x,t)\rangle\propto\text{e}^{i(4J-h)t}/\sqrt{vt}$. Putting together the two contributions (17) and (18) for this case we arrive at

[TABLE]

with some unknown amplitude $A$. In the zero temperature limit, in particular, eq. (19) predicts that the correlation function asymptotically consists of a uniform part which decays as $1/\sqrt{vt}$ and part which decays as $1/vt$ and oscillates with frequency $4J-h$.

Crucially, the oscillating part contains a factor $\cos(k_{F}x)$. This means that for $h\to 0$ ($k_{F}\,a\to\pi/2$) the oscillating contribution is absent if the spatial distance $m$ in the two-point correlation function is odd!

To check the predictions of Eq. (19) we present in Figure 21 low-temperature numerical data for the auto- and nearest-neighbor correlation function at $h/J=0$ (left panel) and $h/J=2$ (right panel).

¿From the log-log plot it is clear that for $tT\ll 1$ the correlation functions decay following a power law in time for all cases. For $h/J=0$ ($k_{F}\,a=\pi/2$) the autocorrelation shows oscillations while the nearest-neighbor correlation function does not whereas both are oscillating for $h/J=2$. A comparison with eq. (19) using the amplitudes of the uniform and the oscillating part as fitting parameters confirms that the long time behavior is described by the non-linear Luttinger liquid result and that the factor $\cos(k_{F}x)$ is indeed the reason for the observed even-odd effect at zero magnetic field. Finally, we note that higher order corrections to eq. (19) do exist, i.e. terms which decay faster than $1/vt$ at $T=0$. Such terms are responsible for the remaining small oscillations observed in the nearest-neighbor correlation function at zero field.

VII Comparison with asymptotic formulae

Starting from a slightly different Fredholm determinant Its et al. Its et al. (1993) reformulated the problem of the asymptotic analysis of the transverse correlation function in terms of a matrix Riemann Hilbert problem. They obtained a closed expression for the leading terms of the asymptotic expansion of the transverse correlation function in $t\rightarrow\infty$ and $m\rightarrow\infty$ with a fixed direction $\phi$ s.t. $\cot\phi=m/4Jt$ at arbitrary $T$ and $h<h_{c}$. The space-like regime (VII.1) corresponds to $0\leq\phi<\pi/4$ while the time-like one (23) does to $\pi/4<\phi\leq\pi/2$. Subsequently their method was applied to the massive case $h>h_{c}$ in a PhD thesis by Jie Jie (1998). An interesting feature of our novel thermal form factor series (29) is that it is more suitable for the long time, large distance analysis in the space-like regime than the Fredholm determinant representation of Its et al. Its et al. (1993) in that the asymptotics can be easily extracted from the first term in the series. For this reason we can obtain the constant term in the asymptotic expansion which was heretofore unknown in the massless phase. In this section, we provide a quantitative comparison of the available asymptotic formulae with our numerical results.

VII.1 Space-like regime in the massless phase $h<h_{c}$

Its et al. Its et al. (1993) use a Hamiltonian $J=-1,h\rightarrow 2h$ in (1) and provide formulae for the conjugate transversal correlation function

[TABLE]

By a unitary transformation and an adaption of parameters, it is related to ours by

[TABLE]

In the space-like regime of the massless phase their asymptotic formula translates into

[TABLE]

after properly recovering the exchange coupling $J$. The constant $C(T,h)$ was not explicitly calculated in Its et al. (1993). Within our framework it is not hard to reproduce this formula and to obtain an explicit expression for $C(T,h)$ from eq. (29),

[TABLE]

where $\lambda_{F}^{-}$ is an inverse image of $-p_{F}$, $p(\lambda_{F}^{-})=-p_{F}$ in (26). See Göhmann et al. (2019b) for details. We supplement convincing numerical support in Figure 22.

We have to adopt the steepest descent paths even in the space like regime for larger values of $m$ in order to achieve sufficient accuracy.

We assumed that the other poles (of $\mu$) are sufficiently away from $\lambda_{F}^{-}$ in the above derivation. When $T\searrow 0$, however, they accumulate towards the Fermi points, and we inevitably have contributions from them which are neglected in the above as exponentially small corrections. Their estimation is an interesting problem which we hope to discuss in the future.

VII.2 Time-like regime in the massless phase $h<h_{c}$

The asymptotic formula of Its et al. (1993) becomes more involved in the time-like regime,

[TABLE]

with

[TABLE]

and

[TABLE]

The authors of Its et al. (1993) did not present the explicit form of $A_{0}(T,h)$ ($C$ in their notation) but commented that the higher-order correction modifies $A_{0}(T,h)\rightarrow A_{0}(T,h)(1+c(t,m))$, where $c(t,m)\sim t^{-\frac{1}{2}}$. We again use (21) and re-introduce $J$. As $I_{m}$ is real-valued for positive $t$, it is natural to consider the ratio

[TABLE]

In view of asymptotic analysis, $r(m,t,h)$ depends on $m$ only through $m/t$.

Based on our numerics (with $n=1536$) we propose

Conjecture 1**.**

For fixed $T,h,m$ and $t\gg 1$, the ratio $r(m,t,h)$ consists of both non-oscillating and oscillating parts,

[TABLE]

The non-oscillating part $r_{\rm non-osc}(m,t,h)$ includes $A_{0}(T,h)$. The period $2\pi/\omega$ of the oscillating part $r_{\rm osc}(m,t,h)$ behaves as $\omega\sim h_{c}+h$ as $h\rightarrow-h_{c}$, and its amplitude is also expanded by powers of $t^{-\frac{1}{2}}$.

We, however, do not have any estimate of the time scale when the above asymptotic form becomes valid. For example, let us look at $r(50,t,h)$ for $T=0.05J,h=0.5J$ with $n=1536$ in Figure 23.

The average value does not seem to saturate, which suggests that $\langle\sigma^{-}_{1}(0)\sigma^{+}_{51}(t)\rangle_{-h,T}$ is yet to reach the region described by the asymptotic formula if $t<100$. There is, however, a limitation on the range of $t$, for numerical accuracy. We thus need to find other ways than dealing with $t\gg 100$. Intuitively, the smaller $m$ is, the shorter we have to wait until $r(m,t,h)$ reaches “equilibrium”. Thus, we focus on the extreme case $m=0$ (the auto-correlation). To be precise, this corresponds to the direction $\phi=0$ and the direct application of the formula (23) may need justification. We however assume its validity in the following.

The real parts of $r(0,t,h)$ are plotted in Figure 24 for $T/J=0.1,0.3,0.5$ with $h/J=0.05$.

While \operatorname{Re}\bigl{(}r_{\rm non-osc}(0,t,h)\bigr{)} is still slowly decreasing for $T/J=0.1$, it already seems to approach “equilibrium” for $T/J=0.3,0.5$ around $Jt\sim 40$. The amplitude of the oscillating part \operatorname{Re}\bigl{(}r_{\rm osc}(0,t,h)\bigr{)} seems to be very slowly decreasing. In order to examine this quantitatively, we fit $r_{\rm non-osc}(0,t,h)$ for $m=0,h=0.05J,T=0.3J$ using the data for $40<tJ<70$,

[TABLE]

The sum of the first two terms in the right-hand side is identified with $A_{0}(0.3J,0.05J)$. The sub-leading terms agree with the remark by Its et al. Taking this for granted, the real part of $r_{\rm osc}(0,t,h)$ is estimated as in Figure 25.

The plot suggests that the amplitude of $\sqrt{Jt}r_{\rm osc}(0,t,h)$ stays constant, that is, $r_{\rm osc}(0,t,h)\sim t^{-\frac{1}{2}}$. This is again consistent with the estimate by Its et al..

There are remarks. One notices that \operatorname{Re}\bigl{(}r(0,t,h)\bigr{)} for $T/J=0.5$ shows a discontinuity at $T/J\sim 70$. This is a numerical artifact due to inaccuracy: $|\langle\sigma^{-}_{1}(0)\sigma^{+}_{1}(t)\rangle|$ already reaches $O(10^{-25})$. Thus, it is too difficult to have precise control over the numerics. The oscillatory part remains observable after sufficiently long time (except for $h=h_{c}$). The origin of the oscillation can be attributed to that of $\Omega(m,t)$. The comparison of the real parts of $r(0,t,-h)$ and $\Omega(0,t)$ is given in Figure 26 for $T/J=0.2,h/J=2$.

We can easily check that their frequencies coincide. One expects from eq. (10) that $\langle\sigma^{-}_{1}(0)\sigma^{+}_{1}(t)\rangle$ is proportional to $\Omega(0,t)$, but it is not so simple. The inverse power of $\Omega$ in the kernel of the Fredholm determinant partially cancels the oscillation, and this may contribute to the non-oscillating part $A_{0}$. The center of oscillation of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{1}(t)\rangle$ is thus different from zero, while it is zero for $\Omega(0,t)$.

VII.3 Consistency with the result from nonlinear Luttinger Liquid theory

Finally, we comment on the relation between the above analysis and the results from the continuum theory, (17) and (19). We shall neglect logarithmic corrections in $t$ and consider only the leading $T\rightarrow 0$ behavior. Let us start from the space-like case, $m,t\rightarrow\infty$ and $t<t_{c}$. One applies a Sommerfeld type argument to the integral in (VII.1) and finds

[TABLE]

This is consistent with (17) for $T$ small but finite and $\pi Tx/v\gg 1$. Next consider the time-like case $m,t\rightarrow\infty$ and $t>t_{c}$. The integral $I_{m}$ in (VII.2) has different forms for $p_{-}<p_{F}$ and for $p_{F}<p_{-}$,333The symbol $p_{-}$ is defined in (13).

[TABLE]

Thus if $p_{-}<p_{F}$, together with Conjecture 1 and relation (21), we conclude that

[TABLE]

One identifies the terms on the right hand side with the dominant and sub-dominant terms in (19) if $tT\gg 1$. This condition is consistent with the asymptotic analysis as $T$ is small but finite. In view of (19), the third and fourth terms in (25) look irrelevant and may be attributed to the weak transient behavior which still exists. We however do not have a conclusive evidence.

The above result suggests that the meaning of the space-like and time-like regimes undergoes a change in the low-temperature limit, or, in other words, that the ‘light cone’ needs to be redefined. At elevated temperatures, it is defined by $m/t=4J$, while in the low-$T$ limit, it is defined by the Fermi velocity, $m/t=v(\leq 4J)$. This means that the space-like regime becomes enlarged in the low-T limit and that it becomes more enlarged if the magnetic field is higher.

VII.4 Time-like regime in the massive case $h>h_{c}$

A similar long time, large distance asymptotic formula for the massive case $h>h_{c}$ was proposed by Jie Jie (1998). Set

[TABLE]

Then Jie argued that444An obvious typo in the second term ($(-Jt)^{i\nu_{1}}$ instead of $(-Jt)^{i\nu_{2}}$) is corrected.

[TABLE]

where \alpha=\pi-\arcsin\bigl{(}\frac{m}{4tJ}\bigr{)} and $I_{m}(t)$ is defined in (VII.2). The coefficient ${\mathfrak{C}}$ is a smooth function of $\alpha,T,m$ and $h$. The explicit forms of $u_{-},v_{-}$ and ${\mathfrak{C}}$ are too complicated to be reproduced here. We only remark that $|v_{-}|\gg|u_{-}|$ as $T\ll 1$. Thus, only the second term survives in this limit, and it results the long-period oscillation $\omega\sim h-4J$ observed numerically. The ratio

[TABLE]

is plotted in Figure 27. It seems to reach a constant value as time evolves.

VIII Comparison with the Pfaffian representation

In this section we compare our results with results obtained by more conventional methods. Among them, the Pfaffian representation of the correlation function is the most well-established for the dynamics of the XX model. It deals with an open chain of $L$ sites. As the translational invariance is broken, we rather consider $\langle\sigma_{j}^{-}(0)\sigma_{j+m}^{+}(t)\rangle$ instead of (3) and typically center the two-point correlation function at $L/2$. After a simple calculation using the Fermion algebra, one represents $\langle\sigma_{j}^{-}(0)\sigma_{j+m}^{+}(t)\rangle$ by the Pfaffian of an anti-symmetric $(2j+2m+2)\times(2j+2m+2)$ dimensional matrix, see e.g. McCoy et al. (1971); Derzhko et al. (2000) for details.

For small $m$ and $t$ the two results, one obtained by our formalism and the other one by the Pfaffian method,555We refer the numerical calculation based on the Pfaffian representation of the correlation function as “the Pfaffian method” for short. match perfectly. See Figure 28 for the autocorrelation function for $m=0$, $h/J=0.05$, $T/J=0.1$ and $L=128$.

One of the advantages of the present approach in comparison with the Pfaffian method lies in the fact that the distance $m$ appears as a mere parameter. This allows us to deal with large $m$ immediately. For illustration, the data for $m=100,200$ and $300$, obtained from eq. (10) are plotted in Figure 29 for $T/J=0.05,h/J=0.5$.666 As $m$ is large enough, we use the steepest descent method also in the space-like regime. It becomes inaccurate around $t_{c}=m/4J$ and some points are omitted near $t_{c}$ in the plots.

As indicated above, the Pfaffian formulation deals with a finite-size system. In order to eliminate the boundary effects, one has to choose the system size $L$, the spin distance $m$ and the time $t$ so that $L/2-m/2>vt$ is satisfied. Here $v$ stands for the sound velocity. Thus, for $m$ and $t$ becoming larger, one has to deal with a larger-size matrix. It is fair to note that one also needs to increase the number of discretization points $n$ in the evaluation of the Fredholm determinant in order to keep the numerical precision for larger $t$. An example is discussed in G. With these tunings and using a stable direct evaluation of the Pfaffian (instead of calculating it as the square root of the determinant), the two results coincide also for large $m$ and $t$. An example with $m=200,T/J=0.05$ and $h/J=0.5$ is shown in Figure 30.

We stress, however, that the numerical error in the present scheme solely originates from the discretization and should be distinguished form the error due to the finite size effect explained above. In the former case the error was estimated in Bornemann (2010). For the type of kernel occurring in our Fredholm determinant the error decreases algebraically with the number of discretization points.

IX Summary and conclusion

In this communication we revisited the equilibrium dynamics of the XX chain. We employed a new scheme for the exact evaluation of correlation functions, based on the QTM and on a thermal form factor expansion. We rewrote the thermal form factor series for the transversal correlation function as an explicit factor times a Fredholm determinant. This representation was the starting point of our numerical analysis, for which we utilized a direct discretization of the determinant as suggested in Bornemann (2010). We demonstrated that the method yields high-precision results for the transverse correlation function for a wide range of time, distance, temperature and magnetic field, if we properly use the freedom in the choice of the integration contour involved in the definition of the Fredholm determinant. The algorithm works for long times and large distances, far in the asymptotic regime, where we could confirm all existing analytic results Its et al. (1993); Jie (1998) and could even provide an estimation of the next order correction from our numerical data. We also checked our results by comparing with a numerical evaluation of the correlation function using the Pfaffian representation. Furthermore, we were able to explain even-odd effects in the observed oscillations in the massless time-like regime by non-linear Luttinger liquid theory.

We naturally expect, based on the success of the numerical study here, that the novel Fredholm determinant representation also provides an appropriate starting point for the analytic study of the XX chain. Although we have just started the investigation in this direction, there is already evidence that this is indeed the case Göhmann et al. (2019a, b).

A most important issue, from our point of view, is a possible extension to a truly interacting system, the XXZ chain. A single Fredholm determinant representation does not seem to exist in this case. But the thermal form factor series exists in the interacting case as well and represents the transversal correlation function as a series of multiple integrals having certain analogies with a Fredholm determinant Göhmann et al. (2017). Our numerical investigation suggests that the first few terms of this series may already provide accurate estimates in appropriate limiting cases. We hope to report concrete results in future publications.

Acknowledgments. F. Göhmann acknowledges financial support by the (DFG) in the framework of the research group FOR 2316 ‘Correlations in Integrable Quantum Many-Body Systems’. K. K. Kozlowski is supported by the CNRS and by the ‘Projet international de coopération scientifique No. PICS07877’: Fonctions de corrélations dynamiques dans la chaîne XXZ à température finie, Allemagne, 2018-2020. J. Sirker acknowledges support by the DFG in the framework of the research group FOR 2316 ‘Correlations in Integrable Quantum Many-Body Systems’ as well as support by the National Science and Engineering Research Council (NSERC) of Canada. J. Suzuki is grateful for support by a JSPS Grant-in-Aid for Scientific Research (C) No. 15K05208, No. 18K03452 and by a JSPS Grant-in-Aid for Scientific Research (B) No. 18H01141.

Appendix A A previous result

In this appendix we recall a previous result, obtained in Göhmann et al. (2017), from which we have derived the theorem in Section II.

In Göhmann et al. (2017) we proposed a method for calculating dynamical correlation functions at finite temperature in integrable lattice models of Yang-Baxter type. The method is based on an expansion of the correlation functions as a series over matrix elements of a time-dependent quantum transfer matrix rather than the Hamiltonian. Staying with the example of the transverse correlation functions of the XX chain, the series is obtained, when we apply, in a slightly sophisticated manner, respecting the integrability of the original quantum chain Klümper (1993); Suzuki et al. (1990); Klümper (1992), a Suzuki-Trotter decomposition Suzuki (1985) to the exponential factors in (3). As a result, the correlation function is represented as a (properly normalized) partition function of an alternating vertex model acting on a fictitious space of size $2N+2$, where $N$ is the Trotter number. This partition function can be interpreted as a trace of a product of many staggered column-to-column transfer matrices, the quantum transfer matrices, with two local insertions corresponding two the two local operators $\sigma^{-}$, $\sigma^{+}$, whose correlation function is considered. Using the ‘solution of the quantum inverse problem’ corresponding to the staggered monodromy matrices appearing in the columns and expanding in terms of the eigenstates of the QTM the thermal form factor series is obtained.

In case of the XX chains all matrix elements and correlation lengths appearing in the series can be calculated explicitly and the sums over classes of excitations (one hole, two holes $+$ one particle, three holes $+$ two particles …) can be turned into integrals. The integrands are composed of a number of functions characteristic of the XX chain. In the first place there are momentum and energy of the one-particle excitations as functions of the rapidity variable,

[TABLE]

The other functions needed in order to write the form factor series are

[TABLE]

and

[TABLE]

The contours $\cal C$ and ${\cal C}^{\prime}$ in the definition of $\cal A$ are such that $\cal C$ tightly encloses ${\cal C}^{\prime}$, and $\cal C$ itself simply surrounds the strip $|\operatorname{Im}(\lambda)|<\frac{\pi}{4}$ counterclockwise with two infinitesimal deformations at the Fermi rapidities $\lambda_{F}^{\pm}$ defined by $\epsilon(\lambda_{F}^{\pm})=0$, $-\pi/2<\Im(\lambda_{F}^{\pm})<\pi/2$, $\Re(\lambda_{F}^{-})<0<\Re(\lambda_{F}^{+})$.

Employing all of the above notation the thermal form factor series for the transverse correlation function (3) can be written in the form

[TABLE]

Here we are using a contour ${\cal C}_{\rm out}$ which tightly encloses ${\cal C}$. Note that as compared to the formula given in Göhmann et al. (2017) we have introduced a few minor simplifications by calculating two of the integrals in our original formula explicitly (one is zero the other one a constant) and by absorbing a constant term into the definition of the function $\Phi$.

Appendix B The analytic properties of $\mu$

We first consider the massive case $h>h_{c}$ and assume that $\operatorname{Im}(p)>0$. We take the integral over a closed rectangular contour ${{\cal C}_{+}}$ (of infinite height and width $\pi$), including the real axis as the bottom part,

[TABLE]

Note that the integrand on the upper edge of the rectangle reduces to a constant $\frac{i}{2}$ and the contributions from left and right sides cancel each other. We thus have

[TABLE]

On the other hand, one can evaluate the integral directly by summing up the contributions from the branch cuts connecting $q^{u}_{2j}$ and $q^{u}_{2j+1}$ and the pole at $q=p$,

[TABLE]

Thus, we conclude

[TABLE]

By remembering (15), one verifies that the double zeros and double poles at $q^{u}_{2k+1}$ cancel and that the single poles at $q^{u}_{2k}$ survive if $\operatorname{Im}(p)>0$.

The argument for the case $\operatorname{Im}(p)<0$ goes almost in parallel. This time we consider a closed rectangular contour in the lower half plane. Using a similar reasoning as above, we obtain,

[TABLE]

Thus, for $\operatorname{Im}(p)<0$, $\mu(p)$ possesses single zeros at $q^{d}_{2k}$ and double poles at $q^{d}_{2k+1}$. The same conclusion can be drawn from (30), though it is derived under the assumption that $\operatorname{Im}(p)>0$.

Next we consider the massless regime ($h<h_{c}$) and assume that $\operatorname{Im}(p)\geq 0$. We again adopt a contour in the upper half plane which includes ${\cal E}$ as a part. As ${\cal E}$ encircles the Fermi point $p_{F}$ from below, we obtain

[TABLE]

thus,

[TABLE]

From this we easily see that $\mu(p)$ has single poles at $q^{\ell}_{-2j}\,(j\geq 1)$ and $q^{r}_{-2j}\,(j\geq 0)$ if $\operatorname{Im}(p)\geq 0$. The expression is valid also for $\operatorname{Im}(p)<0$. Thus, we see that $\mu(p)$ has single zeros at $q^{\ell}_{2j}\,(j\geq 0)$ and $q^{r}_{2j}\,(j\geq 1)$, while it possesses double zeros at $q^{\ell}_{2j+1}\,(j\geq 0)$ and $q^{r}_{2j+1}\,(j\geq 0)$.

Appendix C Static correlations for small $m$ at finite $T$

The static correlation functions of the $S=\frac{1}{2}$ XXZ model were studied in Boos et al. (2008) in the framework of a QTM approach. Explicit formulas for the reduced density matrix of a short chain segment in an arbitrary magnetic field at arbitrary temperature were obtained. They can be used to express the two-point functions in terms of two fundamental functions $\omega(x,y)$ and $\omega^{\prime}(x,y)$ (note that the prime does not mean a derivative). See eqs. (59) and (62) in Boos et al. (2008) for their definition. Their arguments are introduced as inhomogeneities which must be sent to zero in the end. The homogeneous (physical) limit produces derivatives with respect to the first and the second argument, which will be denoted by, e.g., $\omega_{x}=\partial_{x}\omega(x,y)|_{x,y\rightarrow 0}$, $\omega_{y}=\partial_{y}\omega(x,y)|_{x,y\rightarrow 0}$ and so on. The formulas behave non-trivially when the deformation parameter $q$ that parameterizes the anisotropy of the XXZ chain as $\Delta=(q+q^{-1})/2$ approaches a root of unity, $q^{n}=1$. The XX chain corresponds to $q=i$, viz. $\Delta=0$. It is conjectured that the expression for the transversal two-point function has a singularity of removable type, whenever one evaluates correlators points that are at distance $km$, for some integer $k$, if $q^{m}=-1$. In the present case we conjecture that l’Hospital’s rule must be applied for the computation of $\langle\sigma_{1}^{x}\sigma_{m+1}^{x}\rangle$ if $m=2,4,6,\dots$.

For $m=2$, for instance,

[TABLE]

where $\eta$ is such that $\Delta=\operatorname{ch}(\eta)$. The first and the third terms are obviously singular at the XX point $\eta=\frac{i\pi}{2}$. We thus introduce a small deviation $\gamma=-i\eta=\frac{\pi}{2}-\varepsilon$ and obtain by straightforward expansion with respect to $\varepsilon$,

[TABLE]

We have verified numerically that the first term vanishes, i.e., the apparent singularities cancel each other. We, however, need to evaluate derivatives with respect to the anisotropy parameter, which adds an extra elaboration.

For generic $\gamma$ the function $\omega$ is expressible in terms of auxiliary functions $\mathfrak{b}(x),\bar{\mathfrak{b}}(x),g_{\mu}^{(\pm)},{g^{\prime}}_{\mu}^{(\pm)}$ defined in Boos et al. (2008),

[TABLE]

A similar expression can be written for $\omega^{\prime}$. In order to evaluate $\omega$ and $\omega^{\prime}$, one needs to solve non-linear integral equations (cf. eqs. (52)-(55) in Boos et al. (2008)) for the auxiliary functions. The XX model is exceptional, however, in that one of the integration kernels, $F$, vanishes. Then $\mathfrak{b}(x),\bar{\mathfrak{b}}(x),g_{\mu}^{(\pm)}$ can be obtained explicitly, and ${g^{\prime}}_{\mu}^{(\pm)}$ are calculated by taking convolutions of $g_{\mu}^{(\pm)}$ with appropriate kernels. Even in evaluating the derivatives with respect to $\gamma$, one does not have to solve the integral equations, but only needs to perform integrations over already known auxiliary functions. For example, one needs to evaluate $\partial_{\gamma}\mathfrak{b}^{-1}=-\mathfrak{b}^{-1}\partial_{\gamma}\ln\mathfrak{b}$. This is evaluated from,

[TABLE]

where $D_{b}$ denotes a known function and $\delta$ is a small quantity introduced for technical reasons. The parameter $\gamma$ must be set equal to $\frac{\pi}{2}$ after taking the derivative on the right-hand side. We thus need only known functions in order to evaluate $\partial_{\gamma}\mathfrak{b}^{-1}$. In a similar manner, one evaluates the $\gamma$ derivatives of $g_{\mu}^{(\pm)}$, ${g^{\prime}}_{\mu}^{(\pm)}$ and then obtains $\langle\sigma_{1}^{x}\sigma_{3}^{x}\rangle$.

On the other hand, since there are no apparent singularities for $m=1$ or $3$, we can take the limits directly in these cases and obtain

[TABLE]

We compare the values of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(0)\rangle$ obtained from our Fredholm determinant representation and those obtained from the above exact formulas777 Half of the values of (31), (32), (33) for $m=2,1,3$ reflecting the symmetry between $xx$ and $yy$ correlators. The static autocorrelation is obtained from the magnetization $\mathfrak{m}(T,h)=\int_{-\pi}^{\pi}\frac{dp}{4\pi}\operatorname{th}\frac{\varepsilon(p)}{2T}$ by $\frac{1}{2}-\mathfrak{m}(T,h)$. in Table 2 (Table 3) for $h=0.1J$ $(<h_{c})$ ($h=4.1J$ $(>h_{c})$). Their agreement up to a reasonable number of digits assures the validity of our formulation in the static limit for small $m$.

Appendix D Static correlations for larger $m$ at $T\ll 1$

While exact correlations are not available for larger $m$ at finite temperatures, the ground state correlation function is obtained explicitly for $h=0$ Kitanine et al. (2002). We find that it is neatly expressed in terms of Barnes’ $G$ function,

[TABLE]

where

[TABLE]

Our formulation has numerical problems in the limit $T,h\rightarrow 0$. Nevertheless, we tried to check the consistency at larger $m$ by taking numerical limits. For this purpose we set $h=0.01$ and extrapolated the zero temperature values by fitting data for $T=0.03\sim 0.06J$ using quadratic curves (see Table 4). In spite of the above difference in details, we find that the results qualitatively agree with the exact values.

Appendix E The modification due to poles

The integration contours ${\cal E},\bar{{\cal E}}$ are originally located near the real axis. As discussed in the main text, we need to shift ${\cal E}$ to $[-\pi,\pi]+i\delta$ and $\bar{{\cal E}}$ to $[-\pi,\pi]-i\delta$, so that they pass through the saddle points, especially when $m$ is large in the space-like regime. In the course of this deformation, the paths cross poles and we have to take account of these contributions and modify (10). This can be easily done by following appendix C of Dugave et al. (2016).

We denote the sets of poles of $\mu(p)$ and $\bar{\mu}(q)$, crossed by the contours, by $\{s_{j}\}_{j=1}^{m}$ and $\{\bar{s}_{j}\}_{j=1}^{m}$, respectively. These sets may contain, in particular, the Fermi points, namely, $p_{F}\in\{s_{j}\}$ and $-p_{F}\in\{\bar{s}_{j}\}$ in the massless case.

The shift of ${\cal E}$ modifies the functions by analytic continuation,

[TABLE]

Here we have used the fact that $q$ is in the lower half plane. Note that $\tilde{V}(q,q^{\prime})$ is also modified due to the expression (11).

We also have to take account of the result of the shift of $\bar{{\cal E}}$. The analytic continuation modifies the Fredholm determinant,

[TABLE]

where we introduced

[TABLE]

In order to simplify this, we define the resolvent Kernel $R(q,q^{\prime})$ by,

[TABLE]

By the simple determinant identity

[TABLE]

one concludes that

[TABLE]

where $(\hat{K}f)(q):=(\hat{P}f)(q)-(\hat{V}f)(q)$.

Appendix F Choice of the contour: massless and time-like regime

For $t$ slightly greater than $t_{c}$, $p_{\pm}$ is near $\pi/2$ from (13) and $p_{F}$ is not necessarily between $p_{\pm}$. In this short time region, however, the straight contours work equally well as in the space-like regime and do not have to be deformed. Since $p_{+}$ (respectively  $p_{-}$) moves towards $\pi$ (respectively [math]), $p_{F}$ soon lies between $p_{\pm}$ (we assume this below) and we need to think about possible deformations.

We do not take the paths depicted in Figure 7 but use a steepest descent path for the hole variable $p$ and a straight line just below the real axis for the particle variable $q$, depicted explicitly in Figure 31. We list reasons for the choice below.

We did not deform the contour for the $q$ variable into the upper half plane, suggested by the steepest descent path, for two reasons. First, as already mentioned, the pole contribution of $\varphi(p,q)$ at $p=q$ brings a divergent contribution for $t\gg 1$. Second, $\bar{\mu}$ behaves irregularly in a region sandwiched by a double pole (in the upper half plane, closest to the real axis) and a zero (at $q=p_{F}$) of $\bar{\mu}$. 2. 2.

On the other hand, we do deform the contour for the $p$ variable into the lower half plane, but in a specific manner. This is due to the fact that there are two closely positioned poles for $\mu$ (the black and the red triangles in Figure 31). In the region sandwiched by these, $\mu$ behaves quite unstable. We thus decided to deform he contour over the double pole of $\mu$ (the red triangle). Note that at the blue triangle (= a pole of $\bar{\mu}$), $\mu$ is null, and around the point it behaves smoothly. 3. 3.

We can shift the blue contour below the red triangle, in principle. This may be sometimes a nice choice as long as $t$ is not so large, since the pole contribution of $\varphi(p,q)$ at $p=q$ behaves irregularly between the black and the red triangle. In a later stage, however, the phase factor ${\rm e}^{-tu_{m,t}(q)}$ shows a divergent behavior if the contour is too far away from the real axis.

We define the displacements $\delta_{1}$ and $\delta_{2}$ as in Figure 31. For reasons 2 and 3, we take $\delta_{1}>\delta_{2}$, as the opposite choice $\delta_{1}<\delta_{2}$ would ruin the stability of the long time runs by reason 3.

The most subtle problem is posed by the question as to how to stabilize both factors ${\rm e}^{-tu_{m,t}(q)}$ and $\mu(q_{i})\mu(q_{j})\bar{\mu}(q_{j}){\rm e}^{tu_{m,t}(q_{i})+tu_{m,t}(q_{j})}$ simultaneously. The former is the phase factor generic for the particles and diverging for $\operatorname{Im}(q)<0$. This factor is multiplied by ${\rm e}^{-tu_{m,t}(p)}$ for which an optimal path is already chosen and it suppresses the divergent behavior. Thus, we conclude that the former component is already harmless.

The latter comes from the $p=q$ pole of $\varphi(p,q)$. The naive answer is to use the blue contour just below the real axis, like in the massive case. This does not necessarily give a stable result: a small $\delta_{2}$ causes divergent behavior of $\mu(q)$. We pay attention to the fact that the product $\bar{\mu}(q_{j})\mu(q_{j})$ remains finite for $q_{j}\sim p_{F}$. We then tune $\delta_{2}$ so that $|\mu(q){\rm e}^{2tu_{m,t}(q)}|\sim 1$ for $q\sim p_{F}-i\delta_{2}$ is satisfied for each given $t$. As a consequence $\delta_{2}$ gets smaller for larger $t$. Therefore the calculation fails eventually at a very late stage when $\delta_{2}$ turns infinitesimally small. Empirically, the correlation functions become too small to be detected before this limitation approaches. Thus, practically, the present method successfully yields a stable calculation for a reasonably long time range as demonstrated in the main text.

Appendix G The choice of $n$

We present results of an experimental study on the choice of the dimensions $n$ of the discretized Fredholm determinant.

We start from proposing the following:

Conjecture 2**.**

Both real and imaginary part of $(-1)^{m}\langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(t)\rangle$ stay positive in the massless time-like regime.

This conjecture is motivated by the empirical fact that the numerical result becomes unstable after the real (or imaginary) part of $(-1)^{m}\langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(t)\rangle$ crosses zero for $t>0$.

Thus, we increase $n$ until the positivity of both parts of $(-1)^{m}\langle\sigma^{-}_{1}(0)\sigma^{+}_{m+1}(t)\rangle$ is satisfied and until the result converges (with variation of $n$). In the following we restrict ourselves to the parameters $m=50,T=0.05J,h=0.5J$ and increase the value of $n$.

Especially, we take a closer look at the late stage and compare the real parts with $n=280$, $400$ and $512$. See Figure 32 (left). The real part become negative for $n=280,400$. The case with $n=280$ shows abrupt change in the oscillation amplitude, and the amplitude is not decreasing, but seems increasing for $n=400$. Both cases seem unreasonable. The case with $n=512$ seems well-behaved. The real and the imaginary part of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{51}(t)\rangle$ is plotted in Figure 32 (right) with $n=512$. It does not show any particular fluctuation around $20<Jt<40$ any longer.

We further increase the values of $n$. Figure 33 supplements the real parts of $\langle\sigma^{-}_{1}(0)\sigma^{+}_{51}(t)\rangle$ with $n=1024$, $1280$ and $1536$. One can hardly see the difference. The zoom-in in the right plot convinces us that the results with $1280$ and $1536$ have already reached numerical convergence. Thus, we conclude that the choice $n=1536$ is a safe choice.

We remark that when $T$ is higher, we do not have to take such large values of $n$. The correlations will anyway decay to very small values $\sim 10^{-7}$ quickly, even with relatively small $n$ (like $n=128\sim 256$).