The equilibrium dynamics of the XX chain revisited

The equilibrium dynamics of the XX chain revisited

Frank Göhmann, Karol K. Kozlowski, Jesko Sirker and Junji Suzuki
Fakultät für Mathematik und Naturwissenschaften,
Bergische Universität Wuppertal, 42097 Wuppertal, Germany
Univ Lyon, ENS de Lyon, Univ Claude Bernard,
CNRS, Laboratoire de Physique, F-69342 Lyon, France
Department of Physics & Astronomy, University of Manitoba,
Winnipeg, Manitoba, Canada R3T 2N2
Department of Physics, Faculty of Science, Shizuoka University,
Ohya 836, Suruga, Shizuoka, Japan
June 2019
Abstract

The equilibrium dynamics of the spin- 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.

1 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 [1]. 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 [2]. The evaluation of the two- and four-spinon contributions to the dynamical structure factor of the massive XXZ chain [3, 5] and the four spinon ones in the XXX limit [4] were important outcomes of this method.

A complementary approach for 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 [6] and a special determinant formula for the scalar product of on-shell and off-shell Bethe vectors [7]. The algebraic Bethe ansatz approach successfully confirmed [8] 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 two- and multi-point functions in the thermodynamic limit [20, 21, 22, 23, 24, 25, 26, 27, 28].

In a parallel development, field theoretical approaches have been extended to include non-linearities in the spectrum [29, 30]. 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 [31, 32, 33, 34]. Particular relevant in the present context are results for the asymptotics of the dynamical longitudinal spin-spin correlation function at finite temperatures [35] and the transverse spin-spin structure factor [36] for the XXZ chain. 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.

It has been observed that the simplest multiple integrals representing the reduced density matrix of short sub-segments of the XXZ chain factorize [10, 11, 12, 13, 14]. 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 [15, 16, 17]. While the structure turned out to be relevant for the analysis of the short-range correlation functions [19, 18], 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[21, 25] of the predictions of conformal field theory [37] 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 [38], 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 [39, 40, 41]. The multiple integrals factorize even at finite temperatures [42, 43, 44], which finds its explanation again in the hidden Grassmann structure [45]. 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 [46, 47, 48]. 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 [48]. 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 [49] 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 [49] with the thermal form factor expansion of [46], we recently proposed [50] 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 [51] for the longitudinal equilibrium correlation function of the spin- XX chain. We also obtained a novel thermal form factor series for the transversal two-point correlation function which will be analyzed in this work.

It may seem rather amazing that for a as simple model as the XX chain, whose Hamiltonian can be expressed in terms of spin-less Fermions [52], 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 [53, 54] can be applied. Results at finite temperatures, that were obtained within approaches based on the use of Fermions, comprise asymptotic formulae for high temperatures [54, 55] 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 [56, 57].

A completely different route to the evaluation of the transverse dynamical two-point functions was taken in [59]. Following the strategy of [60] the authors of [59] 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 belongs to the integrable type. This formulation was then directly suitable for an asymptotic analysis at long times and large distances by means of an associated matrix Riemann-Hilbert problem amenable to a ‘non-linear steepest descent method’ [69]. In [58] it was shown that the long-time, large-distance behavior of the two-point function for a fixed ratio in a weak magnetic field () goes like . The explicit expressions for and were obtained in the space-like () and time-like () asymptotic regimes. The result was subsequently extended to strong magnetic field () [61].

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 that was derived in a previous communication [50]. In [62, 63] 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 the amplitude in the space-like regime at weak magnetic field (). In this communication, we shall concentrate on the numerical evaluation of the series. We will remind the reader in Section 2 that the form factor series can be neatly re-summed in terms of a Fredholm determinant different from the one in [59]. The kernel of the integral operator that defines the Fredholm determinant is a strongly oscillating function for large values of and . We shall demonstrate that, nevertheless, the idea presented in [67] 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 . Consequently, we are able to numerically check the asymptotic results of [58, 61]. We are also able to suggest a higher-order correction for 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 organised as follows. In Section 2, we start from a brief review of the results in [50] 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 3. The steepest-descent paths in the space- 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 4. 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 5. In Section 6 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 [59] is given in Section 7. Thanks to the high precision calculation, we identify a higher order correction in the massless phase. In Section 8 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 9. Some of the technical details as well as a comparison with the exact static values are supplemented in several appendices.

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

