Effective free-fermionic form factors and the XY spin chain

We introduce effective form factors for one-dimensional lattice fermions with arbitrary phase shifts. We study tau functions defined as series of these form factors. On the one hand we perform the exact summation and present tau functions as Fredholm determinants in the thermodynamic limit. On the other hand simple expressions of form factors allow us to present the corresponding series as integrals of elementary functions. Using this approach we re-derive the asymptotics of static correlation functions of the XY quantum chain at finite temperature.


Introduction
Exactly solvable one-dimensional quantum mechanical systems of interacting spins, bosons, and fermions provide a unique platform for studying non-perturbative effects. The algebraic and coordinate Bethe ansatz allow one to find the wave functions [1] and analytically address the thermodynamic properties of these systems [2]. The matrix elements of physical operators can also be found analytically in many cases [3][4][5][6], but the computation of correlation functions still remains quite challenging. For the vacuum correlation functions there are effective numerical methods based on integrability [7]. The asymptotic behaviour of the correlation functions can be investigated by means of effective field theory (Luttinger liquid) [8]. The origin of this behavior has been linked to the finite-size scaling of the matrix elements computed by means of the Bethe Ansatz [9][10][11][12][13]. For dynamical correlation functions based on this approach the corresponding effective field theory bears the name of non-linear Luttinger liquid [14][15][16].
At finite temperature, or more generally at finite entropy (density of states), both the numerical and field theory approaches experience some difficulties. In the numerical approaches one has to scan a much larger portion of Hilbert space to saturate the sum rules, as the form-factors (matrix elements of the physical operators) decay exponentially with the systems size contrary to the power-law decays at zero temperatures (see for instance [17,18]). The field theory approach is based mainly on the linear spectrum for the soft modes (low-energy excitations) which is valid only for very low temperatures [19,20].
A more rigorous approach was developed to evaluate finite temperature correlation function in integrable lattice models of Yang-Baxter type, based on the Quantum Transfer Matrix (QTM) [21]. The notion of the thermal form factor was introduced [22], which turned out to be useful for the asymptotic analysis of twopoint functions [22][23][24]. In the scaling limit, thermal form factors also arise axiomatically in the context of Integrable Quantum Field Theory [25][26][27][28][29][30][31][32][33]. Less rigorous but numerically accurate approaches are based on the thermodynamic limit of the form-factors and restricting summation to a finite number of particle-hole pairs [34][35][36][37].
However, the complete understanding has not been achieved and recently the QTM methods were revisited to address correlation function for the XX spin-chain [38][39][40]. Moreover, new systematic approaches for correlation functions in the Ising model for low density [41] and in the Lieb-Liniger model for the strongcoupling expansions [42] have been proposed. The generalization of Smirnov's form factor axioms for the thermodynamic states has been formulated in Ref. [43] and successfully applied to the reconstruction of the generalized hydrodynamic description of the correlation functions [44,45].
In this work we develop a heuristic approach to address correlation functions at finite density of states. Our main motivating example is the XY spin chain in a transverse magnetic field. On one hand it allows to get exact answers for spin-spin correlation functions in terms of Toeplitz or Fredholm determinants and on the other hand matches the complexity of generic systems. As we have mentioned above, this complexity is combinatorial in nature and reflects the fact that each form factor for the thermal states is exponentially small so the number of relevant form factors is exponentially large. This makes direct computation of the corresponding sum for the correlation functions notoriously difficult and force researchers to focus at most on the two particle-hole excitations [34][35][36][37], consider semiclassical approximations [46] or develop other approximation schemes [41,42].
We deal with this problem in a different manner. Namely, to describe the spin-spin correlation function evaluated on a state with finite density of entropy (energy) we introduce effective form factors for the fermions with the modified phase shift that absorbs information about the state and significantly simplifies combinatorics of excitations making it essentially analogous to the zero-temperature case. Here we have to emphasize that the expressions of form factors was inspired by the XX spin chain [47][48][49][50], rather than genuine spin form factors in the Ising/XY models [51][52][53][54][55].
We focus on the static correlation functions for which we demonstrate that after complete summation of the effective form factors series and taking the thermodynamic limit the answer can be presented in the form of Fredholm determinants. For the proper choice of the phase shift the kernels in the Fredholm determinants differ from the exact ones [50] by the exponentially small (in distance between spin operators) terms. Conversely, by first taking the thermodynamic limit of the effective form factors and then performing their summation we manage to present the Fredholm determinants as integrals of elementary functions. This kind of asymptotic behavior for models in the continuum (not the lattice) arises similarly from the solution of the Riemann-Hilbert problem for operators acting on the whole real line [56]. This asymptotics was conjectured to be universal for correlation functions of any gapless model of statistical mechanics at any temperature and for an arbitrary coupling constant [57].
An important ingredient for our asymptotic analysis is the winding number of the state-dependent phase shift ν(q) defined as the difference across the Brillouin zone, namely ν(+π) − ν(−π) = δ ∈ Z. (1) We recover the correlation length in the lattice version of the asymptotics in Ref. [57] at δ = 1 and additionally give an analytic expression for the prefactor. For δ = 0 and δ = ±1 we derive asymptotic behavior for the correlation function in the XY spin chain at finite temperature and compare it with the known answer [58]. Different winding numbers correspond to different values of the magnetic field and anisotropy. The winding number |δ| ≥ 2 does not have a direct physical interpretation in this model, but we perform the asymptotic analysis anyway and find the results consistent with the generalization of Szegő formulas [59]. Moreover, we have observed a peculiar identity between Toeplitz determinants and Fredholm determinant of sine-kernel type with finite rank, which, to the best of our knowledge, is new where the operatorsŜ ν and δV ν are generalized sine-kernels that act on L 2 ([−π, π]) and are defined by their kernels and Notice that the right hand side of this formula can be also be presented as a Fredholm determinant but with modified shifted ν(k) [60]. For ν(q) corresponding to the XY spin chain Eq. (2) gives different representations of the spin-spin correlation function, the left hand side was obtained in [50] and the right hand side in [58,61]. The paper is organized as follows. In Sec. 2 we define the tau function together with the effective form factors and outline the derivation of the Fredholm determinant presentation resulting from the summation of form factors. The details of this derivation are presented in Appendix A. In Sec. 3 we study the thermodynamic limit of the form factors and argue for an explicit presentation of the form factors series as integrals of elementary functions. All necessary technical results are given in Appendices B and C. Sec. 4 deals with the application of the general formulas to the XY spin chain. In Sec. 5 we discuss connection of the general result to the Toeplitz determinant and relations such as Eq. (2). Sec. 6 concludes the paper and offers an outlook.

