The propagator of the finite XXZ spin-$\tfrac{1}{2}$ chain

We derive contour integral formulas for the real space propagator of the spin-$\tfrac12$ XXZ chain. The exact results are valid in any finite volume with periodic boundary conditions, and for any value of the anisotropy parameter. The integrals are on fixed contours, that are independent of the Bethe Ansatz solution of the model and the string hypothesis. The propagator is obtained through a lattice path integral, which is evaluated exactly utilizing the so-called $F$-basis in the mirror (or quantum) channel. The final expression is similar to the Yudson representation of the infinite volume propagator, with the volume entering as a parameter. The contour integrals involve an amplitude describing the particular propagation process; this amplitude is similar to (but not identical with) the Bethe Ansatz wave function. An important feature is that it depends on two sets of coordinates (initial and final), and it is manifestly periodic for an arbitrary set of rapidities.


Introduction
Quantum integrable models are special theories that can describe strongly correlated many body systems, such that their Hamiltonians can be diagonalized using exact methods. Their study goes back to the solution of the Heisenberg spin chain by Bethe [1]. Whereas the largest part of the literature is devoted to the study of the state functions and correlation functions in the ground state or at finite temperatures, recently considerable interest was also devoted to the study of out-of-equilibrium situations. This is motivated by experimental advances that make it possible to measure the dynamical properties of isolated quantum systems, thus prompting for a theoretical understanding of the observed phenomena. Two particular areas that have been actively investigated in recent years are the equilibration/thermalization of isolated integrable systems, and the description of their transport properties.
Regarding equilibration the main paradigm is that of the Generalized Gibbs Ensemble (GGE) [2]: it is now widely believed that in homogeneous situations isolated time evolution leads to a steady state that is characterized by a complete set of (local and quasi-local) conserved charges of the model. The imprecise notion of the "complete set" of the charges can be defined rigorously by focusing on particular models; for example in the spin-1/2 XXZ chain the complete GGE was established using Bethe Ansatz techniques in [3]. Furthermore, in this model there are exact methods to compute the mean values of local observables in the highly excited states that emerge after the quantum quenches [4,5]. On the other hand, the GGE has not yet been established for models with higher rank symmetries.
The second area, namely the description of transport in integrable systems has been treated by the so-called Generalized Hydrodynamics (GHD) [6,7,8]. In the first approximation this theory captures the physics at the Euler-scale (a combined large time and long distance limit), and it rests on the dissipationless scattering in integrable models, guaranteeing ballistic propagation of the quasi-particles. However, it has been shown recently that diffusion can also be described within the same framework [9]. The GHD has been established for a number of models and it is extremely successful: it provides predictions that agree with DMRG calculations up to many digits [7,10,9]. The theory can also describe two-point or higher point correlation functions in the Euler scale limit [8].
It would be interesting to go beyond the GHD and derive exact results for time evolution in interacting integrable models. One the one hand, this could provide a microscopic derivation of the GHD. On the other hand, it could give access to physical effects beyond its reach. In the following we review some of the approaches towards exact treatment of non-equilibrium dynamics, with a special focus on the XXZ chain.
First of all, the direct method consists of the insertion of (one or two) complete sets of eigenstates in finite volume and of the computation of the resulting spectral series. Generally, such a task is quite hopeless, at least from an analytic point of view. Nevertheless, depending on the problem it can be treated numerically [11], or it can yield analytic expressions for power-law correction terms in the long-time limit [12,13] using the Quench Action logic [14]. A necessary ingredient in any such calculation is to have exact formulas for the overlaps with the initial states; these are known in a number of cases [15,16,17,18,19,20]. It was also argued recently in [21] that one should focus on a sub-class of initial states (called integrable states) where factorized overlaps can be expected. We should note that for such states the time evolution of the von Neumann and Rényi entropies have been computed in [22,23,24,25].
An independent approach is that of the Quantum Transfer Matrix (QTM) method, which was originally devised to compute the thermodynamics of the spin chain [26]. The main idea is to build a lattice path integral for the partition function, which can be evaluated in the socalled quantum (or rotated, or mirror) channel by exchanging the space and time coordinates. This way the summation over all the eigenstates of the system is replaced by the focus on a single leading eigenstate of the QTM. The method was generalized in [27,28,29] to yield the Loschmidt amplitude in certain homogeneous quenches. It was already argued in [27] that even time dependent local correlators could be computed with the QTM, somewhat analogous to the determination of finite temperature correlators (see [30] and references therein). However, the computations have not yet been carried out and are expected to be considerably more involved.
A further idea towards exact treatment of real time dynamics is through the Yudson representation for the propagator of integrable models. Originally developed in [31,32] and worked out for the Lieb-Liniger and XXZ models in [33,34,35,36], this method computes the propagator of a finite number of particles in an infinite volume system. It is built on two basic ideas. First, it uses the fact that the Bethe wave functions form a complete set in infinite volume, and instead of a summation over the Bethe roots (solutions to the Bethe equations in finite volume) one needs to integrate over the rapidities with appropriate weight functions [37,38]. Depending on the model one can have remarkable simplifications, for example in the XXZ chain or the Gaudin-Yang model the Yudson representation involves single integrals over certain contours instead of a sum of integrals over all string states [36,39]. The second idea is more technical: in the resulting multiple integral formula one can replace one side of the propagator (corresponding either to the "in" or the "out" configuration) by a free wave function, leading to a further considerable simplification. We should also note that the basic ideas of the Yudson representation appeared independently in [40].
It was already demonstrated in [34,35,36] that the Yudson method can yield concrete predictions for real time evolution of observables, nevertheless it has severe limitations. First of all, the number of integrals in the propagator is always equal to the number of particles involved. Thus it is quite difficult to take the physical thermodynamic limit. Second, it is an infinite volume method, therefore it can only describe physical processes where the particles disperse into infinity after some initial interaction.
In the present work we compute the propagator of the XXZ chain in finite volume. The advantage of our approach is that the finite volume propagator can describe both spatially homogeneous and inhomogeneous situations, and it could also be used to study finite volume effects. Our derivations are independent from the works on the Yudson method: we use the QTM approach to develop a lattice path integral for the propagator, which we evaluate exactly for any finite volume. Our final formulas are similar, but not identical to the Yudson representation.
The paper is organized as follows. In 2 we construct a lattice path integral for the propagator and we express it in the so-called mirror or quantum channel using the elements of the quantum monodromy matrix. In 3 we introduce the so-called F -basis, which provides an essential simplification for the subsequent calculations. In 4 and 5 we treat the one and two particle propagators in detail. These sections are more technical, but both include a closing subsection where we summarize and interpret the corresponding formulas. Section 6 includes our main result for the general m-particle propagator. We conclude in 7, and some of the more technical calculations are detailed in the appendices A-D.