We consider the spin- XX chain,

 H=JL∑i=1(σxiσxi+1+σyiσyi+1)−h2L∑i=1σzi (1)

with periodic boundaries in the thermodynamic limit . We assume .
There are two critical values

 hc:=4J,tc:=m4J (2)

of magnetic field and time which will occur frequently in the following. The system is in the massless phase if and in the massive one if .
The quantity of interest in this work is the transverse correlation function

 ⟨σ−1(0)σ+m+1(t)⟩J,h,T (3)

for arbitrary and . The subscripts will often be suppressed below. We divide the space-time plane into two regimes: the space-like regime with and the time-like regime with . 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 form factor expansion to this model was already described in detail in Section 3.5 of [50]. We start from a brief summary of the result obtained there.

2.1 A summary of previously obtained results

Eq. (3) is written more explicitly as

 ⟨σ−1(0)σ+m+1(t)⟩=limL→∞tr1,⋯,L(e−(1T+it)Hσ−1eitHσ+m+1)tr1,⋯,L(e−1TH). (4)

We then apply a Suzuki-Trotter decomposition [64] to the exponential factors in a slightly sophisticated manner, respecting the integrability of the original quantum chain [38, 65, 66]. As a result, we introduce an alternating vertex model acting on a fictitious space of size , where is the Trotter number, and the associated QTM [50]. We then substitute the solution of the inverse problem to represent all quantities in (4) by entries of the quantum monodromy matrix with a twist angle . In terms of the physical parameters the twist is given by . We denote the QTM by , its -th eigenvalue by and its eigenstate by . We assume that the numbering is such that is the dominant state, i.e. the unique eigenstate of whose eigenvalue has maximal modulus. In order to deal with the Trotter limit , we follow the general idea of the form factor expansion and represent the correlation function by inserting a complete set of eigenstates between the operators and then summing up this expansion by alternative methods.

Every state is characterized by the associated auxiliary function . In general, is given as solution to a non-linear integral equation. For the XX model, however, due to its non-interacting nature, we know its explicit form given by

 an(λ)=(−1)se−ϵ(λ)T,ϵ(λ)=h−4iJsh(2λ). (5)

Here takes values or , depending on the state, and the dominant state has . The Bethe roots satisfy

 an(λa)+1=0. (6)

All Bethe roots of the dominant state are located in the strip . Some roots of the excited states are distributed in the strip which are referred to as “particles”. There are additional solutions of (6) in which do not belong to the set of Bethe roots and which will be called “holes”. Note that they are defined modulo in the imaginary direction due to the periodicity of the functions . The hole-particle excitations with holes and particles exhaust all possible excitations of the QTM of the XX chain.

The matrix elements between and are readily evaluated by using the quantum inverse scattering method (QISM) algebra and Slavnov’s formula [7]. They are expressible using the above auxiliary functions. The non-vanishing elements are given by excitations with and in (6) for the transversal correlation function.