Effective form factors
We start with the formal definition of the static correlation function (tau function), as a form-factor series here the ordered set k = {k 1 , . . . k N +1 } consists of N + 1 distinct shifted momenta inside the Brillouin zone (k i ∼ k i + 2πn, n ∈ Z) each being a solution of the transcendental equation for a smooth function ν(k). This function plays the role of the phase shift and is assumed to be compatible with the Brillouin zone structure, i.e. ν(π) − ν(−π) = δ ∈ Z.
The integer δ is the winding number (index). One can easily argue that the number of solutions of Eq. (6) is L + δ, each root defined up to O(1/L 2 ) terms. The set q = {q 1 , . . . , q N } is an ordered set of N distinct solutions of equation Motivated by the spin form factors for quantum XY chain written in the XXO basis [50], we postulate the following form-factor where det D is a trigonometric Cauchy type determinant that can be presented in two equivalent forms Furthermore, since we do not specify the specific operator we will sometimes refer to Eq. (9) as to the overlap, and use this term interchangeably with form factor. We assume that the index is of order δ ∼ O(1) as both the system size and the number of particles are approaching the thermodynamic limit N → ∞, L → ∞ such that N/L = 1. In this case, the summation over q can be performed exactly (see Appendix A) and the result for the tau function reads where the determinants are taken in the space L 2 (S 1 ) and the corresponding operators are defined by their kernels with The diagonal terms k = q are understood as in L'Hopital's limiting procedure. Further, we impose the relation to present tau function as withŜ ν being a generalized sine-kernel This way, the remainderR =V −Ŝ ν consists of integrals in Eq. (14), which are exponentially suppressed 1 for large and positive x. Let us call the tau function with discardedR as τ S , namely This particular generalization of the sine-kernel is contained in the prefactor (e 2πiν(k) − 1) and allows one to describe a modification of the system from the vacuum state for which ν(q) is constant within the arc k ∈ [−k F , k F ] and zero everywhere else, to the finite-entropy state, where, for instance, for the thermal state of the fermionic system the prefactor would be proportional to the single-particle Fermi distribution function 2 . In Sec. 4 we relate this type of kernel to the static spin-spin correlations in the XY chain. Then ν(k) will depend not only on the state but also on the parameters of the model.
3 Thermodynamic limit and direct summation of form factors 3.1 Winding number δ = 1 In the previous section, we considered the summation of the form factor series and subsequent taking of the thermodynamic limit. This leads to the presentation of the tau function as a Fredholm determinant. The essence of this derivation, which is outlined in Appendix A, is that each momentum q i ∈ q was treated independently. In this section, we focus more on the detailed structure of the ordered sets q in the sum of Eq. (5). The total number of solutions of the equation e iqL = 1 inside the Brillouin zone is L, which can be presented as As we have already pointed out above, the number of solution of Eq. (6), depends on the winding number δ.
In particular, for δ = 1 there exist exactly L + 1 solutions inside the Brillouin zone If we choose the set k = {k 1 . . . , k L+1 } in Eq. (5) then summation over q will only involve one term q = {q 1 , . . . q L }. In the large L limit the corresponding overlap reduces to a constant which is slightly counterintuitive from the orthogonality catastrophe point of view [63]. The explicit value of this constant is given by Eq. (210). The difference of momenta in Eq. (5) can be evaluated in the large L limit as Combining these observations together we obtain explicit equality for the Fredholm determinants in Eq. (11), and approximation for the generalized sine-kernel It is interesting to note that only the periodic (i.e. having winding number δ = 0) part of ν(q) has entered the final answer.

