Non-analytic behavior of the Loschmidt echo in XXZ spin chains: exact results

We address the computation of the Loschmidt echo in interacting integrable spin chains after a quantum quench. We focus on the massless regime of the XXZ spin-1/2 chain and present exact results for the dynamical free energy (Loschmidt echo per site) for a special class of integrable initial states. For the first time we are able to observe and describe points of non-analyticities using exact methods, by following the Loschmidt echo up to large real times. The dynamical free energy is computed as the leading eigenvalue of an appropriate Quantum Transfer Matrix, and the non-analyticities arise from the level crossings of this matrix. Our exact results are expressed in terms of"excited-state"thermodynamic Bethe ansatz equations, whose solutions involve non-trivial Riemann surfaces. By evaluating our formulas, we provide explicit numerical results for the quench from the N\'eel state, and we determine the first few non-analytic points.


Introduction
Despite their long history, in the past decade the theory of integrable models has witnessed a series of unexpected developments. Among these, the most prominent one is arguably the realization that analytic techniques from integrability, traditionally tailored for ground-state and thermal physics, provide powerful tools also out of equilibrium [1].
A simple but physically interesting protocol which has proven to be within the reach of integrability has been the one of quantum quenches [2]: a system is prepared in some well-defined state |Ψ 0 and left to evolve unitarily with some Hamiltonian H. This problem has been extensively studied in the past few years, since it represents an ideal and simplified setting for exploring several important questions of many-body physics out of equilibrium [3][4][5][6][7]. Among these, the problem of relaxation has represented a main motivation in the study of quantum quenches [8]: based on the knowledge of the initial state |Ψ 0 , can we predict the local stationary properties of the system at large times?
A complete answer to this question has been obtained: while in the generic case the postquench stationary properties are thermal [8], in integrable systems they are locally captured by a Generalized Gibbs Ensemble (GGE) [9][10][11][12][13][14][15][16]. The latter is analogous to the Gibbs statistical ensemble, but it is constructed by taking into account, in addition to the Hamiltonian, all higher local conserved operators [17]. Besides this established conceptual picture, recent research has also provided us with quantitative means to make explicit predictions in concrete cases; a relevant example is the Quench Action approach [18,19], which allows us to compute, in several cases of interest, the stationary values of local correlations at large times [20][21][22][23][24][25][26][27][28].
Partly motivated by this problem, an analytic computation of the so-called Loschmidt echo in the XXZ Heisenberg chain was initiated in [52,53]. The latter is not a local quantity: it is defined as the squared absolute value of the overlap between evolved and initial states. However, it is of experimental relevance being accessible, for example, by nuclear magnetic resonance [54,55]. Most prominently, the Loschmidt echo is a central object in the study of dynamical phase transitions , and quantum revivals [79][80][81][82][83], and as such, it has received increasing attention over the past few years. Furthermore, the calculation of the Loschmidt echo represents an intermediate step towards the more ambitious goal of computing the time evolution of local observables [52,53].
A promising analytical approach to its computation was proposed in [52], which is based on the so called Quantum Transfer Matrix (QTM) formalism [84][85][86]. This method can be applied quite generally to an infinite family of initial integrable states, which have been introduced and studied in [87]. Building upon the results of [52], a full solution to the problem of computing the Loschmidt echo for imaginary time and arbitrary initial integrable states was presented in [53], while partial results were obtained for real times. In this work we complete the programme initiated in [53] and provide a full solution to the real-time problem, for which a significant amount of additional techniques has to be introduced.
Differently from [52,53], in this work we will focus on quantum quenches to the gapless regime of the XXZ Hamiltonian. The reason to do this is two-fold: on the one hand, some technical simplifications occur, which allow us to reduce the amount of unnecessary complications. On the other hand, for the particular initial states considered, the Loschmidt echo displays, in the gapless regime, non-analytic points at relatively short times [66]; this allows us to show explicitly that our method is perfectly capable to capture them. Since we are only interested in presenting the general methods, we will focus uniquely on quenches from the Néel state, which has already served many times in the recent literature as a prototypical case of study [22,23,30,88]. We emphasize, however, that the method detailed in this work is general, applies for arbitrary values of the anisotropy, and can be carried out for arbitrary integrable states [87].
The organization of this work is as follows. In Sec. 2 we present the XXZ Hamiltonian and the quench protocol, while the QTM approach is reviewed in Sec. 3. We tackle the computation of the Loschmidt echo in Sec. 4, where the small-time dynamics is addressed: in this case, no additional complication arises with respect to imaginary times. The calculation of the Loschmidt echo for arbitrary time is presented in Sec. 5 and Sec. 6, where all the new analytic techniques are introduced. Explicit results for the quench from the Néel state are also reported and discussed. Finally, our conclusions are presented in Sec. 7. The most technical part of our work is consigned to the appendix.