Lattice path integral for the propagator
We consider the spin− 1 2 XXZ model, described by the following Hamiltonian: Here σ x j , σ y j , σ z j are the usual Pauli matrices, acting on the jth subspace of the tensor product space ⊗ L j=1 C 2 . We assume periodic boundary conditions and use the following parametrization of the anisotropy parameter ∆: We do not restrict ourselves to any regime in ∆.
Our goal is to compute the real space propagator of the spin chain, which is defined as follows. First we choose the reference state to be the state with all spins down 1 We denote the basis state with m spins up at positions a j ∈ {1, . . . , L}, j = 1 . . . m by For the coordinate variables we always assume a j < a k for j < k.
The real space propagator is then defined as Note that the in and out states have the same magnetization; the spin-z conservation of the Hamiltonian implies that all other matrix elements of the propagator are identically zero. The propagator depends on the volume L, but for simplicity we omit this in the notation. 1 The reason for choosing the down spins as the vacuum is simply to conform with the conventions of the F -basis, to be introduced in Section 3. To put forward: A vacuum state with down spins will correspond to a product of D-operators, which are diagonal in the F -basis and thus easy to treat.
The propagator satisfies the Schrödinger-type equations and the initial condition In (2.7)Ĥ a,b are operators that act as the Hamiltonian on the corresponding coordinates. The equality between the second and third expressions in (2.7) follow from the fact that the Hamiltonian is a symmetric matrix in the spin basis: it can be seen from (2.1) that all its matrix elements are real, and a real Hermitian matrix is symmetric. In the following we discuss the symmetry properties of the propagator. It follows from the definition that However, a stronger condition also holds, the propagator is a symmetric matrix in the spin basis: This follows from the second equality in (2.7), or alternatively, from the fact that the exponentials of the Hamiltonian are also symmetric. Translational invariance and space reflection invariance lead to the conditions where we introduced the short-hand notations and similarly for {b}.
The XXZ model is well known to be exactly solvable by the different versions of the Bethe Ansatz [41]. In principle the propagator could be computed if we construct all eigenstates of the system. In the present work we will follow a different route, nevertheless we review this direct approach as well.
Denoting by |Ψ j a complete set of states we get where E j are the energy eigenvalues. The wave functions appearing in the sum above can be constructed by the coordinate Bethe Ansatz. Due to spin-z conservation it is enough to consider the L m eigenstates in the sector with m down spins. The states are characterized by a set of rapidities (also called Bethe roots) {u} m , that describe the pseudo-momenta of the interacting spin waves.
There are different conventions for writing down the wave function, here we use the formula (2.13) where the function P (u) describes the propagation of a single magnon and it is given by (2.14) and for later use we also introduced In order for such a state to be an eigenstate the rapidities need to satisfy the Bethe equations, that follow from the periodicity of the wave function: In these conventions the one particle solutions for ∆ > 1, η ∈ R are purely imaginary, whereas for ∆ < 1, η ∈ iR they are purely real. The energy eigenvalues are In these conventions the norm of the wave function for the eigenstates is [42] , (2.20) and G is the so-called Gaudin matrix: where Q j are the logarithms of the Bethe equations defined in (2.16) and .
The representation (2.12) could serve as a starting point to derive exact results for the propagator. The standard idea is to transform the summation over the Bethe states into multiple integral formulas, by using the Gaudin determinant as the multi-dimensional residue for the contour integrals around a particular set of Bethe roots. However, there are several difficulties with this approach. First, one needs to know all possible positions of the Bethe roots, and then to combine and/or transform the resulting integrals into some manageable form. Second, one needs to prove the completeness of the Bethe Ansatz, which is a notoriously difficult problem, especially for the homogeneous chain. Typically the spin chain has a number of singular solutions too, which need to be taken into account [43,44]. Some of these problems could be circumvented by introducing a twist and/or inhomogeneities along the chain, but certain technical difficulties (for example the treatment of the contours) would remain.
Instead, in the present work we pursue a completely different approach. We compute the propagator by a lattice path integral, in close analogy with the study of the thermodynamical state functions of the model [45]. We obtain multiple integral formulas as anticipated by the insertion of the complete set, but our calculation immediately gives us the contours for any particle number and any finite volume L. The relation with the explicit summation over intermediate states can be studied after the main results are found, an example of this will be discussed in Section 4.

Trotter decomposition and quantum transfer matrix
As usually, we consider the Hilbert space of the chain and an auxiliary space C 2 denoted with the index 0 and construct the monodromy matrix of the spin chain as The transfer matrix is given by τ (u) = Tr 0 T (u).