Winding number δ = 0
For δ = 0 we proceed similarly to the previous subsection. This time however the maximal possible number of the k i ∈ k is L, so the maximal set q consists of N = L − 1 momenta. There are exactly L such sets and they can be parameterized by the position of the "hole" q (a) = {q 1 , . . . , q a−1 , q a+1 , . . . , q L }, a = 1, 2, . . . , L.
The overlap is given by The derivation can be found in Appendix C.3 along with the explicit expression for A[q a ] (see Eq. (224)). For a ∼ L and L − a ∼ L the last two factors cancel each other, which yields the following explicit expression On a technical side, we have used a variation of Sokhotski-Plemelj formula to transform the integral in the exponential. For a ∼ 1 and L − a ∼ 1, the overlap is still O(1/L), so we can replace the sum in tau function Eq. (5) by an integral with and Equivalently we may re-write Y 0 (x) as a contour integral in the variable z = e ik where the contour C > is a circle centered at the origin with slightly larger than unit radius and We assume that ν(q) is non-singular in the region of integration, thus S(z) is holomorphic outside the unit circle on the Riemann sphere, so the asymptotic for large positive integers x is defined by the analytic behavior of ν(k) in the upper-half plane. For example, if e −2πiν(k) is a meromorphic function (of z = e ik ) outside the unit circle in the complex plane having simple poles at z 1 , z 2 , . . . , with 1 < |z 1 | < |z 2 | < · · · , then for large x the leading contribution comes from the smallest pole Applying this formula together with Eq. (28), we have an asymptotic expression for the sine-kernel Fredholm determinant for δ = 0 Remark. Notice that even for δ = 1 one could have chosen k = {k 1 , . . . k L }. This would not affect the derivation of the Fredholm determinants, but instead of one term in the form factor series as in the previous section, we still get a sum of L terms. Using Appendix (C.3) and specifically Eq. (226), we obtain Here, by τ δ=1 (x) we mean the r.h.s of Eq. (22). Notice that contrary to the δ = 0 scenario, the middle parts a ∼ L and L − a ∼ L, are suppressed as 1/L 2 , so their contributions are negligible as L → ∞. The soft-modes at the edges a L and L − a L now start to play more important role because the corresponding overlaps are O(1). The prefactor in front of the Gamma functions simplifies to one and the whole series reads In order to compute this sum in the L → ∞ limit we expand it at the edges and then perform the summation of the simplified expression extending the upper limit to infinity. Namely, the asymptotics leads to This way we restore the correct result even in the different formulation of the form-factor series.

Negative winding number δ < 0
Let us consider δ = 1 − n for positive integers n ∈ Z > . The maximal number of solutions of Eq. (6) is = L + δ. We choose all of them to comprise our set k The set q a1,...an is obtained from the complete set q in Eq. (19) by the omission of the "particle" (creating a "hole") at positions q ai q a1,...
The total difference of momenta for such a state reads The corresponding overlap is analyzed thoroughly in Appendix C.4 for a i ∼ L, L − a i ∼ L. It gives the following contribution to the tau function (5) where and In order to evaluate the tau function we proceed similarly to Ref. [64,65] (see also [66]). First, we notice that one can present the product of sines in (41) as Vandermonde determinants where ε j1...jn is a completely antisymmetric tensor; the summation over repeated indices is implied. This expression is an almost factorized, so in the second step we render summation over q ai to be independent, namely This immediately allows us to write the tau function (5) in the thermodynamic limit where ν δ (q) ≡ ν(q) − (q + π)δ/(2π) has zero winding number and Y δ (x) stands for The integral has been transformed using identities such as Eq. (26) to facilitate finding asymptotic behavior at large positive x. Indeed, the exponential is an analytic function, so after the proper deformation of the integration contour it can be dropped. This way, we demonstrate that, in fact, Y δ (x) depends only on ν δ (x), namely Let us emphasize that Eq. (46) gives the exact answer for the Fredholm determinants (11). The asymptotic behavior for large x will give also asymptotics for the generalized sine-kernel determinants (18). Similar to the treatment of the δ = 0, the asymptotic expansion of Y δ (x) is connected with analytic properties of ν(q). Let us assume that the first n leading terms are given by Then the leading order of the determinant reads det 1≤j,k≤n For n = 1 we reproduce the results from the previous subsection.