The model
We consider the XXZ spin-1/2 chain where we take J > 0, while σ α j are the Pauli matrices. We assume periodic boundary conditions, σ α L+1 ≡ σ α 1 , and take the length L to be an even integer. We indicate the associated Hilbert space as H = h 1 ⊗ . . . ⊗ h L , where h j C 2 is the local Hilbert space corresponding to the site j. In this work we focus on the gapless regime of the model with γ ∈ R. As a technical hypothesis, we restrict to the special case of anisotropies corresponding to the so-called root of unity points, where γ is a rational multiple of π. We focus in particular on the simplest case where p > 1 is an integer number.

The quench protocol and the Loschmidt echo
In this work we are interested in quantum quenches from a special class of integrable initial states. These have been introduced and defined in [87] to be the states annihilated by all the local conserved operators of the Hamiltonian which are odd under space reflection. They include two-site product states and matrix product states with arbitrary bond dimensions. In order to illustrate the main ideas, we will focus on the simplest example, the well-known Néel state From our derivation, detailed in the following, it will be clear that our approach could be directly applied more generally to all the integrable states defined in [87]. The Loschmidt echo is arguably the simplest quantity to compute after a quantum quench. Given the initial state |Ψ 0 , it is defined as and is a measure of the probability of finding the system close to its initial configuration. For a global quench, L (t) decays exponentially with the volume L, and it is natural to introduce the Loschmidt echo per site or, alternatively, the return rate For global quenches, the analytical computation of the Loschmidt echo is in general extremely hard, and most of the results in the literature are restricted to either free or conformal systems, and local quenches (see however [81] for an analytical computation starting from a domain wall state). An analytical approach was introduced in [52,53] in the gapped regime of the XXZ Hamiltonian, where an exact calculation was presented for the partition function In particular, the starting point of [53] was considering real values of w (imaginary time evolution), for which the exact solution for Z(w) was obtained in terms of analytic formulas. Subsequently, analytic continuation was performed to obtain the real-time evolution of the Loschmidt echo. This procedure was shown to provide the correct result only up to a finite time t * . In this work we go beyond and present a complete derivation of the real-time Loschmidt echo, by considering complex values of the parameter w from the beginning.
In the next section, we review the Quantum Transfer matrix construction introduced in [52,53], which can be carried out straightforwardly also in the massless regime of the Hamiltonian (1).

General idea
The idea behind the Quantum Transfer Matrix approach to the Loschmidt echo relies on interpreting it as a particular boundary partition function; this is a natural identification which has been exploited many times in the literature [57, 61-63, 65, 66, 70-73, 80]. Here, we only briefly review the main formulas for later reference, referring in particular to [53] for a detailed and pedagogical treatment.
We start by introducing the building blocks of our algebraic construction, which is based on the so-called algebraic Bethe ansatz method [89]. The latter is a powerful set of techniques which allows us, among many other things, to analytically diagonalize the Hamiltonian (1). The central object is the R-matrix, from which, one can define the transfer matrix By means of the identity [86] τ where γ is defined in (2), the subsequent action of transfer matrices can be interpreted as a discrete approximation to the unitary evolution, as it follows from the well-known Suzuki-Trotter decomposition Indeed, using (11), one has where we defined As we have already stressed, our approach to the computation of the partition function (8) applies to all integrable states introduced in [87]; these constitute a large family which include matrix product states of arbitrary bond dimension. For the sake of presentation, in the following we will however restrict to two-site product states of the form Following [53], it is straightforward to simplify the partition function (8) by means of the above identities. After a few steps which are not reported here, one obtains where we introduced the boundary quantum transfer matrix and where T QTM (u) is the corresponding monodromy operator In these formulas we have omitted the subscript w in β w . Eq. (16) is the starting point for our derivation, along the lines of [53]. Assuming that the limits of large N and large L can be exchanged, Eq. (16) amounts to compute the w-dependent leading eigenvalue Λ 0 of T ; indeed we have This problem has been completely solved in [53] for arbitrary choices of the state |ψ 0 and w ∈ R. The idea is to relate the operator T to an integrable transfer matrix with open boundary conditions, which can be analyzed by means of the so-called boundary algebraic Bethe ansatz [90][91][92][93]. The form of T explicitly depends on the initial state considered, and for generic integrable states [87] one needs to resort to the non-diagonal version of the latter [94][95][96][97][98][99][100][101], as explicitly worked out in [53]. In this work, to avoid unnecessary complications, we will restrict to initial states which only require us to deal with the simpler diagonal boundary algebraic Bethe ansatz. Furthermore, as repeatedly stressed, we will work in the gapless regime of the Hamiltonian (1), contrary to [53]. In the next section we will thus review the technical aspects of the boundary algebraic Bethe ansatz in the gapless case, referring the reader to the literature for a more comprehensive treatment [90][91][92][93].