(2.24)
Here R jk (u) is the so-called R-matrix acting on the spaces j and k. We use the normalization where b(u) is given in (2.15) and c(u) = sinh(η) sinh(u + η) . (2.26) The function b satisfies the relation b −1 (u) = b(−u − η) which will be used often in this work. The R-matrix satisfies the Yang-Baxter (YB) equation and the unitarity and crossing relations where t 1 denotes transposition in the first vector space. Further, the R-matrix satisfies the initial conditionR(0) = P , where P is the permutation matrix of It follows from the YB relation that the transfer matrices with different spectral parameters form a commuting family: The transfer matrix satisfies the initial condition τ (0) = U , where U is the shift operator by one site to the left (towards decreasing lattice site indices). Also, it generates the Hamiltonian (2.1) through the relation . (2.29) For future use we also introduce the space reflected transfer matrix, which is defined as It satisfies the initial conditionτ (0) = U −1 , and it also generates the Hamiltonian by the same relation as (2.29). It follows that the two transfer matrices τ (u) andτ (u) can be used to generate the time evolution operator through a Trotter approximation. For any s ∈ C For our purposes it is convenient to use the crossing relation (2.28) to relate the two transfer matrices to each other: This leads to the following form of the Trotter decomposition: (2.34) The advantage of this representation is that it only uses the same transfer matrix. For technical reasons it is better to use a set of non-coinciding inhomogeneities. Therefore we define β j with j = 1, . . . , N such that |β j − β k | = O(1/N ). Focusing on real times we thus write The propagator is thus expressed as (2. 37) Due to the correspondence between the XXZ model and the six vertex model, at any finite N the above expression is equal to a six vertex partition function of an L × 2N square lattice, with particular boundary conditions. The individual weights of the six vertex model are depicted on Fig. 1, whereas an example for the partition function is shown in Fig. 2.
Here the presence of the transfer matrices leads to periodic boundary conditions in the space direction (chosen as the horizontal direction), whereas in the time direction (chosen as the vertical direction) we have fixed boundary conditions given by the initial and final states |a 1 , . . . , a m and b 1 , . . . , b m |.
The six vertex model is invariant under a reflection along the North-West diagonal. Due to this reflection symmetry of the R-matrix, any such partition function can be evaluated alternatively in the so-called quantum channel (also called the rotated or mirror channel). We can thus build an alternative monodromy matrix that acts on an inhomogeneous spin chain of length 2N , such that the inhomogeneities are determined by the spectral parameters of the transfer matrices in (2.37). To be precise we define (2.38) The initial and final states of the propagator become boundary conditions in the quantum channel, and (due to the periodic boundary conditions in the spatial direction) the partition function can be evaluated as a trace of a particular ordered product of the monodromy matrix elements A(0), B(0), C(0) or D(0). The explicit relation is where s a,b j are the the spin components at position j in the initial and final states, respectively. Explicitly and similarly for s b j . The normalized expression for the propagator is then As an example, we consider L = 6 and a specific matrix element of the two particle propagator: According to (2.39), this is proportional to Tr B(0)C(0)A(0)D(0)D(0)D(0) . We note that the ordering of the inhomogeneities in (2.38) does not influence the traces that we intend to compute. On the one hand, this follows from the commutativity of the transfer matrices τ in (2.34). On the other hand, this symmetry will be explicit after the introduction of the F-basis in Section 3.
The symmetry (2.10) of the propagator can be observed at finite Trotter number as well. Starting from the expression (2.38) for the quantum monodromy matrix we can perform a series of crossing transformations on the R-matrices leading to Note that in (2.44) the order of the inhomogeneities has been modified as a result of the crossing, but this does not effect the traces. Also, the action of the S operators also drops out due to S 2 = 1. Thus it follows from (2.43) that the initial and final states can be exchanged, and an alternative formula for (2.41) is Our goal is to compute the traces (2.41)-(2.46) at finite N using exact methods, and to take the N → ∞ limit afterwards. We will start with low particle numbers, but we will consider arbitrary L volumes. Typically there will be a large number of D operators in the product (corresponding to the vacuum with the down spins), and a smaller number of B, C, and A operators depending on the initial and final positions of the excitations.
We stress that in the usual QTM method one deals with the powers of Tr T QT M (0) = A(0)+D(0), and the partition function is computed only in the L → ∞ limit. This leads to the major simplification that only the leading eigenvalue of Tr T QT M needs to be considered. In contrast, here we are dealing with the individual matrix elements of the monodromy matrix, we keep the volume L at a fixed finite value, and evaluate the product (2.41)-(2.46) exactly.
The direct evaluation of the traces of the product of operators in (2.41)-(2.46) would be quite cumbersome, due to the complicated forms of the A, B, C, and D operators. A remarkable simplification can be achieved by performing an appropriate basis transformation. We compute the traces in the so-called F -basis, whose big advantage is that the D operators are diagonal and the B, C, and A operators take also sufficiently simple forms. This way we obtain manageable expressions for the propagator.
In the remainder of the paper we will work with (2.46) instead of (2.41). The reason for this is simply that in the conventions that we are using it leads to more transparent prescriptions for constructing the propagator.

F -basis
The factorizing F -matrices of Maillet and Santos were introduced in [46] and later used in [47] for the calculation of the form factors in the spin chain. From a computational point of view, their main advantage is that in the basis generated by the F -matrices (the so-called F -basis) the monodromy matrix elements take very simple forms and their expression is completely symmetric with respect to the ordering of the quantum spaces of the chain. In the following we give a brief overview on the F -basis of the XXZ model. We will not repeat the derivations, rather we will just cite the formulas necessary for our computations. For a diagrammatic introduction to the F -matrices we refer the reader to [48], and we also note that an interesting alternative derivation is also given in [49].
With full generality we consider a spin chain of length n with a set of inhomogeneities ξ j , j = 1 . . . n. As usually, the monodromy matrix is For the main goal of this paper, namely the computation of the propagator, the F -basis will be established in the QTM channel and we will identify n = 2N , and the inhomogeneities will be set to Here we re-ordered the set of inhomogeneity parameters, but this does not effect the computation of the traces, as already remarked earlier. In the present section our goal is just to introduce the necessary formulas for the F -basis, therefore we keep n and ξ j as arbitrary parameters.
3) The factorizing F -matrix is an invertible matrix F 1...n (ξ , . . . , ξ n ) which satisfies the following condition for any π ∈ S n : In other words it factorizes the (composite) R-matrix. It was shown in [46] that the F -matrices can be constructed recursively, by increasing the length of the spin chain at each step. Alternatively, they can be computed by a summation over the matrices R π [49]. However, these formulas will not be needed in the present paper, therefore we refer the reader to the original papers.
The important property of the F -matrices that will be used in our calculations is, that the monodromy matrix elements take especially simple form in the basis generated by them. Let us perform a basis transformation in the physical space (and keep the auxiliary space unchanged), and let us denote byT the monodromy matrix in the new basis: It follows from the relation (3.4) that the matrix representation in the new basis is completely symmetric with respect to the sites and inhomogeneities. It was shown in [46] that the F -matrices given there lead to the following diagonal form for theD-operator:D Note, that this expression is completely symmetric under simultaneous permutation of the vector spaces i and the corresponding spectral parameters ξ i . Analogously, the other elements of the monodromy matrix are given bỹ (3.7) (3.8) Here σ ± i are the usual spin raising and lowering operators acting on site i: σ ± = 1 2 (σ x ± iσ y ). Note, that whileD,B,C are formally the same as in [46] ,Ã is formally different. The reason for this is that in [46] only the SU (2)-symmetric XXX− 1 2 case was considered, where there are some special identities for the rational functions involved. On the other hand, the formula forÃ can be computed using the quantum determinant (this was already remarked in [47] and the computation was performed earlier by other researchers [50]) or by taking the limit of the formulas for the XYZ model published in [51]. For the sake of completeness we present the detailed derivation in the Appendix A, together with an alternative form (A.5).
The advantage of the F -basis is that it provides polarization-free expressions for the B, C and A operators. By this we mean that the action of the particle creation and annihilation is dressed only diagonally. In contrast, in the original spin basis we would be dealing with expressions of the type where Ω i , Ω ijk , . . . are diagonal operators on all sites but on site i, site i, j, k, etc., respectively. Multiplying such sums would be a practically unfeasible task. On the other hand, the computation of the products of diagonally dressed operators is relatively straightforward. We also remind that theÃ,B,C,D operators satisfy the same commutation relations as A, B, C, D.
In the next two sections we employ the F -basis to compute the traces of products of monodromy matrix elements to obtain the propagator as given by (2.46).