Following [50] we set

 e(λ)=2sh(2λ), (7) z(λ)=12πiln(1+a0(λ)1+an(λ)), (8) Φ(x)=e(x)2×exp{2∫Cdμcth(x−μ)z(μ)}, (9) D({xj}nhj=1,{yk}npk=1)=[∏1≤j

and

 A=exp{2∫Cdμcth(2μ)z(μ)−∫C′⊂Cdλ∫Cdμcth′(λ−μ)z(λ)z(μ)}, (11a) A(m)=A×exp{−∫Cdμz(μ)me(μ)}. (11b)

The contour tightly encloses in (11a), and it simply surrounds counterclockwise the strip in the massive phase (). In the massless phase () the contours must be slightly deformed due to the existence of two Fermi points. A precise prescription for the deformation will be given below after a suitable change of integration variables.

The mentioned summation over all the excited states of the quantum transfer matrix, or equivalently, all possible Bethe roots configuration is replaced, in the Trotter limit , by a series of contour integrals which takes the form

 (12)

The rapidities are associated with holes and the with particles. The symbol denotes the boundary values of from inside the contour. The particle contour is chosen as .

In order to improve the behavior of integrand in the exponent of at , we slightly modify the definition (9),

 ~Φ(x):= Φ(x)exp{2∫Cdμcth(2μ)z(μ)} =e(x)2×exp{2∫Cdμe(μ)z(μ)sh(x+μ)sh(x−μ)} =Φ(x)th(h2T), (13)

where we used a fact that () is even (odd) on the line . Thus, one can equivalently write

 ⟨σ−1(0)σ+m+1(t)⟩=(−1)m+1th(h2T)A(m)∞∑n=11n!(n−1)!∫Cn∏r=1dxrπi~Φ−(xr)ei(mp(xr)−tϵ(xr))1−eϵ(xr)T×∫Coutn−1∏s=1dysπie−i(mp(ys)−tϵ(ys))~Φ−(ys)[1−e−ϵ(ys)T]D({xr}nr=1,{ys}n−1s=1). (14)

Here we use a contour which encloses by utilizing the periodicity in the imaginary direction. Note that the factor is absorbed by reversing the integration direction. This is the main result in [50] for the transverse correlation function of the XX chain.

2.2 A representation by a Fredholm determinant

We find it useful to adopt the momentum variable for the numerical calculation,

 p(λ)=1ilnthλi. (15)

We choose the branch cut to be modulo . The bare and the dressed energy are represented in familiar forms,

 e(λ)=2cosp(λ)iϵ(λ)=h−4Jcosp(λ).

Hereafter we use the same symbol to denote the dressed energy as a function of the momentum variable .

The substitution in (14) leads to

 ⟨σ−1(0)σ+m+1(t)⟩= (−1)m+1th(h2T)A(m)∞∑n=1(−1)n−1n!(n−1)!∫En∏r=1dprμ(pr)ei(mpr−tϵ(pr)) ×∫¯En−1∏s=1dqs¯μ(qs)e−i(mqs−tϵ(qs))D2p({pr}nr=1,{qs}n−1s=1), (16) μ(p)=eσ+(p)2π(1−eϵ(p)/T),¯μ(q)=e−σ−(q)2π(1−e−ϵ(q)/T), (17) σ(p)=∫Edq2πi1tan(p−q2)ln(1+e−ϵ(q)/T1−e−ϵ(q)/T), (18) Dp({pr}nr=1,{qs}n−1s=1)=∏1≤j

where , resp. is the boundary value as approaches from the above, resp. below. We call hole momenta and particle momenta. The contour is a straight line in the massive phase, . Similarly, is a contour just below the real axis. When , the contours make a little detour in order to avoid the Fermi points 111The Fermi momentum with the proper dimension will be introduced later in Section 6. (defined by ), as shown in Figure1. The contours and are images of and under the transformation (15).

We will now rewrite (16) by a Fredholm determinant. For this purpose, let us prepare some notations in advance,

 ~V(qi,qj) =∫Edpμ(p)φ(p,qi)φ(p,qj)e2tum,t(p),V(qi,qj)=e−tum,t(qi)~V(qi,qj)e−tum,t(qj), (19) ~v(qi) =∫Edpμ(p)φ(p,qi)e2tum,t(p),v(qi)=e−tum,t(qi)~v(qi), (20) Ω(m,t) =∫Edpμ(p)e2tum,t(p), (21) P(qi,qj) =v(qi)v(qj)Ω(m,t), (22) φ(p,q)=eiq−p2sin(q−p2),um,t(p)=i2(mtp−ϵ(p)). (23)

We further define integral operators,

 (^Vf)(q) =∫¯Edq′V(q,q′)f(q′)¯μ(q′),(^Pf)(q)=v(q)Ω(m,t)∫¯Edq′v(q′)f(q′)¯μ(q′).

The operator is obviously a one-dimensional projector.
After these preparations, we have [62],

Proposition 1.

The transverse correlation function of the spin- XX model is represented by a Fredholm determinant using the integral operator ,

 ⟨σ−1(0)σ+m+1(t)⟩=(−1)m+1th(h2T)A(m)Ω(m,t)det¯E(1+^P−^V). (24)

The symbol means the Fredholm determinant where the integral operator is acting on functions supported on .

Although the right hand side of eq. (24) contains the factor , it is canceled by the same factor in . Thus, the limit gives a finite number.

We also remark that belongs to the class of integrable integral operators. This is easily seen from the identity

 (25)

A Fredholm determinant with a similar but different kernel, which nevertheless belongs to the family of integrable integral operators, was analyzed in [58]. The matrix Riemann-Hilbert approach provides a powerful tool for the analysis of such operators. The leading terms of the large space and time scale behavior of the transverse correlation function of the XX model have been successfully derived within this framework.

Here we are interested in a quantitative study on an arbitrary space and time scale. We shall employ a numerical method for this purpose.

3 The analytic properties of integrands and the steepest descent paths

The contours are optimal choices in the static limit. As time evolves, the phase factors 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 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.

3.1 The saddle points and the steepest descent paths

We need to take account of the steepest descent paths for in the asymptotic region . 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.

3.1.1 The space-like regime t<tc

The saddle points in the space-like regime lie at

 p±=π2±iarch(m4Jt). (26)

The red points in Figure 2 denote them. The dashed curve corresponds to the locum of points such that , and the shaded region satisfies . Thus, the steepest descent path for the hole momentum passes through horizontally and never enters into the shaded region. We have an upside-down figure for the path of the particle momentum which passes through , as it has the conjugate phase . As a result, the optimal path for the hole (particle) momentum in the space-like regime is similar to (): it stays above (below) the real axis.

3.1.2 The time-like regime t>tc

In this regime the saddle points are on the real axis,

 p±=π2±arccos(m4Jt). (27)

They are depicted by red dots in the right panel of Figure 2. The shaded region indicates . Thus, for the hole momentum, the steepest descent path runs through and , while staying inside of the un-shaded region. Accordingly, it cannot stay above the real axis, but passes below the real axis for . The figure for the particle momentum is obtained by turning the figure for the hole momentum upside down.

3.2 The analyticity of μ and ¯μ

We only have to know () in (17) as functions of 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 () phase and for the massive () phase.

3.2.1 The massive regime h>hc

We set

 qj=arccos(h+ijπT4J)(j∈Z), (28)

so that

 1+e−ϵ(q2j−1)/T=0,1−e−ϵ(q2j)/T=0. (29)

We further define “upper roots” and “lower” roots by

 quj={qjj≤0,−qjj≥1,qdj={−qjj≤0,qjj≥1,

so that and .

The roots closest to the real axis are and .

The following analytic properties of are derived in Appendix A  .

Lemma 1.

The function has
i) simple poles at in the upper half plane,
ii) double poles at in the lower half plane,
iii) simple zeros at in the lower half plane.
There is a strip including the real axis which is free of zeros and poles (of width ).

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 are obtained by taking the mirror image w.r.t. 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.