The boundary algebraic Bethe ansatz
The central object of the boundary algebraic Bethe ansatz approach is the boundary transfer matrix where Here R ij (u) is the R-matrix introduced in (9). The inhomogeneities ξ j are parameters which can be chosen arbitrarily, and the trace in (20) is performed over the auxiliary space h 0 C 2 .
The 2 × 2 boundary matrices K ± (u) have to be chosen to satisfy appropriate non-linear relations known as reflection equations [90]. In the diagonal case of interest in this work, the general solution to the latter reads To simplify the discussion, we restrict from the beginning to quantum quenches from the Néel state (4). On the one hand, it has already served many times in the recent literature as a prototypical case of study [22,23,30,88]; on the other hand, it is straightforward to apply the techniques introduced in the following to treat more general integrable states, so that this is by no means a restrictive choice. Specifying the initial state to be the Néel state amounts to choosing |ψ 0 = | ↑↓ . Then, following [53], one can show that in this case one has provided that the inhomogeneities and boundary parameters in τ B (u) are chosen as and Eq. (25) is a key relation, which allows us to resort to integrability based methods to compute the Loschmidt echo. Indeed, the problem of computing the partition function (8) is reduced to finding the leading eigenvalue of the boundary transfer matrix (20), where the dependence on w is encoded in the inhomogeneities (26) and (27). The diagonalization of the boundary transfer matrix (20) can be performed analytically within the framework of the boundary algebraic Bethe ansatz [90][91][92][93]. Note that for diagonal boundaries, the transfer matrix (20) commutes with the operator counting the number R of down spins. Then, the eigenstates of the open transfer matrix (20) are constructed in terms of a set of complex numbers {λ j } R j=1 , which are the so-called Bethe roots or rapidities. They are obtained as the solution of a non-linear set of equations (Bethe equations) [91], which in our case read Each set of rapidities λ ≡ {λ j } R j=1 associated with the different eigenstates uniquely specifies the corresponding eigenvalue τ λ (u) of the boundary transfer matrix τ B (u) in (20). First, given so that the corresponding eigenvalue reads [91] where we have defined Eq. (31) is sometimes referred to as the T − Q relation [93]. We note that the Bethe equations (29) can be rewritten in terms of the functions introduced above as Eq. (31) provides a formal solution to the problem of diagonalizing the transfer matrix (10) for finite N . From Eqs. (19) and (25) we see that our goal is to compute the leading eigenvalue in the limit N → ∞. Two main difficulties, accordingly, arise: the first one consists in the determination of the Bethe roots λ ≡ {λ j } R j=1 corresponding to the leading eigenvalue at finite N ; the second pertains the computation of the limit N → ∞ of the expression (31).
The configuration of Bethe roots depends on w [defined in (19)]. For each "time" w, eigenvalues which are close to each other might correspond to very different sets of rapidities; as w varies each set of Bethe roots also varies continuously. However, it might happen that two eigenvalues undergo a crossing: accordingly, the set of Bethe roots corresponding to the leading eigenvalue might change abruptly as w varies smoothly, which makes the computation of the Loschmidt echo non-trivial. It turns out that for w ∈ R no crossing occurs, and the Bethe roots associated with the leading eigenvalue have a similar qualitative behavior for all values of w ∈ R [52]. This is not the case for imaginary times w = it (t ∈ R), considered in this work. In order to set the stage, we will start in Sec. 4 by treating the case of small times t, where no crossing arises. This allows us to focus on a single eigenvalue, for which the configuration of Bethe roots is relatively simple. In Sec. 5 the same eigenvalue is computed for arbitrary values of t, for which the technical treatment becomes necessarily more sophisticated. Finally, for large times crossings arise, as it will be discussed in Sec. 6: in this case, our strategy will consist of computing, for each t, also higher eigenvalues of τ B (u) and to follow their evolution continuously, keeping track of all the subsequent crossings. i=1 associated with the leading eigenvalue of the boundary QTM in the complex-λ plane, for 2N = 8, γ = π/3 and w = it = 1i. The xand yaxes correspond to real and imaginary parts of the rapidities.

The Loschmidt echo at small times
We set w = it in Eq. (8) with t ∈ R, and t sufficiently small. Following [53], we start with a preliminary numerical analysis at finite values of N of the eigenvalues of the boundary QTM (20), which can be obtained by exact diagonalization. It is found that for small values of t the leading eigenvalue of the boundary QTM is unique, with a finite gap with respect to the higher ones. Furthermore, it lies in the sector of zero magnetization, and therefore is associated with R = N Bethe roots. For small values of N , these can be identified numerically following a standard procedure, by comparing the formal T −Q relation (31) with explicit diagonalization of τ B (u), as already done in [53]. An example of a configuration of Bethe roots for the leading eigenvalue of τ B (u) is displayed in Fig. 1.
The Bethe roots do not arrange, in the limit N → ∞, according to a smooth rapidity distribution function. In order to take the infinite-N limit, then, two routes can be followed. The first one consists in writing down a single nonlinear integral equation for an appropriate auxiliary function, as done in [52]. While this can be easily done for the Néel state, serious complications arise in the study of more general integrable states, corresponding to nondiagonal boundary conditions [53]. Accordingly, in order to keep the discussion as general as possible, we will follow the second approach, detailed in [53], which is based on the socalled Y -system relations [102], and for which no additional complication arises in the case of generic states. In the following subsections we will introduce the Y -system and review how the latter can be exploited to obtain directly the spectrum of the boundary QTM.