Positive winding numbers δ > 1
For δ > 1, similar to δ = 1, we can keep the maximal available number of k i ∈ k, so the r.h.s sum in Eq. (5) consists of one term, which is of order O(1/L). This means that the corresponding Fredholm determinants in Eq. (11) vanish identically. The reason for this can already be seen before going into the thermodynamic limit. Namely, first we notice that the matrix A ij in Eq. (126) can be considered on the full set of momenta k = k 1 , . . . k L+δ , which will not change the determinant's limiting value But since A ij has the form of Eq. (124) which can be schematically written as for some functions ϕ and φ. This means that the rank of this matrix is maximally L and addition of the rank one matrix δV can increase the rank to at most L + 1. Therefore, for δ > 1 The corresponding determinants with sine-kernels are not zero i.e. det(1 +Ŝ ν ) = 0. In this case we see that even though the difference betweenV andŜ ν is exponentially small, it cannot be neglected, contrary to the cases for δ ≤ 1. To address this issue we modify the definition of the tau function by considering summation over k in Eq. (5) instead of q, namely where the overlaps keep their form (9) but with the modified relation between ν(q) and g(q), namely In the thermodynamic limit this sum transforms into Fredholm determinants (see Appendix A) For large x > 0, we notice that Γ is exponentially suppressed and E − (q) ≈ ie −ixq , so τ − (x) transforms into a generlized sine kernel Fredholm determinant Eq. (18) up to terms in exponentially small x. The corresponding asymptotic can be obtained in a way similar to δ < 0, however instead of summation over "holes" q a we will have summation over extra particles k a . We demonstrate how it works for δ = 2. In this case the set q consists of L elements and the set k of L+1 elements, which we parameterize by the omission one of the L+2 momenta from the all possible solutions of Eq. (6). Namely, The relative momentum of this state in the thermodynamic limit reads as The corresponding overlaps are given in Appendix C.2. Taking the thermodynamic limit we obtain the following presentation suited for the asymptotic analysis when x → +∞ For δ > 2 one can obtain similar determinant representation as for τ − (x) as in Eq. (46). Even though we have constructed τ − (x) to address positive indices δ > 1 it is possible to describe δ < 1 by the previous choice of g(k) Eq. (15) and considering x < 0.

Quantum XY spin chain and its correlation functions
In this section we consider an application of the general results obtained in the previous sections to the derivation of large distance asymptotics of thermal spin-spin correlation functions of the quantum XY spin chain. The quantum XY spin chain in a transverse field is defined by the Hamiltonian [61,67] where a periodic boundary condition for the spin operators is assumed σ α L+1 = σ α 1 , γ is an anisotropy parameter, and h is the strength of the magnetic field. The Hamiltonian H XY of XY model can be considered as an anisotropic deformation of the Hamiltonian H XX of XX model, corresponding to γ = 0. The Hamiltonian H XX can be diagonalized in two steps: Jordan-Wigner transformation to fermionic operators and a Fourier transform to momenta representation. To diagonalize the Hamiltonian H XY an additional Bogoliubov transformation is needed [50,61,67] specified by the angle θ(p): Here E(p) stands for the spectrum of the effective Dirac fermions A p , and the Hamiltonian H XY reduces to the free-fermionic one, namely We skip the details of the fermions boundary conditions as they are not important in the thermodynamic limit. We focus on the following spin-spin correlation function at finite temperature .
It is the most interesting two point correlation function as the others are either trivial in the thermodynamic limit: σ x m+1 σ y 1 T = 0, can be expressed in terms of elementary functions as σ z m+1 σ z 1 T (see Ref. [58]), or related to G(m) after the change γ → −γ as σ y m+1 σ y 1 T . We follow Ref. [50] to present G(m) in the thermodynamic limit as Fredholm determinants (m > 0): where the operatorsŴ , δW are integral operators on L 2 ([−π, π]) with the kernels given by In this form we immediately observe that G(m) can be identified with τ S (m) defined in Eq. (18) with the appropriate choice of ν(q), which can be read off from the prefactor in front of the sine-kernel This way, to find the large m asymptotics we approximate G(m) by τ (m) from Eq. (11) and use results for form factor series obtained in the previous sections. The analysis depends on the winding number δ = ν(π)−ν(−π), which can be read off from the following form of the phase shift The winding number is governed by the Bogoliubov angle θ(π) − θ(−π) = 2πδ. The possible values of δ are δ = 0, ±1 depending on the anisotropy parameter γ and the magnetic field h (see Fig. 1 for the typical behaviour of the Bogoliubov angle and the phase diagram). Notice also that Eq. (72) implies that the integral entering the asymptotic formulas can be presented as In Fig. 2 we plot exact values for the correlation function G(m) (red dots) and compare them with the asymptotic formulas written explicitly below (blue curves). We see that large m asymptotics gives reasonable approximation even for m ∼ 1. In fact, to get any visual discrepancy we had to consider large negative anisotropies in the ferromagnetic phase (δ = −1 and γ < −1). It turns out that in this case the asymptotic formulas for non-integer m acquire nonzero imaginary part, which is discarded in the plot. For integer points the imaginary part is equal to zero. Below we analyze each case separately and present analytical formulas for the asymptotics. These expressions turn out to be in accordance with the results of Ref. [58] but have a more compact form.