The propagator: One particle case
In this section, we compute the one particle propagator. This is a very simple object, which could be calculated even without using any methods of integrability. Nevertheless we perform the detailed computations using the F -basis, which serves as a good warm up for the more complicated cases.
According to the rule (2.46) there are two different cases, that need to be treated separately: • The particle moves by ℓ sites with l = 1 . . . L − 1. In this case we are dealing with a trace of the form The simplest case is when ℓ = 1, this is treated in 4.1, whereas the generic case of ℓ = 2 . . . L − 1 is considered in 4.2.
• The particle stays at its position. In this case we are dealing with a trace of the form For the computations of the traces we will use the F -basis; the original spin basis of the QTM will not be used anymore. Therefore, in the rest of the paper it is understood that all A, B, C, D operators are given by their concrete matrix representations in the F -basis. Also, all monodromy matrix elements will be evaluated at the spectral parameter u = 0, therefore we will drop this from the notation and understand that A ≡ A(0), etc.
The computations of the traces will be performed algebraically using arbitrary inhomogeneities ξ j , j = 1 . . . n. We will use the following notations: .
When turning to the Trotter limit, we will set n = 2N and the inhomogeneities will be specified according to (3.2). Finally, the β j parameters will be sent to β = 2 sinh(η)t with t being the physical time parameter.

One particle propagation by one site
In our framework the simplest case is when one particle hops one site. According to (2.46) this is evaluated as First, we compute the D L−2 BC product. We can denote the structure of B and C as are diagonal matrices in the jth space, depending on index i, through its parameters. The product of these two is: where we separated the terms depending on whether the off-diagonal matrices are on the same site (first sum) or not (second sum).
Hence the BC product: .

(4.4)
As D is diagonal, multiplying by D L−2 is simple: . (4.5) Taking the trace gives Note, that the second sum from (4.5) dropped out automatically due to its tracelessness. We substitute this expression back to the propagator: It is easy to see, that this expression is singular due to b(0) = 0, and the homogeneous limit of ξ j → ξ k can not be taken directly. This was the main reason behind the introduction of the set of non-coinciding β j parameters. Nevertheless the expression can be evaluated as a contour integral, by the use of the following identity.
Consider the meromorphic functions f (u) and g j (u), and a contour C such that each g j has a simple pole at z j within the contour and f does not have poles inside C. Then In our case the set of inhomogeneities ξ j , j = 1 . . . 2N is given by (3.2). It follows that the corresponding functions will have poles in the neighborhood of z = 0 and z = η. Therefore we define C to be an union of two small contours around 0 and η: C = C 0 ∪ C η . We apply the above identity with the choice One can easily see that f is indeed free of poles in C. This leads to the integral representation (4.10) After substituting (3.2) we are now free to take the homogeneous limit. This leads to (4.11) To take the Trotter limit, note that each term involving sinh L (β/N ) is sub-leading, irrespective of the value of u under the integral. Therefore we get Taylor expanding the last factor to first order leads to (4.13) Substituting this back to the propagator, we get the final contour integral expression (4.14) An interpretation of this formula together with the remaining cases will be given in 4.4.

One particle propagation by ℓ sites
We consider the propagator describing one particle moving ℓ sites, and the corresponding partition function in the Trotter decomposition: The computation of the product is very similar to the previous case, leading to the following result: . (4.16) For the trace we get: This sum can be turned into a contour integral using the method of the previous subsection. We can apply the identity (4.8) with the same g j functions, but the following f : (4.18) It can be seen that f (u) is always regular at u = 0, and it is regular at u = η if ℓ + 1 ≤ L. Therefore the contour integral representation based on (4.8) is valid only for 1 ≤ ℓ ≤ L − 1. This constraint is in agreement with our general picture about the lattice path integral: the propagator with ℓ = L corresponds to a particle staying at its position, due to periodicity. In this case a different trace needs to be evaluated which includes an A-operator, and this is treated in the next subsection.
Repeating the steps of the previous subsection we obtain for ℓ < L

One particle not moving
Here we consider the case, when the particle stays in place. For this matrix element we have from (2.46): The A operator takes the following value at u = 0: . (4.20) Out of the three terms of the A-operator, two contribute to the trace: Note that the second term coincides with the expression of the previous subsection with ℓ = 0. This suggests a close relation between the two cases, namely that the case of a particle staying at its place should be given directly by same formula, where ℓ = 0. In the following we will see that this is indeed true, with the addition that the first term in (4.21) is also needed to produce the correct contour integral.
Consider again the product f (u) j g j (u), but in the present case assume that f (u) has a simple pole atz within the contour C. Then we get We choose It can be seen that f has a simple pole at u = 0 with residue one. Therefore, the identity (4.22) immediately gives the right hand side of (4.21).
Performing the Trotter limit as before we obtain the final contour integral (4.24) As anticipated, this is formally identical to the result of the previous subsection with displacement ℓ = 0.