The Y -system
It is an established result that the boundary transfer matrix τ B (u) can be used to build an infinite family of transfer matrices {t j (u)} ∞ j=0 via the so-called fusion procedure [103][104][105], in complete analogy with the well-known case of periodic boundary conditions [102]. In the following, we outline the aspects of this construction which are relevant for our work.
The fused transfer matrices t j (u) act on the same space as the transfer matrix τ B (u) and form a commuting set of operators, namely They can be obtained recursively as Here f (u) is defined as where the functions φ(u), ω 1 (u) and ω 2 (u), are defined in (33), (34) and (35) respectively. Let us point out that the results presented in this paragraph hold for generic values of the inhomogeneities ξ i and of the boundary spectral parameters ξ ± (note that they can also be extended to the case of non-diagonal boundary reflection matrices [53]). Finally, from Eqs. (38) one can derive the relation where The set of relations (40) is usually referred to as T -system. From the latter one can derive a new set of functional relations, the so-called Y -system, which is expressed in terms of the operators with the choice Y 0 ≡ 0. From this definition, and the T -system (40), the following Y -system is readily derived

Truncation of the Y − system at the root of unity
In general, the Y -system (43) consists of an infinite number of functional relations. This is not an issue, as for practical purposes of numerical evaluation of the Loschmidt echo it can be truncated to a finite number n MAX of them, introducing an error which decreases rapidly as n MAX is increased [53]. There exists, however, a particular case where an exact truncation takes place, and the infinite system is exactly equivalent to a finite one: namely when the parameter q = e iγ is a root of unity. In this work we will restrict to this case, in order to reduce the number of unnecessary complications. As an additional simplification, we will impose another restriction to the values of γ, which we will choose to be of the form (3), with p > 1 integer. This makes the final form of the Y -system particularly simple. Generalization to the case γ = qπ/p, with q, p > 1 integers is possible, but will not be discussed here.
For the values of γ in Eq. (3), an exact truncation of the Y -system takes place due to an additional relation between the fused transfer matrices t p+1 and t p−1 , which can be traced back to the representation theory of the underlying quantum group U q (sl 2 ) with q = e iγ . Such a relation was originally observed for the periodic chain in [106] (see also [107,108]), and for general integrable open boundaries in [94,109]. Recasting the results of the latter in our notations, for the particular case of diagonal boundary reflection matrices, we find where we have defined Noticing now that making use of (44) and of the definitions (42), the following truncated Y −system can be obtained where The Y -system can be further simplified once we impose the boundary parameters to take the values of interest in the present problem, namely ξ ± = ∓γ/2, cf. Eq. (28). Indeed, in this case one has α(u) = 0, and consequently we obtain We note that this form of the Y -system, which is particularly convenient from the computational point of view, holds whenever ξ + = ξ − or ξ + = −ξ − .

From the Y -system to the Loschmidt echo
The importance of the above construction for the computation of the Loschmidt echo is that it allows us to express the latter in terms of the solution of a system of non-linear equations; in turn, these can be easily evaluated numerically to yield the exact value in the infinite-N limit. The same idea was already exploited in [53], and is reviewed in the following. First, from Eq. (37), we see that the operators Y j commute with one another, with the transfer matrices, and by construction with the global magnetization S z (since we are restricting to diagonal boundary conditions). Accordingly, the set of functional relations (43) can be understood at the level of individual eigenvalues of the operators Y j (u). In the following, we indicate as y j (λ) the eigenvalue of Y j (iλ) (note that a rotation of π/2 in the complex plane of the argument of y j (λ) has been performed for convenience); we will refer to y j (λ) as the Y -functions. They satisfy the Y -system where κ(λ) denotes the eigenvalue of the operator K(iλ). It will be useful to know the asymptotic behavior of the Y -functions on Bethe states as λ → ±∞. First, we make use of the simple relation for the magnetization |S z | of a given eigenstate in terms of the number R of the corresponding Bethe roots, Then, from the T − Q relation (31), we can deduce Finally, the asymptotic behavior of the higher Y -functions y j , j ≥ 2, is easily obtained by recursion using (54).
Following [53], we define the normalized boundary transfer matrix where so that T (0) coincides with the operator T in Eq. (25). For each time t, we indicate with {Λ t (λ)} ∞ =0 the set of eigenvalues of the corresponding operator T (λ) [the dependence on t is through the parameter β, defined in (14)]. With each eigenvalue Λ t (λ), is associated a set of Y -functions {y ( ) j (λ)} p j=1 , and the following relations can be derived [53] 1 + y ( ) where Eq. (64) is a functional relation between the eigenvalue Λ (λ) and the Y -function y ( ) 1 (λ). In order to compute the former, and thus the Loschmidt echo (19), we need a final step, namely to cast the functional relations (57)- (59) and (64) into a set of integral equations. This programme will be followed explicitly in the next sections.