Paramagnetic phase h > 1 (δ = 0)
We start our consideration with relatively large magnetic field h > 1. For zero temperature such values of h correspond to the paramangetic phase, while for finite temperature the corresponding ν(q) has zero winding number δ = 0. This way, we use formula (33) to find asymptotic behavior of the correlation function G(m) at large m, namely where T 0 (m) and S(z) are given by (28) and (31), respectively. The point z 1 is the position of the pole of e −2πiν(k) outside unit circle with minimal absolute value. To find z 1 we factorize The exponent of the angle θ(k) of Bogoliubov transformation can also be presented in a factorized form, which leads to It useful to present Q(z) coth β 2 Q(z) as an infinite product Notice that in such a form the branch cut singularities disappear manifestly. Moreover, the analysis of the poles of e −2πiν(k) is now a straightforward task, from which we conclude that the smallest (by the absolute value) pole outside the unit circle is z 1 = y + for all non-zero temperatures. Therefore using Eq. (77) and Eq. (78) we obtain Finally, taking into account (73) for δ = 0, the asymptotics reads where The sign of the magnetic field h is irrelevant since it can be flipped by the conjugation of the Hamiltonian with σ x acting in each site. Therefore below we consider 0 < h < 1.

Ferromagnetic phase
In the ferromagnetic phase h < 1 with positive anisotropy γ > 0 we use the asymptotics (22) and the integral (73) for δ = 1 to obtain with For particular values of the parameters we plot exact correlation function G(m) and its asymptotics (83) in the left panel of Fig. 2. Note that even though formulas for the correlation length in different parameter regions Eq. (81) and Eq. (84) look different, the transition h < 1 and h > 1 is analytic in h. The same is true for prefactors A given by Eq. (82) and Eq. (85) (see the corresponding plots in Fig. 3). This reflects the fact that at finite temperature in one dimensional systems with short-range interactions phase transitions are absent and the physical observables are smooth functions of system parameters. This observation was used in Ref. [68] to obtain correct expressions for the correlation length and prefactor for the Ising model in the scaling limit.

Ferromagnetic phase
In this region of parameters, the correlation function G(m) is given by Eq. (46) for n = 2. We will need the large m asymptotics of Y −1 (m), which for h 2 + γ 2 = 1 is given by where, as seen from Eq. (76), Therefore Eq. (50) becomes Finally, the large distance asymptotic for G(m) following from (46) is where In the case when h 2 + γ 2 = 1 we have x + = y + and the derivation is changed slightly (in particular, Y −1 (m) ≈ (B + Cm)e −m log x+ ) however the final formula for the asymptotic of G(m) is the same. Notice that for noninteger values of m the right hand side of Eq. (91) becomes a complex valued function. We plot the typical behaviour of G(m) and the real part of its asymptotics in (91) in the right panel of Fig. 2.