Summary and interpretation of the one particle formulas
The interpretation of the contour integrals can be given if we introduce a shift of −η/2 in the integration variable. Correspondingly, we also introduce a contour C ± that is a union of two small circles around the points ±η/2: C ± = C −η/2 ∪C η/2 . Then the results of the previous subsections can be summarized as follows: Here q(u), P (u) and ε(u) are defined in (2.20), (2.14) and (2.18), respectively. These functions are the same ones that appeared in Section 2 in the coordinate Bethe Ansatz representation of the eigenstates of the model: P (u) describes the propagation of a single magnon, whereas ε(u) describes the one-particle energy (note that we substituted the relation β = 2 sinh(η)t).
The denominator above is simply the one particle Bethe equation. In this simple case this denominator is the only part of the formula that depends on L.
The appearance of the elements of the coordinate Bethe Ansatz is of course not a coincidence: in the following we show that (4.25) can be transformed into a sum over the physical eigenstates, just as written in (2.12). We also note that a representation similar to (4.25) is also obtained in the Yudson-representation [36], but without the presence of the Bethe equation.
In the remainder of this subsection we restrict ourselves to ∆ > 1, η ∈ R, but a similar analysis can be performed also for ∆ < 1. We consider the analytic structure of the integrand in (4.25). This function is periodic in the direction of Im u, with period π, hence, we focus on the strip −π/2 ≤ Im u ≤ π/2. There are singularities at the following locations within this strip: 1. Essential singularities of exp(−iε(u)t) at u = ±η/2. These arose in the course of the Trotter, and they are inside the contour C ± . 2. Poles at u j = iλ j , j = 1 . . . L with λ j ∈ R corresponding to the solutions on the one-particle Bethe equations These are outside the contour.
There are no other poles within this strip. We modify the contour as follows: We blow up the two small circles around ±η/2, until they meet in the middle around the Bethe roots. In the imaginary direction, they are enlarged until Im u = ±π/2, and we take the remaining side of the contours to ±∞. This contour modification is shown in Fig. 5. The integrals on the horizontal segments cancels each other due to periodicity. The integrals on the vertical segments [−u − iπ/2, u − iπ/2] and [−u + iπ/2, u + iπ/2] with u → ∞ tend to zero, because q(u) → 0 and the remaining functions remain regular. Hence the integration on C ± can be rewritten as a sum over integrals around the Bethe roots: 27) iπ/2 −iπ/2 −η/2 η/2 Figure 5: A contour deformation that transforms the integral (4.25) into a sum over Bethe roots. The solid circles make up the original contour C ± , whereas the crosses denote the positions of the one-particle Bethe roots.
Here C uj 's denote small contour enclosing the u j one particle Bethe roots. Note that due to change of direction of the contour, there is an overall (−1) factor. It follows from (2.20) that at each Bethe root we have This gives (4.28) Consider the one particle Bethe vectors denoted by |u j , j = 1 . . . L. According to the general formula (2.13) their wave functions are simply with the energy E j = ε(u j ) and their norm given by Ψ|Ψ = L. It follows that the propagator is indeed expressed as k|u j e −iε(uj )t u j |k − ℓ u j |u j . (4.29) In this simple case we could have started with the representation (4.29) to derive (4.25) afterwards. However, the advantage of the seemingly long derivation is that formulas analogous to (4.25) will be obtained also for the higher particle cases, using the same contours. This would not have been evident had we started only from the particle picture, especially due to the complications with the string (and modified string) solutions.
For the sake of completeness we show that the original contour integral (4.25) is compatible with the space reflection symmetry of the model. Consider the transformation u → −u for the contour integral. Under this transformation C ± is fixed but he integration measure du gets an overall (−1) factor. Also, Therefore we get Thus we have indeed obtained the space reflection symmetry. For the one-particle propagator the relation (4.31) coincides with the symmetry property (2.10), but for higher 5 The propagator: Two particle case In this section we derive the two particle propagator. The computation is similar to the one particle case, nevertheless it poses some additional difficulties. We put forward that the final result displays all features of the generic multi-particle propagator.
In order to simplify the notations and later computations we introduce generalized operators that involve the action of D-operators from the left: Their explicit forms are We will also use a simplified notation for the generalized operators, whenever we want to suppress the notation of the number of inserted D's: When considering products of generalized matrices in this notation, we assume, that the number of inserted D's is generic and it can be put back to the formulas whenever needed. In order to compute the propagator a, b|e −iHt |c, d , we need to distinguish different cases depending on the relative position of the coordinates. According to the rule (2.41) and the cyclicity of the trace there are five different possibilities, that correspond to specific traces as follows.
• For c < a < d < b one needs to compute Tr B (ℓ1) C (ℓ2) B (ℓ3) C (ℓ4) , where • For a = c < b < d one needs to compute Tr A (ℓ1) C (ℓ2) B (ℓ3) , where • For a = c < d < b one needs to compute Tr A (ℓ1) B (ℓ2) C (ℓ3) , where • For a = c < b = d one needs to compute Tr A (ℓ1) A (ℓ2) , where For the following computations we introduce the notion of the contracted traces, which are similar to (but not identical with) the contractions used in the Wick theorem in free field theory. The motivation for the definition comes from the form of theX operators: they are sums of products of operators such that in each product there are only one or two operators that act non-diagonally. After multiplying allX and expanding the product, one can keep track of these non-diagonal operators by listing their indices, i.e. on which C 2 subspaces of the QTM they act. The contracted trace is a particular sum of the traces of these products, where we sum over a specific index pattern while excluding coinciding indices. In the following we give examples for this idea.
For simplicity, we first consider the B (ℓ) and C (ℓ) operators. Only those terms of the product of B (ℓ) 's and C (ℓ) 's have non-vanishing trace, where an equal number of σ + and σ − share the same indices, and are ordered alternatingly. The simplest contracted traces are (using the notation (4.2)) (5.10) Here we take those n terms from the expansion of theBČ (orČB) products where the non-diagonal part acts on site i, and afterwards we sum over these traces.
In the case of four operators we have more possibilities, for example k .

(5.11)
In the first case we take those n terms where each non-diagonal part acts on site i, whereas in the second case the non-diagonal pieces act on sites i and j and we require i = j. This distinction is important: if we evaluate such a term for generic i, j and set i = j afterwards, we obtain a result that is different from the direct evaluation of the first case. Our strategy in the calculations is that we keep track of all possibilities explicitly and transform the sum of all terms into a contour integral. A more complicated example of a contracted trace is l .