Short-time dynamics
The procedure to cast the functional equations into integral ones is well-known in the literature (see e.g. [110]). We summarize it here, for convenience, referring to the literature, and in particular to our previous work [53], for a detailed explanation.
In order to write down the integral equations, one needs to take first the logarithmic derivative of both sides of Eqs. (57)-(59). Next, one performs a Fourier transform of both sides, obtaining a number of integrals along segments with non-zero imaginary parts. Finally, one moves the integration contours back to the real axis; in order to do so, one needs to take into account all the singularities of the logarithmic derivatives, which correspond to the zeros and poles of the functions y j (λ) and κ(λ) inside the so called physical strip. This is the region of the complex-λ plane defined by After performing these steps, one is left with a set of equations of the form where the . . . denotes additional contributions coming from the poles and zeros of the functions y j (λ) and κ(λ). Here, the following notation for the Fourier transform has been so that its inverse reads Eq. (67) can be transformed back to real space, yielding the desired set of non-linear integral equations.
As it is clear from the above discussion, the only piece of information needed in this calculation is the location of poles and zeros of the functions y j (λ) and κ(λ). It turns out that this can be determined analytically for small times, where no additional difficulty arises with respect to the case of imaginary-time evolution. Our analysis of the analytic structure is based on numerical inspection at finite Trotter numbers N of the functions y j (λ) and κ(λ). These can always be obtained implementing the operators Y j (λ) and K(λ) for finite N . Numerical inspection for Trotter numbers up to 2N = 8 reveals the following analytic structure, which is found to be always present for small times t: • y 1 displays the following structure y 1 has a zero of order 2 at λ = 0; y 1 has poles at λ = ±i γ 2 , of order 2 for p even, of order 1 for p odd; y 1 has poles of order 2N at λ = ±i γ 2 + β 2N ; y 1 has zeros of order 2N at λ = ±i γ 2 − β 2N ; • for j ≥ 1, the only poles or zeros in the physical strip (except possible pairs at ±iγ/2) of y j , j ≥ 2, are a double zero at λ = 0 (resp. double pole) for j odd (resp. j even); • the only poles or zeros of κ in the physical strip (except possible pairs at ±iγ/2) are a double pole at λ = 0 for p even.
Note that additional pairs of zeros or poles of the auxiliary functions y j (λ) and κ(λ) at ±iγ/2 do not give contributions to the integral equations and can be neglected. Using the above information and the procedure outlined above, we obtain easily the integral equations corresponding to the functional relations (57)- (59). For p odd, we have ln y j = s * ln(1 + y j−1 ) + s * ln(1 + y j+1 ) + (−1) j 2 ln coth πλ 2γ ln κ = s * ln(1 + y p−1 ) + 2 ln coth πλ 2γ .
Similarly, from Eq. (64), one gets the following relation between Λ and y 1 ln Λ = s * ln where The equations above are exact at finite Trotter number N . It is straightforward to compute the Trotter limit using where we set for convenience J = 1 in (14). The resulting equations can be solved numerically by iteration, and their validity holds until the analytical structure of the Yfunctions remains as outlined above. Note that they are the same that one would obtain by analytic continuation of the imaginary-time result, namely for w ∈ R. Indeed, it was already observed in [53] that the correct real-time Loschmidt echo could be derived in this way for small times. It was already observed in [53], however, that these equations hold only up to a given time 0 < t * < ∞, after which they do not provide anymore the correct prediction for the Loschmidt echo. In the following, we show explicitly that this is due to the fact that at t = t * additional zeros of the Y -functions enter the physical strip, and new source terms of the integral equations have to be considered. This allows us to go beyond the results of [53], and compute the Loschmidt echo for intermediate and large times.