Relation to Toeplitz determinants
The traditional approach to the correlation functions in the XY spin chain is in presenting them via Toeplitz determinants [58,61]. Asymptotic analysis of these structures can be performed by means of the Szegő theorem [69,70] and its generalization 3 by Hartwig and Fisher [59]. Let us comment on how similar structures can appear within our effective form factors approach. In addition to tau functions (5) and (54) that contained different number of "particles" in bra-and ket-states, we define where the quasi momenta q are solutions of e iqL = 1, while p are solutions of the following equation Here for convenience we have chosen a different notation for the phase shift. We focus on the case of nonpositive winding numbers for this function i.e. ω(π) − ω(−π) ≤ 0. The corresponding form factors read The summation in Eq. (94) can be performed using techniques developed in Appendix (A), which together with the identification e −gω(p) = e −2πiω(p) − 1 (97) leads to the Fredholm determinant expression of τ 0 whereŜ Notice that definitions of the kernels ofV andŜ ν differ from their analogues introduced in Sec. 2 by the conjugation with diagonal matrices, which does not change the value of the determinant. Comparing overlaps (9) and (96) (see Appendix C.5) we conclude that imposing the following relation between ν(q) and ω(q) we obtain exact equality for the tau functions, namely det 1 +V ν + δV ν − det 1 +V ν = det 1 +V ν1 .
Here the finite rank contribution is modified due the conjugation with the diagonal matrices Similar relations can be obtained between τ − (x) and τ 0 (x) for δ > 1. For large positive x, functions r ω (x) are exponentially small, so Eq. (104) holds for the generalized sine-kernelsŜ ν .
In fact we can easily demonstrate that this relation is true for any positive integer x. To do so we will clarify the relation between Fredholm and Toeplitz determinants (cf. [60,71]). It is convenient to slight deform of the kernel by the set of functions a 0 (p), a 1 (p), ..., a x−1 (p) S a ν (p, q) = e 2πiν(p) − 1 2π x−1 n=0 a n (p)e in(q−p) .
For a i (q) = 1 one can easily see that we recover the kernelŜ ν up to conjugation with diagonal matrices, which does not affect the value of the determinant det 1 +Ŝ ν = det 1 + S a a0=a1=...ax−1=1 .
Furthermore, we can treat S a as a product of two rectangular matrices Then using the fact that det (1 + AB) = det (1 + BA) we obtain a relation between the Fredholm determinant and determinant of matrix x by x, namely For a n = 1 the matrix T nm transforms into the Toeplitz one, namely In order to account for the finite rank we notice that because rank-one contributions are at most linear in the determinant expansion, we can present To account for the finite α one must choose a 0 (q) = 1 − αe −ixq and a n (q) = 1 for n ≥ 1, therefore Since we are looking only to the terms linear in α, we can leave only terms that are proportional to α in the first row. Moreover we can replace this row with the last one. This way we obtain wherec Here we see the shift ν(q) → ν 1 (q) as predicted from the finite size scaling of the form factors in Eq. (101). This shift together with Eq. (109) completes the proof of Eq. (104).
Let us also comment on how results of Sec. (3.3) reproduce Hartwig and Fisher asymptotic behaviour (Theorem 4 in Ref. [59]). As ν δ (q) has zero winding number we can expand it as Then the integral in the exponential in Eq. (48) can be evaluated as In this derivation we used that |e i(q+i0) | < 1 and expanded the denominator as a geometric series. Substituting this result back into Eq. (48) we immediately see that Y δ (x) = l x , where the Fourier modes l m are defined through the relation Finally, expressing double integral in the asymptotic expression Eq. (46) we obtain the statement of Theorem 4 in Ref. [59].

Summary and Outlook
In this work we have introduced the form factors (overlaps) to simulate the static correlation functions for the states with finite entropy. The state was determined by the phase shift function ν(q). For the traditional approaches dealing with the finite entropy states is notoriously difficult but for our approach it is rather advantageous situation, since almost all available quantum numbers are occupied which tremendously simplifies the computation of form factor series. This allows us, in particular, to re-derive known asymptotics for the static two point correlators in the XY spin chain and present them in a more compact form. We hope that the simplicity of this approach will make it possible to obtain the full asymptotic expansion at large distances. Apart from the thermal state we can apply our approach to the states resulting from the long time evolution after a quench [72][73][74][75], to models of 1D anyons [76][77][78][79][80], or mobile impurity models [62,81,82]. This can be done by the appropriate modification of the phase shift function. We will discuss it elsewhere.
It is interesting to note that ν(q) is apparently connected with the auxiliary functions that appears in the Quantum Transfer Matrix (QTM) approach and specifies the Bethe roots for QTM [38][39][40]83]. It would be interesting to completely clarify connection between these two approaches.
The correlation functions at zero temperature (entropy) can be formally accounted by the jump discontinuities in ν(q), which can also be treated by the form factor summation developed for the critical models [11,13]. In this case the role of the lattice is not essential and the exponential asymptotic behaviour is expected to be replaced by a power-law, which can be obtained from the proper modification of the generalized sine-kernels (see section 9 in Ref. [84]). To address dynamical correlation functions we must modify appropriately the form factors and the spectral factor e −i kix → e −i (kix− (ki)t) . The detailed constructions and extraction of the asymptotic behavior is another intriguing direction for future research. However we can already anticipate that for the space-like region, i.e. when the saddle point of the expression kx − (k)t is outside the Brillouin zone, the asymptotic analysis remains largely unchanged, which can be immediately seen in the asymptotics of Ref. [57]. For the time-like region the main problem will be that a suitable ν(q) might have a jump discontinuity which leads to additional power-law behavior (c.f. Ref. [85]). Finally, there will be extra 1/ √ t terms connected to the saddle point contributions indicated by the non-linear Luttinger theory [14][15][16].