(5.12)
Incorporating theǍ operators into this framework is straightforward based on (3.8):Ǎ can be regarded as a sum of a diagonal term (first line in (3.8)) and aBČ product (second and third lines in (3.8)). When considering a contracted trace involvingǍ we have to specify, whether we refer to the diagonal or theBČ part. For terms that involve the diagonal piece we attach an empty set of indices to theǍ operator and it will be denoted asǍ () . For the non-diagonal terms there will be two indices attached toǍ that specify the sites on which the two non-diagonal operators act. Examples for this will be shown in the subsections below and in Appendices C-D.
After these considerations we can give the definition for the contracted trace: The contracted trace is a sum of particular terms in the expansion of the trace of the product ofX's in the F -basis, and it is characterized by a specific index pattern. We sum over those terms in the expansion, where the off-diagonal factors of the various operators act on the subspaces designated by the index pattern. Furthermore, we exclude those terms from the sum where different indices would take coinciding values.
In this section we consider two configurations for the two particle propagator in detail (those corresponding to TrBČBČ and TrǍČB); these two cases showcase all the computational nuances, which arise in the two particle case. Detailed calculations in the remaining three other cases are presented in B-D. A summary of the two-particle propagator is given in subsection 5.3.

The TrBČBČ case
In the case of TrBČBČ there are three non-vanishing contracted traces: These are computed as (5.14) Here we used ℓ 1 + ℓ 2 + ℓ 3 + ℓ 4 + 4 = L. For their sum we thus get The Trotter limit is taken after the identifications (3.2). As in the one particle case, this expression suffers from singularities: The b −1 ij type terms are singular in the β i → β, ∀i limit.
To overcome this, we use the same trick as for the one particle case. Consider the following double contour integral: , (5.16) where f and h i are meromorphic functions, which do not have poles inside C, and 1/g i have simple poles at ξ 1 , . . . , ξ n located inside the contour C. The integral can be evaluated by the successive application of the uni-variate residue formula: .
The second summand of this expression is equal to TrBČBČ (iijj) + TrBČBČ (ijji) if we specify the functions as and S(u) given by (2.17). It is easy to see that f is indeed free of poles, given that −1 < x < L − 1 and −1 < y < L − 1.
If we substitute back these functions into the first summand of (5.17) then we obtain the first term of (5.13): Thus the double integral reproduces all three terms of (5.13) and we can write .
The propagator is obtained after we include the normalization factors: , (5.23) valid for c < a < d < b.
The Trotter limit is taken analogously to the one particle case, leading to the result Note, that taking the Trotter limit is completely independent of the function f . A more transparent result is obtained if we express the ℓ i , i = 1 . . . 4 with the original coordinates a, b, c, d according to (5.5). Similar to the one-particle formula (4.31) we introduce a shift of η/2 and use the functions q(u), P (u) and S(u) (given by (2.20), (2.14) and (2.17)) to express the propagator as This is a remarkably simple form, with a clear physical interpretation to be discussed in 5.3.

The TrǍČB case
Here we compute the case TrǍČB describing the propagator for a = c < b < d. The following contracted traces give non-vanishing contributions: (5.26) As explained above, the empty indices () for theǍ operator denote that for term we take the diagonal part ofǍ, whereas in the other two cases we take the non-diagonal part and the two indices denote the sites on which the non-diagonal factors act. These three contracted traces are computed as follows. We will use the relation ℓ 1 + ℓ 2 + ℓ 3 + 3 = L relevant to this case.
In the first equation we used b ℓ1+ℓ2+ℓ3+2 The full partition function is thus This expression is similar to (5.15), but we notice a few differences: • There is no TrǍČB (iiii) term, because this contracted trace is automatically zero. By its structure this term would correspond to TrBČBČ (iiii) from the previous case.
• On the other hand, there is an extra term TrǍČB (()ii) which does not have corresponding term in the TrBČBČ case. Based on the one particle case (Section 4.3) we can anticipate the role of this term on the contour integral side: The f function will not be free of poles within C, and considering the pole of f will lead to this term.
• In the previous case only symmetric products of b −1 ij b −1 ji occurred, for both TrBČBČ (iijj) and TrBČBČ (ijji) . However, in this case there is also a non-symmetric b −2 ij present. All these differences find an explanation as we transform the sums into a common contour integral.
Consider again the integral gi(u1)gi(u2) with the assumption that f (u 1 , u 2 ) has a pole at u 1 =z: .
The second summand can be identified with TrǍČB (iijj) + TrǍČB (ijij) if we specify sym. (u 1 , u 2 ) and S(u) given by (5.19), (5.21) and (2.17) respectively. The difference between f non−sym. and f sym is merely the extra factor of S(u 2 − u 1 ). This factor has a pole for u 2 − u 1 + η = 0, and such points are included in the double contour integral. However, this pole is canceled by the same factor appearing also in the denominator, therefore it does not give any new terms.
After this identification the first term of (5.31) is evaluated as Here we used the identity S(0) = −1. The main reason for the vanishing of this expression is that for this particular contracted trace the functions f non−sym. and f sym. cancel each other. Finally we treat the third term in (5.31). The function f (u 1 , u 2 ) is singular at u 1 = 0 with the residue given by Collecting all factors we can see that the third term in (5.31) reproduces the term TrǍČB (()ii) in (5.26). Thus we get the full equality: .
Taking the Trotter limit follows the same steps as previously. We introduce the shift of −η/2 once again and obtain the final result (5.38)

Summary of two particle case
In the previous subsections we have computed two out of the five cases for the two-particle propagator. The other three cases can be treated similarly, and the detailed computations are presented in appendices B to D. The common properties of these calculations are the following: The partition function at finite N is transformed into a double contour integral. The h k and g k functions are the same in all cases, and f depends on the specific configuration. The various contracted traces correspond to various terms in the sum over residues. Up to now we have performed the identification on a case by case basis; a general combinatorial proof for higher particle cases will be presented elsewhere.
The derivations lead to the following general form: .