Full time dependence of transfer matrix eigenvalues
The Y -system encoded in Eqs. (57)-(59) is valid at any time t. However, as we already stressed, the integral equations derived in the last section hold only up to a critical value  Fig. 2.
Importantly, we found that the position of the additional zeros and poles for t > t * can not be determined analytically. In order to overcome this issue, we employ a procedure which was initially introduced within the framework of the so-called excited-state thermodynamic Bethe ansatz (TBA). The kinds of techniques that we will employ were first introduced in the context of thermal physics in one-dimensional solvable models [111,112] and integrable quantum field theories [113,114], and will be illustrated in the following.
For simplicity, we consider the case p = 2, for which the Y -system reads From numerical inspection, we see that no additional poles enter the physical strip, and only zeros of the Y -functions appear, which always come in pairs of opposite value. Suppose that additional zeros of κ(λ) enter the physical strip for a given time t. The contributions of zeros and poles are clearly additive, so we can consider a single pair of zeros ±δ (κ) . Note that we label arbitrarily one of them δ (κ) and the other −δ (κ) . Define in the following Up to a global constant, applying the usual trick of integration in the complex plane we get the following term in the r.h.s. of (67) This could be integrated to give where Note that (86) is symmetric under δ (κ) → −δ (κ) , as it should. The calculations for additional zeros of y 1 (λ) are exactly the same. One should only pay attention to the fact that now zeros of y 1 are of order 2: if this was not the case, the function 1 + κ(u) would display a point of non-analyticity. Again, up to a global additive constant, we get the additional term We can collect these calculations and provide the final result for the integral equations in the presence of additional zeros. Suppose that y 1 (λ) and κ(λ) have respectively n y and n κ additional zeros in the physical strip; then, in the Trotter limit N → ∞, we obtain the following set of TBA equations ln κ(λ) = 2 ln coth 3λ 2 + while C 1 and C 2 are two constants which should be fixed for the particular eigenvalue investigated. Indeed, noticing that lim λ→∞ G(λ, δ) = 0, and defining we obtain from Eq. (61) In the same way, the equation for the eigenvalue of the transfer matrix has to be modified in the presence of additional zeros as Since the values of the zeros {δ (y/κ) j } are not known analytically, they need to be determined self-consistently. In particular, using the Y -system relations, they are immediately seen to satisfy κ δ (y) These equations complement those in (89) and (90), and finally allow us to compute the realtime evolution of a given eigenvalue Λ t (λ). In order to obtain explicit numerical results, one can proceed as follows. First, one starts with an initial guess on the position of the additional zeros and poles {δ (y/κ) j }. Using this guess, one solves the integral equations (89) and (90), yielding an approximation for y 1 (λ) and κ(λ). Next, employing the latter, one solves Eqs. (97) and (98) for {δ (y/κ) j }, which serve as an improved guess for the next iteration. The additional zeros follow non-trivial trajectories in the physical strip, as displayed in Fig. 3. Furthermore, their number can also vary in time. In Fig. 4 we display the leading eigenvalue Λ t 0 for different values of the anisotropy as a function of time. In each plot, we also specify the number of additional singularities which enter the physical strip, and which have to be consistently determined from Eqs. (97), Eqs. (98).
From the numerical point of view, there is an additional non-trivial complication. Indeed, the driving term in Eq. (89) is imaginary and one needs to be careful with the determination of  Figure 4. Time-evolution of the logarithm of the eigenvalue Λ t 0 (0), which is the leading one at short times. The two plots correspond to different values of the anisotropy. In the figures, we explicitly indicated the number of additional zeros of κ(λ) entering the physical strip for each time interval (no additional zeros of y 1 (λ) are seen to appear). the branch of the logarithm. In fact, in order to obtain a continuous solution to these equations, one can not avoid to consider the logarithm as a function defined on a multi-sheeted Riemann surface. Accordingly, y j (λ) and κ(λ) need to be thought of as functions taking value in this surface. For completeness, a detailed discussion on this issue is reported in Appendix A, together with other technical aspects of the numerical solution to the non-linear integral equations.

The full spectrum of the Quantum Transfer Matrix
In the last sections, we have solved the problem of computing a single eigenvalue of the boundary transfer matrix for real times. In particular, we have followed the evolution of the eigenvalue Λ t 0 which at t = 0 is the leading one. As we have already mentioned, however, a crossing of eigenvalues will in general occur after a certain timet: for t >t the eigenvalue Λ t 0 will not be the leading one anymore. As it should be clear from our discussion in the previous sections, the Bethe ansatz method allows us to follow the dynamic of a single eigenvalue continuously, starting from a given time t. Ideally, then, one should compute for a given time the full spectrum of the transfer matrix, so that one could keep track of each crossings of the eigenvalues at later times. One can summarize the procedure to do so, as follows: • diagonalize the transfer matrix at finite Trotter number; • for each excitation, find the location of the additional zeros of the functions y j and κ; • use these as an input for the "excited-state" TBA procedure described in the previous section.
Let us follow these steps in detail for the first few leading states at time t = 0.7. While the procedure works in principle for states with arbitrary values of the magnetization S z , the leading QTM eigenvalue always appears to lie in the sector S z = 0 and we will therefore   (66)] associated with the first leading eigenvalues of the boundary QTM in the zero magnetization sector, at time t = 0.7 and for p = 2. The multiplicity of zeros, when different from 1, is indicated by brackets. The last column is obtained from the self-consistent solution to Eqs. (97) and (98).
restrict to the latter in what follows. The time t = 0.7 lies prior to any crossing, and the leading eigenstate is that studied in the previous section. The next leading states are characterized by a set of additional zeros in the physical strip, which we sum up in Table 1.
The location of these additional zeros is quite stable upon increasing N , and can be reliably used as an input for the iterative scheme described in Sec. 5. The resulting eigenvalues are plotted in Fig. 5. Importantly, these allow us to observe a first crossing between the levels 1 and 2 at time t 3.05, yielding a first point of non-analyticity of the Loschmidt echo.
In order to observe crossings at later times, one could in principle follow the same approach for further excited states. However, it turns out that this is not convenient, as the states involved in crossings at later times are found to lie rather deep in the spectrum of the relevant boundary transfer matrix for t = 0.7, and are therefore difficult to identify systematically. Accordingly, we proceed in a more pragmatic way. In particular, we study the boundary QTM spectrum as a function of t for finite sizes 2N = 6, 8, 10, and identify the states such that the associated eigenvalues become the leading one within a finite given time window. In this way we managed to identify the states involved in two subsequent crossings, for which we characterized the additional zeros at a time t = 4. Next, we solve the resulting integral equations to arbitrary times. The final result of this procedure is shown in Fig. 5.
By selecting at each time the leading eigenvalue, one is left with the final exact result for the return rate, and hence the Loschmidt echo per site. This is shown in Fig. 6, for different values of the anisotropy. Our results were tested against iTEBD simulations [115], and calculations from exact diagonalization at finite system size, displaying perfect agreement. As time is increased further, several additional points of non-analyticities are expected to arise; these should correspond to eigenvalues lying deeper and deeper in the spectrum at smaller time. This is in fact a limitation of our method, as these states become increasingly difficult to track in the spectrum of the QTM at finite N . In order to tackle arbitrary time and, in particular, the problem of the asymptotics, it would be much more satisfactory to have at hand a set of integral equations incorporating in a self-consistent way the analytical properties of the leading eigenvalue throughout its crossings. While we have not been able to achieve this goal at present, we hope to return to these issues in future works. Before closing this section, we point out a rather remarkable feature that we have observed, namely that all of the crossings involving the leading eigenvalue seem to coincide exactly with a change in the analytic structure of the y and κ functions. For instance, precisely at the location of the first crossing between the levels 1 and 2, the number of additional pairs of zeros of κ for the level 1 changes from zero to 1. We were not able to provide a theoretical justification for this phenomenon, and at this stage we report it as a simple observation.