A Summation of form factors and determinant formula
In this appendix, we derive formula (11) presenting tau function in the thermodynamic limit as a difference of two Fredholm determinants.
We consider solutions in the large L limit and choose k to fill a Fermi Sea, namely where ν i ≡ ν(k i ). For simplicity, we choose N to be even. First, we identically rewrite the overlap as (note, det D = detD) Then using standard linear algebra manipulations we rewrite the static tau function as where summation over q is happening over the whole Brillouin zone For i = j we present This sum can be rewritten as a contour integral and evaluated at large L, namely, choosing contour γ running around q i and avoiding any other singularities of the integrand we obtain Further, we deform the contour into the rectangle that encapsulates interval [−π, π]. The vertical parts of this rectangle cancel and we are left with two lines above and below the real axis along with the contribution from the pole at q = k i Here we assume that the imaginary shift i0 is chosen to be larger then Im k i = O(1/L). In this form we immediately see that, in the limit L → ∞, the values of c(k) at points k i are equal to the values of the E(k i ) for the analytic function E(k) given by Using E(k) we can obtain values also for some vicinity of k i , which allow us to effectively "omit" solving Bethe equations (6). Performing similar computation for the diagonal components we arrive at The sum can be evaluated in the same way as in Eq. (128): dq π e −g(q)+iqx sin 2 q+i0−ki 2 . (132) Equivalently, using definition (130), we can present So recalling definition of (121) we obtain for generic i and j where for i = j the second term is understood in the L'Hopital rule sense. Similarly we obtain for the finite rank contribution In this form we are at the position to take limit L → ∞, and taking into account that k i is quantized in the units 2π/L, arrive at the Fredholm determinants (11). Similarly, we can perform summation for τ − (x) defined in Eq. (54). Instead of Eq. (122) we obtain the following τ − (x) = det(A + δA) where

Submission
SciPost Physics Here k means sum over all L + δ nonequivalent (mod 2π) solutions of Eq. (6), which can be presented as a where the contour C runs around poles of the denominator only and avoids and singularities of f (k). Then the derivation goes along the lines as for τ (x). Namely, for i = j we present where now instead of Eq. (127) In the thermodynamic limit this function can be replaced by E(q), which does not depend on the system size For positive x it is much more convenient to rewrite this function as Now if we relate e −g(q) = e 2πiν(q) − 1, this function transform into For large positive x the integral can be neglected. For diagonal components we obtain Function Γ can be written as It is also exponentially suppressed for x → +∞. The finite rank contribution is easily evaluated taking into account that After all these transformations one readily obtains the result Eq. (56) in the thermodynamic limit.

B Lemmas about products
In this appendix, we study products that appear in the overlaps. In this section we assume that ν(q) is a smooth function on the segment [−π, π] and assign its values in specific points as ν j , namely First we consider constant function ν(q) = ν = const.
Lemma B.1. The following product formula is valid In the limit L → ∞ this product simplifies to The denominator is equal to Proof. We can rewrite identically the left hand side as Taking into account that we obtain Further, we proceed with the generic function ν(q).
Lemma B.2. For an integer 0 ≤ A ≤ L − 1, the following asymptotic approximation in the limit L → ∞ is valid where Proof. First, we introduce the modified product . (160) Due to this modification, it is enough to expand logB A,L [ν(q)] up to the linear terms in ν j since higher orders will be of order O(1/L), namely Taking into account (150) we transform the sum into an integral The rest of the product can be evaluated in a similar manner. First, we identically transform The logarithm of the remaining product can be expanded only up to linear in ν terms to capture finite terms in L → ∞ limit, namely Combining this result with (162) and using Stirling's formula we obtain the desired result (157).