(5.39)
Here Ψ {a,b},{c,d} (u 1 , u 2 ) is an amplitude, with the explicit form in the five different cases being These five formulas emerged from the concrete computations, and they can be given a common interpretation as follows.
We can see that the amplitude Ψ {a,b},{c,d} (u 1 , u 2 ) is reminiscent of the Bethe Ansatz wave function, but it is not identical to it. Most importantly it depends on two sets of coordinates: the initial and the final ones. It is given by a sum over permutations, where we permute the final positions of two particles, which are started from positions c and d and are indexed with their rapidity parameters u 1 and u 2 , respectively. For each permutation there is an assigned phase, which consists of the one-particle propagation phases and it can also include a scattering phase. We observe that the phases P ℓ (u) with some ℓ associated to the one-particle propagation are such that each particle always travels to the right, and if its final position is to the left of the initial one, then it is required that the particle travels around the volume. It is important that this rule is not imposed by any Ansatz, rather it emerges naturally from the computation. Also, this rule of "moving to the right" is not physical and it is only used to construct Ψ {a,b},{c,d} (u 1 , u 2 ). Regarding the scattering phases we can observe that a factor of S(u 2 − u 1 ) is inserted if during that particular displacement process if there is a crossing of the world lines of the two particles. A pictorial interpretation of these displacement processes and the construction of the amplitude is given in Fig. 6.
The case of c < a < d < b: The case of c < a < b < d: The case of a = c < b < d: We remind that the five cases listed above do not exhaust all possible initial and final positions, and all other possibilities follow from periodicity, which is implied by the periodicity of the traces (2.41).

The multi particle propagator
Here we present our main result for the m-particle propagator in the form of an mfold multiple integral. The formula is a straightforward generalization of the two particle propagator, which already includes all building blocks of the general case. We constructed a general proof of this result, which will be presented elsewhere due to its length and technical nature.
The following is our main formula for the m-particle propagator: . (6.1) Here the contour C ± consists of two small circles around the points ±η/2, the functions q(u), P (u), S(u) and ε(u) are given by (2.20), (2.14), (2.17), (2.18). Notice that the denominator is simply the product of the m-particle Bethe equations. The amplidute Ψ {b},{a} (u 1 , . . . , u m ) is a quantity analogous (but not identical to) a Bethe Ansatz wave function; it is defined below. First we give a set of rules for constructing Ψ {b},{a} (u 1 , . . . , u m ), these rules describe a displacement process of m particles. Afterwards we also give an explicit formula for it.
Consider the process of moving m particles with rapidities {u 1 , . . . , u m } from their initial positions {a 1 , . . . , a m } to a set of final coordinates {b 1 , . . . , b m }. The particle with rapidity u j starts from position a j , and (in analogy with the Bethe wave function) the amplitude Ψ {b},{a} (u 1 , . . . , u m ) will be given as a sum over permutations corresponding to the different possibilities of how the particles can occupy the final positions. For each particle with rapidity u j there will be an associated one-particle propagator of P ℓj (u j ), where ℓ j is the distance that this particle has traveled. However, here we require that every particle only moves to the right (in the direction of increasing x coordinates along the chain) or it can also stay at its place. If a particular final position is to the left of the initial one, then we require that the particle travels "around the volume". It is important that this rule of "moving to the right" is used only for the construction of the amplitude, and it is not to be confused with the real speed of some physical particles, which can be negative as well. Finally, to each of these processes we associate a scattering phase, which is defined analogously as in the Bethe Ansatz: if there is a crossing of two particles during the displacement process then we add a factor of S(u j − u k ), where particle j was initially to the left of particle k. Once the initial and final positions of the individual particles are specified, then the total phase factor is also unambiguously determined. Note however, that in the typical case this factor differs from the usual Bethe Ansatz scattering phase, due to the requirement of moving always to the right, with periodic boundary conditions.
These rules lead to the following formula: and the generalized coordinates d j reflect the requirement of the particles moving to the right: The definition ofS and the fact that we take the product over j < k in (6.2) guarantee that we obtain the correct scattering factors. We stress that P (u), S(u) and ε(u) are the same functions that appeared in the coordinate Bethe Ansatz description of the model in Section 2. Furthermore, the denominator of (6.1) is simply the product of the m-particle Bethe equations. However, the amplitude Ψ {b},{a} (u 1 , . . . , u m ) differs from the original coordinate space Bethe wave function (2.13), and it also differs from the amplitudes used in the Yudson-representation [36]. Its main features are the following: • There is a direct dependence on the initial coordinates {a 1 , . . . , a m } through the factors of P −aj (u j ). Note that these factors could be decoupled, because they do not depend on the permutation or the extra S factors. In this sense, the "in" side of the amplitude has the form of a free wave function. This is similar to the Yudson-representation [36], where the two sides of the propagator are decoupled and one side is simply a free wave function. • The dependence on the final coordinates {b 1 , . . . , b m } is reminiscent of the usual Bethe wave function, with the slight modification that the particles always move to the right, and this way the relative phases have a second, indirect dependence on the initial coordinates. To be precise, the phases depending on the permutation σ and b j can be P (σb)j (u j ) or P (σb)j +L (u j ), for (σb) j ≥ a j or (σb) j < a j , respectively. This is a new feature, which is not present in the earlier wave function amplitudes found in the literature. • The amplitude Ψ {b},{a} (u 1 , . . . , u m ) is manifestly periodic for arbitrary sets of rapidities. This is also a new feature, that arose naturally from our computation, which is manifestly periodic at each step.
If we set the rapidities (u 1 , . . . , u m ) to be a solution to the Bethe equations (2.16), then the amplitude Ψ {b},{a} (u 1 , . . . , u m ) factorizes into a product of a regular Bethe wave function (in the b j -coordinates) and a free wave function (in the a j -coordinates). This follows simply from the compatibility of the P and S factors in our construction. The Bethe equations mean that as a particle moves around the volume, the phases coming from the single particle propagation and the two-particle crossings add up to zero. If this holds for all particles, then it does not matter whether we construct the relative phases by the rule of "always moving to the right" or in the usual way. We stress however, that the contour integrals over C ± always avoid the solutions to the Bethe equations.
We also remark that for generic rapidities the amplitude Ψ {b},{a} (u 1 , . . . , u m ) is not an eigenstate of the Hamiltonian when viewed as a function of the b-coordinates. It is manifestly periodic, but it has discontinuities at b j = a k for some j, k. The Bethe equations could be recovered by recovering continuities at these special points. This can be compared with the usual situation where the Bethe equations follow from the periodicity requirement.
All of these similarities with the usual coordinate Bethe Ansatz are of course not a mere coincidence. We expect that the multiple integral can be transformed into a sum over all Bethe states using specific contour deformations, in analogy with the example of the one-particle case presented in 4.4. This provides the physical interpretation of our contour integrals: they replace the summation over the eigenstates for all m and L. As remarked above, the amplitude Ψ {b},{a} (u 1 , . . . , u m ) factorizes for Bethe roots, thus we would obtain a summation similar to (2.12), where one side would be replaced by a free wave function. This would be analogous to a "finite volume Yudson representation". The study of such contour deformations is left to further research.
To conclude this section we give a technical remark about the A operators. In the calculations above we used the representation (3.8) of A in the F-basis. We could have performed the computations with the alternative representation (A.5), and we would have obtained a different rule for constructing the Ψ {b},{a} amplitudes. In this modified representation the propagation phase associated to a particle staying at its place would be P L (u) instead of P 0 (u) = 1. This corresponds to moving around the volume, and some of the scattering factors would also change according to the modified particle crossings. We chose to work with the form (3.8) because having a phase of 1 for a particle not moving is more natural. Nevertheless both representations are correct and they become equal if we specify the set of rapidities to Bethe roots.