3.2.2 The massless regime h<hc

We denote by , the (right) “Bethe roots”, which are identical to the in eq. (28). Note that . The (left) “Bethe roots” are defined by .

In Appendix A   we show that () has the following analytic properties.

Lemma 2.

The function has
i) simple poles at and at () in the upper half plane,
ii) a simple pole at and a simple zero at on the real axis,
iii) double poles at and at () in the lower half plane,
iv) simple zeros at and at ().
The zeros and poles of 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 and become as . Thus, singularities will approach towards the real axis as and the choice of the contours becomes difficult. This poses a technical problem in the massless regime.

4 Numerical study of the massive regime h>hc

In this and in the next sections we present the results of our numerical study of the transverse correlation function based on (24). We follow the proposal in [67] and evaluate Fredholm determinants by the Gauss-Legendre integration with discrete points .

We shall start from the simpler case . 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.

4.1 The space-like regime t<tc

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 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 and were obtained, e.g., in [44] (for the XXZ model in general). We compare them with the static results obtained from (24) in Appendix B  . The results match with reasonable precision. We have to choose, of course, the appropriate number of discretized points to achieve agreement. For example, fixing and we find

 ⟨σ−1(0)σ+1(0)⟩=⎧⎨⎩0.0186123688⋯exact,0.0186094⋯n=64,0.0186123689⋯n=512.