Conclusion
We addressed the computation of the Loschmidt echo in the XXZ spin-1/2 chain, for a special class of integrable initial states [87]. By employing a QTM approach, we have provided an analytic solution at real times, completing the programme initiated in [52,53]. Within our method, the Loschmidt echo is obtained from the leading eigenvalue of an appropriate boundary QTM. As the time increases, crossings occur giving rise to points of non-analyticity, which are fully captured by our techniques.
Although our approach could be in principle used to study the full spectrum of the boundary QTM, one is practically limited in the number of eigenvalues which can be computed. This results in a limitation in the time interval which can be considered: indeed, as  Figure 6. Exact return rate as a function of time, as computed using the techniques detailed in Sec. 5 and 6. The plots correspond to anisotropies ∆ = 1/2 and ∆ = 1/ √ 2. The points of non-analyticities are explicitly highlighted with arrows, and correspond to crossings of the eigenvalues of the boundary transfer matrix, cf. Fig. 5. Different colors correspond to the fact that the return rate is determined by different eigenvalues which become the leading ones at different times. time increases, more and more crossings are expected, and a very large number of eigenvalues should be taken into account. In particular, using our method, we do not have access to the regime of infinite times, and the study of the asymptotic behavior of the Loschmidt echo remains an interesting open question to be investigated.
Our calculations show that TBA techniques, traditionally tailored for thermal physics, can be successfully used to obtain explicit analytic results also for real-time dynamics. Hence, our approach constitutes a promising direction towards the ambitious goal of computing the time evolution of local observables after a quantum quench. Applications of the techniques employed in this work to this very important question are currently under investigation.

Acknowledgments
We are very grateful to Mario Collura for providing us with iTEBD data, to test our numerical scheme for evaluation of the analytic formulas derived in this work. We acknowledge useful discussions with Pasquale Calabrese, Patrick Dorey, Fabian Essler, Gábor Tákacs and Roberto Tateo. B.P. acknowledges support from the "Premium" Postdoctoral Program of the Hungarian Academy of Sciences, the K2016 grant no. 119204 and the KH-17 grant no. 125567 of the research agency NKFIH. E.V. acknowledges support by the EPSRC under grant EP/N01930X.