Conclusions and Discussion
We have obtained a multiple integral representation for the propagator of the finite volume XXZ chain. Our main result is formula (6.1), which is a compact expression that uses welldefined functions and contours for every L and m, such that the volume enters simply as a parameter. As we explained in the previous section, our solution bears many similarities with the Bethe Ansatz representation of the eigenstates, but there are important differences too. Perhaps the most interesting feature of our result is that the amplitude Ψ {b},{a} (u 1 , . . . , u m ) appearing in the multiple integral is manifestly periodic for an arbitrary set of rapidity parameters. It is important that our solution is not an Ansatz, rather it appeared naturally from the lattice path integral and the use of the F -basis to evaluate it.
In the present work we treated the one and two particle cases in detail, and the main result was presented as a straightforward generalization. We have also constructed a rigorous proof of (6.1) for any m > 2 building on the same starting point (2.41). However, the general proof turned out to be quite lengthy and technical, and it will be published elsewhere.
Our derivations display a number of seemingly miraculous coincidences, that allow for the contour integral representation of the various traces. We observed that already at finite Trotter number N each partition function (corresponding to the specification of the initial and final coordinates) is given by a multiple integral of the form (6.1) with the same amplitude Ψ {b},{a} (u 1 , . . . , u m ). The only effect of the Trotter limit is to produce the correct timedependent factors of e −it j ε(uj ) . The actual computation of the traces in (2.41) and the intermediate sums varies from case to case, but the final contour integral is always of the same form. This suggests that the multiple integral formula could be proven by alternative means. A promising idea is to consider an inhomogeneous finite spin chain at finite N and to set up recursion relations for the partition functions in the spirit of [52]. We leave this to further research.
Having found a compact representation for the propagator, it is also important to discuss the practical applicability of the result. The m-particle propagator is given by an m-fold integral, and numerical implementations become very quickly unfeasible as we increase m. At present it seems that for most numerical purposes the known approximate methods (for example t-DMRG or ABACUS) or even exact diagonalization would perform better than the numerically exact evaluation of the multiple integral. Advantages of our representation could show up in the study of the long time limit or finite size effects at large L.
The propagator is an intermediate object that can be used to compute the time dependent local observables through the sum This is an alternative to the usual spectral representation: here the (double) sum over the Bethe states is included in the exact propagator, and the remaining task is to perform the real space summations. In many practical applications (in concrete quench protocols) the initial states take simple forms in the real space representations, thus the only challenging task is to perform the inner sums that connect the two propagators to the local operator. Here the explicit form of the amplitude Ψ {b},{a} (u 1 , . . . , u m ) appearing in the multiple integral needs to be analyzed, possibly in connection with known form factor calculations. We hope to return to this problem.

A The A operator in the F -basis
Here we compute the A-operator in the F-basis. As it was already remarked in [47], this is easily achieved by using the so-called quantum determinant. We believe that this result is not yet published in the literature, nevertheless the calculation is by no means new, it was already performed by other researchers [50].
The quantum determinant is a specific combination of the A, B, C, and D operators which is proportional to the identity operator. There are in fact four different relations leading to the same scalar factor [41]: All four relations could be used to yield a formula for A(u). The four relations lead to two different expressions for A.

First expression for A in the F -basis
First we choose to compute it using the last relation: Here We compute the CB product as Multiplying this with D −1 is relatively straightforward, as D −1 is diagonal: where we used twice that c(u + η, ξ)b −1 (u + η, ξ) = c(u, ξ).
Adding both terms in (A.2) leads to the first expression for A in the F -basis: (A.5)

Second expression for A in the F -basis
One can also use the relation D(u)A(u + η) − B(u)C(u + η) = n j=1 b(u − ξ j ) · I (A. 6) to compute A in closed form. We get First we compute D −1 (u − η): Hence the diagonal term: .9) and the non-diagonal part: (A.10) Using the identity c(u − η − ξ i ) = c(u − ξ i )b −1 (u − ξ i ) we obtain (3.8).

B The TrBČČB case
Here we treat the two-particle propagator in the case of c < a < b < d, when the trace to be computed is TrBČČB. The following contracted traces are non-vanishing: (B.1) These leads to the following f function (the g k and h k functions are the same): f (u 1 , u 2 ) = f (ℓ2,ℓ1+ℓ2+ℓ3+2) non−sym.

C The TrǍBČ case
Here we treat the two-particle propagator in the case of a = c < d < b, when the trace to be computed is TrǍBČ. The following contracted traces are non-vanishing: This case is similar to the detailed TrBČBČ case, with one further term (TrǍBČ (()ii) ). The corresponding f function is the following: sym. , f is singular, and the residue of f corresponds to the extra term TrǍBČ (()ii) . This leads to the final expression: a, b|e −iHt |c, d = C±×C± du1du2 (2πi) 2 P b−d (u 2 ) + P b−c (u 1 )P L+a−d (u 2 ) × × q(u 1 )q(u 2 ) exp (−i(ε(u 1 ) + ε(u 2 ))t) (1 − P L (u 1 )S(u 1 − u 2 ))(1 − P L (u 2 )S(u 2 − u 1 )) .

D The TrǍǍ case
Here we treat the two-particle propagator in the case of a = c < b = d, when the trace to be computed is TrǍǍ. This case has the most non-vanishing contracted traces, namely the following ones: This case is similar to the previous ones, however, it brings a minor novelty: the TrǍǍ (()()) contracted trace is new. The f function is the following: f (u 1 , u 2 ) = f (−1,−1) sym.
Note, that due to the product structure in all the relevant places, every residue can be taken successively.