Thus, we can say that with , eq. (24) practically reproduces the exact value. For larger segments, some indirect evidence is presented in Appendix C  . 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 .

Figures 4 and 5 show examples for , and for various . The horizontal axes indicate . Although the space-like regime is limited to when and when , the straight contours produce stable values even after .

The real part stays almost flat in the space-like regime. This can be better seen in the case of in Figure 5. The amplitude is enhanced around .

For larger values of (typically ), we need to adopt contours which take into account the steepest descent paths. The saddle points in eq. (26) are located away from the real axis, and we shift the contour ( ) to a straight line passing though (). The integration contours cross poles of and if , which is clearly seen in Figure 3 (left). This modifies our formula (24) slightly, following the argument in [48], which will be summarized in Appendix D    for the reader’s convenience.

4.2 The time-like regime t>tc

The calculation with the straight integration contours gradually becomes unreliable for larger . 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 () and for the particle momentum may be the ones depicted in Figure 7, which respect the steepest descent paths discussed in Section 3.1.

The paths are determined only by the phases and . One, however, needs to be careful: Since and are swapped, there appears an additional term in ,

 ~v(q)=∫Cpdpμ(p)e2tum,t(p)φ(p,q)+4πie2tum,t(q)μ(q),

due to a pole of at . Here denotes the red path in Figure 7. As we interpret , resp. , as the momentum of a hole, resp. a particle, we refer to the singularity as the kinematic pole in analogy with scattering theory. Consequently, the kernel of the integral operator contains a term

 −16π2Ω(m,t)etum,t(qi)+tum,t(qj)μ(qi)μ(qj).

This becomes very large for if and makes the calculation again unstable. We thus choose the red contour in Figure 7 for the variables while we adopt a straight contour for the variable such that and for .
Still, there is a problem in dealing with the intersection points . We simply exclude them from the set of discretized sampling points in the integrals over the variable.
After the above modification and by increasing the number of sampling points for the Fredholm determinant (, typically), a large scale calculation is possible. As an illustration, the plots for are shown in Figure 8.

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

The amplitude of 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 .

Unlike on temperature, the period of oscillation obviously does depend on the magnetic field, and it grows as approaches (Figure 11). It may be easier to check this by looking at the imaginary part of (cf. right panel).

The oscillatory period of is easily understood in the asymptotic region and . A saddle point analysis yields

 Ω(m,t)∼√π2∣∣u′′m,t(p+)t∣∣μ(p+)ei(mp+−ϵ(p+)t)+√π2∣∣u′′m,t(p−)t∣∣μ(p−)ei(mp−−ϵ(p−)t). (30)

Note that , and . Due to the asymmetry of the integrand in (18) for or , one can show that so that . On the other hand . Thus,

 |μ(p−)|=∣∣∣eσ+(p−)2π(1−eϵ(p−)/T)∣∣∣≫|μ(p+)|

and we can safely drop the first term in (30). Thus, we are left with a single oscillating term and conclude that the period of the oscillation diverges as , for some constant , when . Note that is an exceptional point as the Fermi points pinch the real axis and the above argument is not valid.