Remark 1.
For A = L − 1, using Stirling's approximation for Gamma functions we obtain . Remark 2. For A ∼ L and L − A ∼ L the prefactor can be simplified as The next lemma is a simple corollary of the previous one.
Lemma B.3. The following asymptotic expression is valid as L → ∞ Z a ≡ sin 2 πν a L L j =a where Proof. First, we identically present this product as Then using Lemma (B.2) and Stirling's formula we obtain Changing variables we obtain the desired statement.
Further, we proceed with double products.
Lemma B.4. For δ ≥ 0 the following asymptotic expansion is valid in the limit L → ∞ where the L independent prefactor A reads with F (π) is defined in Eq. (168) and G(x) stands for Barnes G-function defined by the functional relation G(x + 1) = Γ(x)G(x). Notice that function ν(q) − δq/2π has zero winding number so the integrals in the exponential are well defined.
Proof. To find the thermodynamic limit of Z we rewrite it as Z = Y 1 Y 2 e R δ with The factors are designed in such a way that terms O(ν n ) for n > 2 do not contribute in L → ∞ case. In particular, we used that So keeping only quadratic terms we obtain and taking L → ∞ Similarly The first part of the product in Y 2 (by grouping terms with the same i − j) can be presented as We consider an additional expression Differentiating it by δ we obtain For large L we can approximate Solving Eq. (183) with initial condition log W (δ = 0) = 0 we obtain Finally, let us find L → ∞ expression for R δ defined in Eq. (176). First we identically transform it into From the form of Eq. (186) one can conclude that as L → ∞ the non-vanishing contributions to the sum will come from indices i = O(L). Therefore, using Stirling's formula we can present S i as Therefore R δ ≈ δR with So far we have proved that with Further, we can use where G(x) is Barnes G-function. The final answer is obtained by tedious but straightforward manipulations with integrals.
In the next lemma we address a similar double product for negative winding numbers δ < 0.
Lemma B.5. Let us define = L + δ, with δ < 0, then the following asymptotic behavior is valid as L → ∞ (here we still assume that |δ| L) Proof. We present this product as a ratio Z = Z 1 /Z 2 with We can identically transform Z 1 as where [ν δ ] i = ν i (1 + δ/L) − δi/L. In thermodynamic limit this expression correspond to the following function This function has zero winding number, so applying the previous lemma, we obtain Similarly, we can evaluate Z 2 . We present it as This corresponds to the positive phase shift ν(q) = −δq/(2π) = |δ|q/(2π), and allows us to use previous lemma once again and obtain Here we used that F (π) = −|δ| log(2π) for ν(q) = |δ|q/(2π) (see Eq. (168)). Finally, Eqs. (197) and (199) immediately lead to the statement of the lemma.

C Orthogonality catastrophe on the lattice
Here using results from Appendix B we evaluate the overlaps in Eq. (9).

C.1 Winding number δ = 1
For δ = 1 there exist L + 1 solutions of Eq. (6) We use all of them in Eq. (9) and set q = {q 1 , . . . q L } with To evaluate Eq. (9) in thermodynamic limit L → ∞, we first we transform identically the determinant (10) as We analyze this expression term by term. The last product can be written down using Eqs. (19) and (20) where ν + = ν L+1 and δ = ν + − ν 1 = ν(π) − ν(−π) = 1 and Notice that Where A is defined in Eq. (173). The rest of the product in Eq. (9) can be evaluated for generic δ where in the last part we have used the relation between g(q) and ν(q) Eq. (15) and assumed that ν(k) has a non-vanishing imaginary part. Combining all factors together in Eq. (9) we obtain C.2 Winding number δ = 2 For δ > 1 computation of the overlaps goes in the similar manner as in the previous section. Namely, first we consider overlap with the setk = k 1 , . . . k L+1 with k j defined in Eq. (200). There instead of Eq. (204) we will have L j=1 with F (π) defined in Eq. (168). Further, Eq. (207) we will replaced accordingly to Lemma (B.4) Taking into account Eqs. (208) and (209), for the corresponding function g(k) (see Eq. (55) ) we find the thermodynamic form for the overlap The overlaps for other sets k can be obtained from this one. We further focus on δ = 2, in this case there are exacly L + 2 sets k parametrized by the omission of one of the solutions of Eq. (6), namely k (a) = {k 1 , . . . , k a−1 , k a+1 , . . . , k L+2 }, a = 1, 2, . . . , L + 2.
With this notations k (L+2) =k. Now let us consider ratio of the excited overlap Using Lemma (B.2) for a ∼ L, L − a ∼ L we obtain Combining this result with Eq. (213) for δ = 2, the overlap can be written as C.3 Winding number δ = 0 Let us study thermodynamic limit of the overlap (9) in the case N = L − 1, which is especially useful for δ = 0. Below, however, for the sake of generality, we will keep δ ≥ 0. Our goal is to evaluate Z a defined via We use notations (20) and (19) to label the momenta and Eq.
We can present it identically as where F (q a ) is given by Eq.