Appendix A. The Riemann sheet TBA
We wish to illustrate the numerical solution to the TBA equation (89) and (90) for t ∈ R. In order to do this, we consider the explicit case of the leading eigenvalue of the transfer matrix.
It is convenient to employ the parametrization Furthermore, for reasons that will be clear later, it is convenient to introduce also the following parametrization All the functions ρ y (λ), R y (λ), ϕ y (λ), Φ y (λ), are real for real values of λ. The integral equations above can be rewritten as ln ρ κ = s * ln R y + ln coth 2 3λ 2 , (A.9) ϕ κ = s * Φ y . (A.10) The system above consists of 4 equations with 8 unknown functions, so it can not be solved unless additional constraints are imposed. Indeed, the functions R κ/y (λ) and Φ κ/y (λ) are not independent from ρ κ/y (λ) and φ κ/y (λ). However, we argue that the dependence is non-trivial and in general "non-local": in order to obtain Φ κ/y for a given λ it is not enough to know the value taken by functions ρ κ/y and φ κ/y at the same λ, but additional, non-local information should be provided.
To understand this better, we consider a very simple example. Suppose we assign the functions As λ varies, one can follow the evolution of 1+y 1 (λ), and compute R y (λ) and Φ y (λ) [defined in (A.6)] accordingly. Note that if we do this naively, for example choosing Φ y (λ) = log (1 + y 1 (λ)) , (A.12) with a fixed branch cut for the logarithm (for example (−∞, 0]), we obtain a discontinuous function Φ y (λ): a discontinuity arises every time 1 + y 1 (λ) crosses the branch cut (−∞, 0]. On the other hand, we assume that the correct solutions to (A.7)-(A.10) are regular functions, so that no discontinuity should arise. In order to obtain continuous solutions R y (λ), Φ y (λ), we introduce an infinitely sheeted Riemann surface, with a branch cut on the semi-infinite line (−∞, 0]. Then, we follow the  Figure A1. Evolution (as a function of λ) of different curves corresponding to (A.11). In subfigure (a) we plot the projection of the curve 1 + y 1 (λ) on a single sheet of the Riemann surface, as λ is increased from λ = 0 to λ = 2. In subfigures (b) and (c) we report the functions Φ y (λ), R y (λ) which are computed by following the evolution of 1 + y 1 (λ) as λ increases.  Figure A2. Evolution (as a function of λ) of different curves corresponding to (A.13). In subfigure (a) we plot the projection of the curve 1 + y 1 (λ) on a single sheet of the Riemann surface, as λ is increased from λ = 0 to λ = 2. In subfigures (b) and (c) we report the functions Φ y (λ), R y (λ) which are computed by following the evolution of 1 + y 1 (λ) as λ increases.
curve 1 + y 1 (λ) on such a Riemann surface, and compute R y (λ) and Φ y (λ) accordingly. This strategy is displayed in Fig. A1. In subfigures (a) we plot the projection of the curve 1 + y 1 (λ) on a single sheet of the Riemann surface, as λ is increased from λ = 0 to λ = 2. By following the curve, one can compute continuous functions Φ y (λ) and R y (λ) which satisfy (A.6). In subfigures (b) and (c) we report the functions Φ y (λ), R y (λ) computed in this was. To see that the value of Φ y (λ) does not uniquely depend on the value of ρ y and ϕ y computed in λ, consider a different curvẽ ϕ y (λ) = 2πλ ,ρ y (λ) = 1 + 1 2 λ . (A.13) Once again, as λ varies we obtain a curve in the infinitely sheeted Riemann surface corresponding to 1 + y 1 (λ). Its projection on a single Riemann sheet is reported in subfigure (a) of Fig. A2, while the corresponding continuous functions Φ y (λ) and ρ y (λ) are displayed in subfigures (b) and (c). Importantly, we see that ϕ y (λ = 2) =φ y (λ = 2) = 4π , (A.14) ρ y (λ = 2) =ρ y (λ = 2) = 2 . From this discussion it follows that, in order to find regular solutions of the system (A.7)-(A.10), one needs to take explicitly into account an infinitely-sheeted Riemann surface. In practice, we use the following numerical scheme: (i) Initialize the values of R κ/y and Φ κ,y .
(iv) Repeat steps (ii) and (iii) until convergence is reached. 1 + y 1 (λ) lies entirely in one single Riemann sheet, and a solution to Eqs. (A.7)-(A.10) is straightforward. As time increases, the same curve eventually enters a new Riemann sheet. For example, for t = 2 [subfigure (b)] the curve crosses the branch cut (−∞, 0] once, and enters into a new Riemann sheet. We see that, as the curve 1 + y 1 (λ) in general wraps around the origin several times, it is not possible to choose a branch cut of the logarithm such that the latter lives on a single Riemann sheet, and necessarily a multi-sheeted surface has to be introduced. This could also be seen from Fig. A4, where we plot the functions Φ y (λ) and R y (λ) at time t = 5.
It is important to note that at each time such that the curve 1+y 1 (λ) enters a new Riemann sheet, a new zero of κ appears in the physical strip, and hence an additional driving term has to be included in the equations. For example, in the particular case of the plot reported in subfigure (b) of Fig. A3 the integral equations become ln ρ y = 2s * ln R κ − ln coth 2 3λ 2 + Re G λ, δ (y) , (A.18) ϕ y = − 2π sin(γ)t s(λ) + 2s * Φ κ + Im G λ, δ (y) , (A. 19) ln ρ κ = s * ln R y + ln coth 2 3λ 2 , (A.20) where G(λ, δ (k) ) is given in (91), while δ (k) has to be computed self-consistently by solving Eq. (98). Note that the numerical iteration to compute the position of the additional zeros is rather stable and efficient, especially when an initial reasonable guess is given. In fact, it only becomes delicate in correspondence of those times for which a new zero enters the physical strip and a new contribution has to be taken into account; in these cases, care must be taken to follow continuously the solution as time is increased slowly, and to provide an accurate initial guess for the position of the additional zeros.