5 Numerical study in the massless regime h<hc

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.

5.1 The space-like regime t<tc

We again transform to a straight contour in the upper (lower) half plane for the hole (particle) momentum if 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 4.1 on the massive phase. The details are explained in Appendix D  . We choose a straight contour for . Although the choice of is in theory arbitrary, in practice one should not take it too small. Table 1 illustrates this with the examples of , and for various choice of . The rightmost column gives the exactly known static values.

This can be easily understood as the contour passes near the singularities of (at ) or (at ). Hence, too small values of spoil the numerical accuracy.

With suitable choices of 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 for small temperatures is given by the effective bandwidth in the Fermionic model, , and goes to zero as approaches as demonstrated in Figure 13.

For larger values of , is determined by the location of the saddle point (26) and one has to the take into account the pole contributions as described in Appendix D  . Some numerical results will be supplemented later (cf. Section 8).

5.2 The time-like regime t>tc

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 . The choice of the paths thus becomes a complicated technical issue, and we decided to summarize it in Appendix E  .
By applying the prescription presented there, the result are now stable as shown in Figure 15, where we have choses the same parameters as in Figure 14. The right panel represents the result up to , and the real part of decays to the order of .

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 [57] calculated in the finite XX-chain (400 sites) using the Pfaffian representation. Their result does not converge numerically for various choices of around for (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 8, 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 (24) 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 of discretization points. We defer the discussion of the right choice of to Appendix F  . With suitable values of , we are able to perform a precise evaluation of on a longer time scale. An example is given in Figure 18. The right panel presents a zoom of the curves for large . This exhibits a slowly decaying oscillating pattern.

6 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 for being odd for small . 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 (24). The behavior of is the same for odd and even. We numerically find that in the massless case the oscillation phase of and that of the Fredholm determinant part almost cancel each other and that this results in the monotonous time evolution for odd . By way of contrast the two contributions do not cancel each other for 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 onto spin-less Fermions with annihilation operator . Next, the dispersion is linearized around the Fermi points with where 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 and where 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 is dominant and is given for small temperatures by [36]

 ⟨σ−(0,0)σ+(m,t)⟩(0)∝(2πTv)1/2eiπxsinh1/4(2πTv(x−vt−iδ))sinh1/4(2πTv(x+vt+iδ)). (31)

Here is the sound velocity. Note, in particular, that for and standard bosonization predicts a decay 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 for , see eq. (27) 222 Here the momentum has a proper dimension . 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 where denotes a particle near the saddle point. It is then a fairly straightforward calculation to show that the additional saddle point contribution is [36]

 ⟨σ−(0,0)σ+(m,t)⟩(1)∝ei¯p+xcos(kFx)⟨σ−(0,0)σ+(m,t)⟩(0)⟨d(0,0)d†(x,t)⟩. (32)

We are interested here in understanding the asymptotic behavior in the time-like regime for and . In this limit we have and we can approximate the propagator for the high-energy particle by . Putting together the two contributions (31) and (32) for this case we arrive at

 ⟨σ−(0,0)σ+(m,t)⟩\lx@stackrelvt≫m∝(−1)m(2πT/v)1/2sinh1/2(2πTt)[1+Acos(kFx)ei(4J−h)t√vt] (33)

with some unknown amplitude . In the zero temperature limit, in particular, eq. (33) predicts that the correlation function asymptotically consists of a uniform part which decays as and part which decays as and oscillates with frequency .

Crucially, the oscillating part contains a factor . This means that for () the oscillating contribution is absent if the spatial distance in the two-point correlation function is odd!

To check the predictions of Eq. (33) we present in Figure 21 low-temperature numerical data for the auto- and nearest-neighbor correlation function at (left panel) and (right panel).

From the log-log plot it is clear that the correlation functions decay following a power law in time for all cases. For () the autocorrelation shows oscillations while the nearest-neighbor correlation function does not whereas both are oscillating for