Exact thermal properties of free-fermionic spin chains

An exact description of integrable spin chains at finite temperature is provided using an elementary algebraic approach in the complete Hilbert space of the system. We focus on spin chain models that admit a description in terms of free fermions, including paradigmatic examples such as the one-dimensional transverse-field quantum Ising and XY models. The exact partition function is derived and compared with the ubiquitous approximation in which only the positive parity sector of the energy spectrum is considered. Errors stemming from this approximation are identified in the neighborhood of the critical point at low temperatures. We further provide the full counting statistics of a wide class of observables at thermal equilibrium and characterize in detail the thermal distribution of the kink number and transverse magnetization in the transverse-field quantum Ising chain.


Introduction
Quantum many-body spin systems that are exactly solvable and exhibit a quantum phase transition have been key to advance our understanding of critical phenomena in the quantum domain. Among them, the one-dimensional XY model and the closely-related transverse-field quantum Ising model (TFQIM) occupy a unique status, and are paradigmatic test-beds of quantum critical behavior [1][2][3][4]. They belongs to a family of models that admit an exact diagonalization by a combination of Jordan-Wigner and Fourier transformations, yielding a formulation of the system in terms of free fermions [5,6]. These family of quasi-free fermion models include as well the Kitaev spin model in one dimension and on a honeycomb lattice [7], among other examples [1,2,4].
Quasi-free fermion models have indeed been instrumental in exploring both equilibrium and nonequilibrium properties. At equilibrium, the study of the ground-state critical behavior was shown to be of relevance to the characterization of the system at finite temperature [8][9][10]. Out of equilibrium, these models have been used to explore the dynamics following a sudden quench (e.g., of the magnetic field). The study of finite-time quenches was key to establish the validity of the universal Kibble-Zurek mechanism in the quantum domain, and confirm the power-law scaling of the number of kinks by driving the ground-state of a paramagnet across the phase transition [11,12], as reported in a variety of experiments [13][14][15][16]. These results have also been extended to nonlinear quenches [17,18] and inhomogeneous systems [19][20][21][22][23], while their breakdown has been characterized in open systems [24][25][26]. More recently, it has been shown that signatures of universality are present in the full kink-number distribution and that all cumulants scale as a universal power-law of the quench time [16,[27][28][29][30][31]. The universal dynamics of defect formation is not always desirable, and a variety of works have been devoted to circumvent it using diverse control protocols [32][33][34][35][36][37][38][39][40][41][42], beyond the use of nonlinear quenches and inhomogeneous driving. In addition, quasi-free fermion models have been discussed in the context of quantum thermodynamics, as a test-bed to explore work statistics and fluctuation theorems [43][44][45][46] and as a working substance in a quantum thermodynamic cycle [47].
Quasi-free fermion models provided an effective description of a variety of condensed-matter systems, where they can be realized with high accuracy in [48]. They are further amenable to quantum simulation with trapped ions [49][50][51][52][53], ultracold gases in optical lattices [54] and superconducting qubits [55]. Digital quantum simulation provides yet another avenue for their study in the laboratory [15,[56][57][58].
In many applications, it is generally desirable to consider a thermal state and analyze the finitetemperature behavior. For a given observable, full information about the eigenvalue distribution and its cumulants can be extracted from the characteristic function. An ubiquitous approximation in such description exploits the parity symmetry of the TFQIM and XY modes, focusing on the positiveparity subspace, while disregarding the rest of the spectrum [1-4, 44, 59-62]. We refer to it as the positive-parity approximation or PPA for short. The PPA is considered to be accurate in the thermodynamic limit [63], invoked in many works [6,64]. However, even in the thermodynamic limit, an exact treatment requires taking into account parity properly and at finite temperature both subspaces are populated. Katsura derived the exact partition function for a finite-size spin chain in 1962 [63]. Kapitonov and Il'inskii provided an alternative derivation of the closed form expression of the exact partition function using functional integrals over Grassmann variables [65]. More recently, Fei and Quan [45] used group theory methods to calculate the exact partition function and quantum work distribution.
In this manuscript, we first elaborate on these results providing an elementary derivation of the exact partition function based on the structure of the Hilbert space. Using this approach, we next provide exact expressions for the eigenvalue distribution (full counting statistics) of a wide class of observables at thermal equilibrium. We present step-by-step worked examples deriving the exact moment generating function for important observables: the kink number and transverse magnetization. In addition, we analyze finite-size effects and illustrate discrepancies between results obtained using the PPA for the partition function and the exact partition function for small systems spins. These discrepancies are of direct relevance to typical system sizes in current experimental realizations of spin systems [66,67]. For convenience of the reader interested in using the final results of a calculation, the corresponding explicit formulas are summarized in boxes that are self-contained and make little or no reference to the rest of the manuscript.

Full Diagonalization of Spin-1 XY Model
We consider the anisotropic one-dimensional XY Hamiltonian for spins 1/2 in a transverse magnetic field g. The Hamiltonian reads: Here, J parameterizes the ferromagnetic (J > 0) or antiferromagnetic (J < 0) exchange interaction between nearest neighbors; we set the energy scale by taking J = 1. The dimensionless anisotropic parameter in the XY plane is given by γ > 0 and L is the number of sites in the chain. For γ = 1, the Hamiltonian (1) corresponds to the Ising model in a transverse magnetic field, which possesses a Z 2 symmetry. The limit γ = 0 describes the isotropic XY model. For the anisotropic case 0 < γ ≤ 1 the model belongs to the Ising universality class, and its phase diagram is determined by the ratio ν = g/J. When ν > 1, the magnetic field dominates over the nearest-neighbor coupling, polarizing the spins along the z direction. This corresponds to a paramagnetic state, with zero magnetization in the xy plane. By contrast, in the regime 0 ≤ ν < 1 the ground state of the system corresponds to a ferromagnetic configuration with polarization along the xy plane. These phases are separated by a quantum phase transition (QPT) at the critical point ν = 1. Finally, for the isotropic case γ = 0, a QPT is observed between gapless (ν < 1) and ferromagnetic (ν > 1) phases. The operatorsX n ,Ŷ n , andẐ n are matrices of order 2 L defined by the relationŝ Here,σ α n denotes the Pauli operator at site n along the axis α = x, y, z,Î n is the identity matrix of order 2 at the site n, and periodic boundary conditions are assumed,σ α L+1 =σ α 1 . A standard way to diagonalize the Hamiltonian in Eq. (1) relies on introducing a new set of Fermionic operators given byσ These expressions represent the well-known Jordan-Wigner transformation [68]. Here,ĉ n andĉ † n are ladder Fermionic operators at site n, which satisfy anti-commutation relations ĉ i ,ĉ † j = δ i, j and ĉ i ,ĉ j = ĉ † i ,ĉ † j = 0. This is in contrast to the Pauli matrices, which satisfy commutation relations σ † n ,σ − m = δ n,mσ z n and σ z n ,σ ± m = ±2δ n,mσ ± n withσ ± n =σ x n ± iσ y n . With periodic boundary conditions in the spin representation, the Fermionic operatorsĉ n andĉ † n satisfy nontrivial boundary conditionŝ Here, the parity operatorΠ is given by (−1)N = exp iπN and has eigenvalues ±1. The parity operator anticommmutes with the creationĉ † n and annihilationĉ n Fermionic operators, (−1)N ,ĉ † n = (−1)N ,ĉ n = 0, and therefore, it commutes with any operator bilinear inĉ † n andĉ n . The Hamiltonian given by Eq. (5) does not conserve the number of Fermionic excitations. However, it is well-known that the TFQIM has a global Z 2 symmetry and, thus, the parity operatorΠ commutes with the Hamiltonian. As a result, the total Hilbert space is split into the direct sum of two 2 L−1 dimensional subspaces of positive (+1) and negative (−1) parity. Using the projectorsΠ ± , the Hamiltonian in Eq. (5) is represented in the form with the reduced HamiltoniansĤ ± being given bŷ A subtle difference betweenĤ + andĤ − is found in the boundary conditions for the Fermion operators. H + obeys antiperiodic boundary conditions (ĉ L+1 = −ĉ 1 andĉ † L+1 = −ĉ † 1 ) whileĤ − satisfies periodic boundary conditions (ĉ L+1 =ĉ 1 andĉ † L+1 =ĉ † 1 ). The Hamiltonian given by Eq. (8) is quadratic in the Fermionic operators and is thus exactly diagonalizable using Fourier and Bogoliubov transformations [64,[69][70][71]. We expand the operatorĉ n via a Fourier transformation in momentum space, The wavevector k takes values in the positive K + and negative K − parity sectors We emphasize that Eqs. (10) and (11) are valid for an even and odd number of particles in the chain.
In the following analysis, we consider even L. In this way, the modes k = 0 and k = π are included in the negative parity sector. For even L, we can rewrite conveniently the momentum values as By direct substitution of Eq. (9) into Eq. (8), the reduced HamiltoniansĤ + andĤ − are expressed in terms ofĉ k andĉ † k asĤ , whereĤ We next make use of a Bogoliubov transformation, and define a new set of fermion operatorsγ k and γ † k given byγ where the real numbers u k and v k satisfy u k = u −k , v k = −v −k and |u k | 2 + |v k | 2 = 1. The canonical anticommutation relations for the operatorsĉ k andĉ † k imply that the same relations are also satisfied byγ k andγ † k , that is, γ k ,γ † k = δ k,k , and γ † k ,γ † k = {γ k ,γ k } = 0. By direct substitution of the Bogoliubov transformations into Eq. (13), after a some algebra, we obtain The terms proportional to γ † k γ † −k and γ k γ −k should vanish for the Hamiltonian to acquire a diagonal form. Writing u k = cos (ϑ k /2) and v k = sin (ϑ k /2), the Bogoliubov angles satisfy .
For numerical simulations, the last condition can be rewritten as γ sin (k) u 2 k − v 2 k +2 (cos (k) − g) u k v k = 0. Finally, the Hamiltonian (13) can be rewritten as a sum of noninteracting termŝ withn k =γ † kγ k denoting the fermion number operator and k (g, γ) = 2 (g − cos k) 2 + γ 2 sin 2 k being the quasiparticle energy of mode k 0, π per particle.

Mathematical tools for the complete Hilbert space
To simplify the presentation, we focus on the positive-parity subspace in this subsection. However, the methods presented are applicable in the negative-parity sector too. In order to keep the notation clear, we use the following conventions: • Hilbert spaces are denoted by letters in blackboard bold style, for example k .
• Operators are denoted by letters with a hat, such asÔ k andĥ k i .
• Operations on tensor products of Hilbert spaces are denoted with calligraphic letters P and N.
To begin with, we note that the positive-parity Hilbert subspace + can be written as the tensor product of subspaces corresponding to each pair of momenta (k and −k) Each subspace k is the linear span of the vacuum and states involving one and two Fermionic excitations with a given momentum Here, |0 k is the vector annihilated by bothĉ k andĉ −k . Each of the subspaces can be divided into the sectors with even (p) k and odd (n) k number of excitations Note that the dimension of the right hand side of equation (19) is equal to 4 L/2 = 2 L , as there are L/2 positive momenta and each corresponding subspace is four-dimensional. However, there is an additional condition in the positive-parity subspace: the parity operatorΠ has eigenvalue +1. Thus, the subspace is only spanned by vectors associated with an even number of quasiparticles. We denote this subspace by P Similarly, we define the subspace spanned by odd number of quasi-particle excitations and denote it by N = N k∈k + k . It is easy to see that both spaces P k∈k + k and N k∈k + k have dimension 2 L−1 and satisfy For the positive-parity subspace only P is relevant; vectors in N have no physical meaning for the system described by the HamiltonianĤ + . However, the spaces P and N (defined for proper momenta) exchange their roles forĤ − ; see Eq. (18). These considerations suggest that to obtain correct results in the positive-parity subspace, it is sufficient to redefine the tensor product to take into account only vectors from P. This can be done for states and observables. Before dealing with observables, we introduce an alternative recursive definition of the spaces P and N, equivalent to Eq. (22). We shall make use of it in deriving the exact partition function and characteristic functions of observables. We start by defining the subspaces for one momentum, see Eq. (21), Next, we specify how to construct spaces P and N when a mode with momentum k n+1 is added: The intuitive meaning of these equations is that in order to obtain an even number of excitations one has to add an even number of excitations to an even number, or an odd number of excitations to an odd number.
We can extend these definitions for operators and density matrices. We assume that operatorsÔ k act independently on each subspace k and eachÔ k can be written as a sum of an even partÔ (p) k and an odd partÔ (n) k asÔ The operatorsÔ (p) k andÔ (n) k act on the total space k , but have a 2 × 2 zero block 0 2 in the respective subspace. The proper restrictions of the tensor product of operatorsÔ k can be defined in a similar way as in Eqs. (24) and (25) k 1 , and are given by Example 2.1: Even and odd parity parts of the Hamiltonian ForĤ k given by Eq. (14), note that for a each mode k n the Hamiltonian can be rewritten aŝ Here,ĥ (n) k n is 4 × 4 zero matrix (with no odd part), andĥ (p) k n =ĥ k n .
As the odd part of Hamiltonian is zero, the description using ordinary tensor products instead of over P is valid for pure states. However, the canonical thermal Gibbs state has a non-vanishing odd-parity contribution:

Example 2.2: Even and odd-parity contributions to the exact Gibbs state
Consider the part of the thermal Gibbs state corresponding to momentum k: Using the expression forĥ k in the the basis Therefore, the even and odd parts read: Using the fact thatĥ k has eigenvalues ± k , we have: Next, we state three propositions helpful in calculating the complete and exact expression of the partition function and the full counting statistics of observables:

Proposition 2.3: Identities for product of operators
Consider two operatorsÔ k andR k acting independently on each subspace k . Then, the following identities are true for operator multiplication The following proposition is useful in calculations involving Gibbs states and time-evolutions:

Proposition 2.4: Identities for exponentials of operators
For every set of operators O k acting on the subspace k , the following identities for exponents of operators hold: Lastly, the use of traces turns out to be essential to determine expectation values of observables, and, more generally, their full counting statistics:

Proposition 2.5: Trace identities
Consider operatorsÔ k that act independently on each subspace k . Then, the traces of the restricted tensor products can be expressed as follows, We present a proof ot Eq. (34) in the Appendix A.
Negative-parity subspace. In the negative-parity subspace, all formulas derived for the positiveparity subspace remain valid. In particular, for all momenta k 0, π expressions from examples 2.1, 2.2 apply. The only difference is that one has to treat carefully the parts of the Hilbert space associated with momenta 0 and π. They are spanned by the following bases: As a result, matrices describing the Hamiltonian and Gibbs state are 2 × 2 instead of 4 × 4. In the following example we give formulas for the even-and odd-parity parts of the Gibbs state in modes k = 0, π: Example 2.6: Even-and odd-parity parts of the exact Gibbs state for 0, π momenta Using equation (16), the explicit form of the Gibbs state of the modes with momenta 0, π, in the bases {|0 0 ,ĉ † 0 |0 0 }, {|0 π ,ĉ † π |0 π }, are respectively given bŷ Thus, the corresponding even-and odd-parity parts read In closing this section, we point out that when L is odd, the momenta 0 and π appear in the positive-parity subspace; the general formulas (24) and (26) are always valid.

The Canonical Partition Function
The partition function is a fundamental object in statistical mechanics from which all equilibrium thermal properties of a system can be derived. It further facilitates the study of critical phenomena through the study of its zeroes in the complex plane, know as Lee-Yang zeros [72].
For its study, we consider a linear spin−1/2 chain described by Eq. (1). The system is prepared in a canonical thermal Gibbs state at finite inverse temperature β and characterized by the initial density operatorρ where Z (β, g, γ) is the canonical partition function given by In a Gibbs state, the system is in a mixture of positive-and negative-parity states and both subspaces should be taken into account. To this end, we consider the operatorρ = exp −βĤ , whereĤ is given by Eq. (1). According to the exact diagonalization (see Sec. 2), the total Hamiltonian can be mapped to a set of independent mode operators in each parity sector. For fixed even L , the operatorρ is given byρ whereρ andρ k are defined in Examples 2.2, with the sets k + and k − given in Eq. (12). For these operators the corresponding reduced partition functions are For simplicity, we calculate Z + and Z − separately, and focus on Z + first. Using the formulas from Example 2.2, one finds Making use of the first identity in (34), we obtain an expression for canonical partition function in the positive-parity sector The computation of the negative-parity part of the partition function proceeds in the same way; we use the second of the trace identities (34) and the expressions from the example 2.1 to find Using (40), the exact partition is the sum of contributions of positive and negative parity: . To sum up, one can rewrite exact partition function in closed-form.
Summary 3.1: Exact partition function for spin-1 2 XY model where In this expression the products run over all momenta, not only those with non-negative values. In general, the total partition function can be represented as the sum of four contributions, where Z ± F (β, g, γ) = k∈K ± 2 cosh (β k (g, γ) /2) and Z ± B (β, g, γ) = k∈K ± 2 sinh (β k (g, γ) /2) are the "Fermionic" and "boundary" contributions. The first term, which takes only into account Fermionic and positive-parity contribution is the only term considered in the PPA, widely used in the literature as the correct approximation in the limit N → ∞ [1,44,59,60,62,64] Summary 3.2: PPA partition function In the isotropic case with γ = 0, the exact partition function admits a more compact expression [73] but this limit lies outside the Ising universality class, our primary focus. The complete expression for the partition function 3.1 was first derived with the aid of creation and annihilation operators by Katsura [63]. An alternative approach has been reported using Grassmann variables, without a numerical characterization [65]; see as well [45].
It is thus natural to analyze the extent to which the PPA Z + F (β, g, γ) provides a valid approximation to the exact partition function. Fig. 1 shows the difference between the ratio Z + F (β) /Z (β) as a function of the inverse of temperature and the magnetic field. The error is negligible away from criticality and at high temperatures. However, prominent discrepancies between the exact partition function 48 and the ubiquitously-used PPA (49) are manifested in the neighborhood of the critical point in the regime of low-temperatures, which is often times the regime studied and of interest. Indeed, in this region errors reach sufficiently large values such that Z + F (β, g, γ) ≈ 0.5 Z (β, g, γ). One can provide a simple and intuitive explanation of the magnitude of this discrepancy by considering the structure of the spectrum. The complete spectrum consists of two disjoint "ladders" of levels, spanning the positive-parity and negative-parity subspaces. In the following analysis we denote by E α g and |g α the lowest energy level and the corresponding eigenstate in the subspace of parity α = ±. The diagonalization procedure of the Ising model yields explicit formulas for these eigenvalues. For even Figure 1: Comparison of the exact and PPA canonical partition functions. The ratio between the total partition function 3.1 and the PPA Eq. (49) is shown in the β-g plane for finite system size L = 50, 100, 5000, 10000, increasing from left to right (anisotropic parameter γ = 1). Significant differences appear close to the critical point g = g c = 1, with the magnitude of Z + F (β, g, γ) deviating by 50% from the exact partition function. The paramagnetic phase is correctly reproduced by the simplified approximation g > g c while errors in the partition function are shown in red in the ferromagnetic phase at low temperatures. number of spins [74] The corresponding eigenstates read where |vac is annihilated by allĉ k for k ∈ K + ∪ K − (including 0 and π modes). In what follows, we restrict ourselves to the TFQIM (γ = 1). In the TFQIM with even number of spins L, the true ground state always lies in the positive-parity subspace (this is not necessary true in the XY model, see [75]). The energy gap δ(g) between these two lowest energy states plays a crucial role. We recall its asymptotic behavior [74] δ(0 < g < 1) = O ∼ exp (−L/ξ(g)) , where ξ(g) denotes the correlation length. In the low temperature regime, the Gibbs state is effectively spanned by the two lowest energy states, |g + and |g − . In this truncation, the partition function and Gibbs state read This low-temperature two-level approximation relies on (51) and disregards the contribution from higher excited states, that are energetically separated from |g + and |g − . The energy gap to the next excited state can be calculated as the energy of a single-particle excitation in the positive-parity subspace, which sufficiently far from the critical point is estimated by while at the critical point, this gap behaves as In the ferromagnetic phase, the first excited state is separated from the ground state by an exponentially vanishing gap and the second excited state lies far away from both of them. Therefore, the correction from high-energy states is negligible in the low temperature limit β∆(g) 1. Similarly, in Figure 2: Ratio between the low-temperature approximation and exact partition functions as a function of the inverse temperature. The accuracy of the two-level approximation (53) is considered for different values of the transverse magnetic field g and two different system sizes. As the energy gaps δ(g) and ∆(g) in the neighbourhood of g c = 1 are comparable, a lower temperature is required to obtain a desired level of accuracy. For given β, the accuracy decreases with increasing system size.
the paramagnetic phase, the ground state is energetically separated from all the excited states. At the critical point the two lowest excited states are separated from the ground state by a comparable gap, However, for large β the error is very small. The accuracy of the the two-level approximation for different phases is shown in Fig. 2. The validity of this approximation (53) explains the magnitude of the errors between the exact and the PPA partition functions shown in Fig. 1. For g < 1, the simplified partition function takes into account only the ground state |g + and can be approximated by e −βE + g , while the complete partition function is approximately This explains the observed error of about 50% between the exact and PPA partition functions.

Full Counting Statistics in Integrable Spin Chains
The characterization of a given observable in a quantum system generally relies on the study of its expectation value. To determine it, experiments often collect a number of measurements, and build a histogram, from which the eigenvalue distribution is estimated. The full counting statistics of an observable focuses on the complete eigenvalue distribution. Its study has proved useful in a wide variety of applications and alternative methods for its measurement have been put forward [76]. A prominent example concerns the counting statistics of the number of fermions (electrons) traversing a point contact in a wire, that is described by the Levitov-Lesovik formula [77][78][79]. Distributions of other observables such as the total energy play a key role in quantum chaos [80] and the statistics of related positive-operator valued measures (POVMs, such as work) are at the core of fluctuation theorems in quantum thermodynamics [81]. In the context of spin chains, the distribution of the order parameter has long been recognized as a probe for criticality and turbulence [82][83][84][85][86][87][88][89][90]. Further, the study of the full counting statistics of quasiparticles and topological defects has been key to uncover universal dynamics of phase transitions beyond the paradigmatic Kibble-Zurek mechanism [16,[27][28][29][30]91].
The full counting statistics is characterized by the probability P (ω) to obtain the eigenvalue ω of a general operatorŴ. It is defined as the expectation value where the δ function is to be interpreted as a Kronecker or Dirac delta function, depending on whether the spectrum ofŴ is point-wise or continuous. The angular bracket denotes the quantum expectation value with respect to a general state characterized by a density matrixρ. We introduce the Fourier transform representation whereP (θ) is the characteristic function given bỹ In cases such as the kink number and the transverse magnetization, the eigenvalues are integers ω ∈ Z and the range of the integral can be restricted from −π to π. The characteristic function is also known as the moment generating function, as it allows to directly compute the mean value and higher-order moments of a given observableŴ according to Further, its logarithm is the cumulant generating function used to derive the cumulants of the distribution through the identity The first cumulant κ 1 is just the mean value, κ 2 is the variance, and κ 3 coincides with the third central moment. Cumulants are useful in characterizing fluctuations in a quantum system. For example, since the only distribution with finite κ 1 , κ 2 0 and vanishing κ m = 0 for m > 2 is the Gaussian distribution, higher cumulants quantify non-normal features of the distribution of interest, e.g., an eigenvalue distribution.
We next derive the general form of characteristic function for a wide class W of observables. This class is defined by the property that any operatorŴ ∈ W, in each parity subspace, can be written in the formŴ = kŴ k , whereŴ and the matrixŵ k has the block-diagonal form Here,ŵ (1) k andŵ (2) k are 2 × 2 are matrices for momenta different from 0, π and 1 × 1 matrices for 0, π momenta. We point out that the notation in equations (65,66) is compatible with matrix expressions from Examples 2.1, 2.2, written in the basis {|00 k , |11 k , |01 k , |10 k }. When off-diagonal blocks vanish, the operatorŴ k can be written as a quadratic form in Fermionic operators. However, there are relevant observables which have components linear in Fermionic operators. For example, the longitudinal magnetizations X i or Y i do not belong to the class W as these observables mix the subspaces with different parities. This leads to severe difficulties. Namely, one has to switch between fermionic operators defined for different sets of momenta. For example, for a given k ∈ K − one has to perform inverse Fourier transform to express c k as a linear combination of fermionic operators in space domain, and then apply Fourier transform with K + as a set of momenta. This will cause, that subspaces with different momenta will be all intertwined, contrary to the basic feature exploited in the systems with periodic boundary conditions -that one can perform calculations independently for every momentum. The treatment of such operators is thus beyond the scope of this paper. The systematic treatment of longitudinal magnetization in zero temperature for Ising model with PBC was conducted in [92], with further extensions including observables involving three fermionic operators in [89]. For other approaches, see for example [93] (open boundary conditions) or [61] (exploiting methods of field theory).
In the following we present the detailed procedure for computing characteristic functionP (θ) of a given observableŴ in the class W.
2. As in the case of the partition function, it is convenient to separate in the full characteristic function the contributions of positive and negative parity: Using Propositions 2.3 and 2.4, we aim at calculating Next, we define the matrixσ Denoting the eigenvalues ofŵ (2) k by µ k and λ k we find tr ρ k exp (iθŵ k ) =σ 11 k e −β k (g,γ) +σ 22 k e β k (g,γ) + e iθµ k + e iθλ k , tr ρ Using Proposition 2.5 we obtain 3. To determineP − (θ) it remains to compute the contributions corresponding to 0, π momenta. Denotingŵ one findsρ Therefore, the negative-parity part of the characteristic function is whereP F (θ) = e β(g−1)+iθw 1 0 + e −β(g−1)+iθw 2 0 e β(g+1)+iθw 1 π + e −β(g+1)+iθw 2 π , Note that this is not the only way to calculate the characteristic function: instead of diagonalizingρ k , one could diagonalize an observableŵ k . However, in our approach the role of the Boltzmann factor set by β k (g, γ), which is usually dominant, is clear from the formulas (74) and (77). In the following sections we apply this method to characterize the full counting statistics of two physically important observables, the number of kinks and the transverse magnetization.

Probability distribution of the number of kinks at thermal equilibrium
We next derive the full generating function for the kink-number operator, which is of fundamental importance in the study of quantum phase transitions [12,16,[27][28][29][30]. Although the relevance of this operator is most apparent in the Ising model, it is also well-defined for the general XY model. In the following, we consider the TFQIM with γ = 1 for simplicity. The explicit form of kink-number operator readsN with eigenvalues n = 0, 1, . . . , L under periodic boundary conditions. Comparing the Ising Hamiltonian Eq. (1), with γ = 1 and g = 0, with the Bogoliubov Hamiltonian (18) at γ = 1 and g = 0, the kink operator takes a simple form as the sum of the number operators of quasiparticles in each momentum [12]. Here, we generalize the kink number operator definition for all values of the magnetic field. First, we rewrite the operator (79) in the following form: By analogy with Eq. (65) and Eq. (66), we define a new set of operatorsn k ,n 0 , andn π ; taking for any mode k 0, π the basis given by {|00 k , |11 k , |01 k , |10 k }, while selecting for 0, π momenta the basis Therefore, we define the operatorŝ and thusn Note that exp iθn (1) k has the simple form Using expressions (68) and (72), one finds This yields the explicit expression of the full characteristic function of the kink-number operator.

Summary 4.1: Full characteristic function for kink number operator
The full characteristic function of the kink number operator Eq. (79) at thermal equilibrium readsP Positive part of characteristic function: Negative part of characteristic function: whereP The exact total partition function is given by Eq. (48), with the eigenenergies k (g, γ) and k=0 given by Eq. (47), and the Bogoliubov angles ϑ k satisfying Eq. (17).
By contrast, in the customary PPA, the characteristic function of the kink-number operator in the thermodynamic limit contains only the first term ofP + (θ):
In Figure 3, we characterize the full counting statistics of kinks as a function of the magnetic field and inverse temperature. By numerical integration of Eq. (60), we find the exact probability distribution function P (n) using Eq. (85). Additionally, we evaluate the PPA probability distribution function using Eq. (89). The use of the PPA partition functions is widely extended in the literature, e.g., to analyze the formation of kinks after non-equilibrium quenches [1,44,59,60,62,64]. For a large magnetic field and low temperature, the PPA works well and reproduces essentially the exact full counting statistics of kinks. By contrast, when thermal fluctuations are suppressed and the magnetic field contribution dominates, the PPA leads to pronounced discrepancies (i.e. see Fig. 3 lower-left panels). The PPA also fails to account for momentum conservation. Under periodic boundary conditions, kinks appear in pairs. In general, the PPA incorrectly predicts a non-zero probability of exciting odd number of kinks: but for large g and β as shown in 3, when P PPA (n = 2 + 1) ≈ 0.
The fact that only even number of kinks in the presence of periodic boundary conditions can be excited is intuitively clear. For a simple mathematical argument, consider the operator L n=1X nXn+1 which is 1 for even kink number and −1 for an odd number. UsingX L+1 =X 1 and X n 2 = L n=1Î n , it satisfies: The PPA characteristic function,P + (θ) andP − (θ) do not exhibit this feature.
In addition, we note that the magnitude of the exact P(n) for even n can be approximated by the coarse-grained PPA approximation, whenever the distribution is symmetric, with tails far from the origin, i.e., P (n) ≈ P PPA (n) + 1 2 [P PPA (n − 1) + P PPA (n + 1)] , as shown in Figure 4. An analysis of the cumulants of the kink-number distribution as a function of the inverse temperature is presented in Fig. 5 for various system sizes. In the paramagnetic phase (g > 1), the mean always exceeds the variance, making the kink-number distribution sub-Poissonian. This need not be the case in the ferromagnetic phase, where the distribution changes from sub-Poissonian to super-Poissonian as the temperature decreases. This behavior is shown to be robust as a function of the system size. The difference between the exact cumulant values and those derived from the PPA is systematically studied in Fig. 6 for a system size of L = 12 spins; the relative error is reduced with increasing system size. The quality of the PPA improves with increasing temperature, in the classical regime, in the ferromagnetic phase. While the dependence of the relative error as a function of the magnetic field g is not monotonic, the bigger discrepancies between the exact results and the PPA are found in the ferromagnetic phase in the low temperature regime, when the relative error can reach 100%. In the paramagnetic phase, the PPA provides an accurate description of the cumulants for different temperatures and values of the magnetic field.
To complete the characterization of the kink-number distribution we consider the limiting cases of the ground-state distribution (β → ∞) and the infinite-temperature case (β → 0) in an exact approach, without using the PPA. The first can be easily described using (83), while in the second we consider a maximally-mixed Gibbs state and apply trace formulas 2.5. For β = 0, the exact result and the PPA coincide. Figure 5: Cumulants of the kink-number distribution as a function of the inverse of temperature β. Using the exact characteristic function given by Eq. (85), the mean kink number κ 1 and the variance κ 2 are shown by red circles and blue squares, respectively. The dashed lines correspond to the numerical results using the PPA characteristic function in Eq. (89). While in the paramagnetic phase the statistics is sub-Poissonian, in the ferromagnetic phase it changes from sub-to super-Poissonian as the temperature is decreased. The magnetic field is increased from 0.0 to 2.0, varying from left to right in steps of 0.5. In the upper panels, the system size is L = 12, while in the lower ones L = 100.

Summary 4.3: Limiting cases of kink number distribution
Exact ground-state characteristic function of the kink-number distribution: Exact infinite-temperature characteristic function of the kink-number distribution: Instances of the corresponding distributions are shown in Fig. 7 for the (pure) ground-state as a function the magnetic field. For g = 0 one finds a Kronecker delta distribution centered at n = 0, with P(0) = 1 and P(n) = 0 for n > 1, as expected. As the magnetic field is cranked up, the distribution broadens and gradually shifts away from the origin, becoming approximately symmetric in the paramagnetic phase.
The right panel in Fig. 7 also shows the corresponding distribution in the infinite-temperature case, that is symmetric, centered at n = L/2 and independent of the transverse magnetic field g, as can be seen from Eq. (94). In fact, full probability distribution for infinite temperature can be found by a combinatorial argument. Working in the basis of eigenstates of σ x i in each site, the probability of obtaining n = 2l kinks is related to the number of basis vectors with 2l spin flips, where we use the fact that an even number of kinks is enforced by boundary conditions. One can choose the location of 2l kinks in the chain in 2 L 2l ways. Therefore, the full probability distribution has the form: The corresponding cumulant values read By keeping the first two cumulants and setting the rest to zero, P β→0 (n = 2l) can be approximated by a Gaussian distribution N(κ 1 , κ 2 ) with mean κ 1 = L/2 and variance κ 2 = L/4. As shown in Fig. 7 this approximation describes the envelope of the distribution with great accuracy.

Probability distribution for the transverse magnetization at thermal equilibrium
We next focus on the derivation of the explicit form of the characteristic function of the transverse magnetizationM behavior [82][83][84][85][86]88,89] and the identification of various many-body states in ultracold-atom quantum simulators [87].
In the Fourier representation, it is the sum of two different contributions: In parallel with Eq. (81), we define a new set of a single-mode operatorsm k ,m 0 , andm π , In addition, in the negative-parity sector, the matrixm k has the same form for the momenta 0, π that is given bym 0 =m π = diag (1, −1). We can easily compute exp iθm (1) k and theσ k matrix to obtain σ (11) k = cos(2θ) + i cos(ϑ k ) sin(2θ),

Summary 4.4: Full generating function of transverse magnetization
The full characteristic function for the transverse magnetization Eq. (97) at thermal equilibrium readsP Positive part of characteristic function: γ)) cos(ϑ k ) + 1) Negative part of characteristic function: withP The exact partition function is given by Eq. (48), with the eigenenergies k (g, γ) and k=0 given by Eq. (47), and the Bogoliubov angles ϑ k satisfying Eq. (17).
By contrast, in the PPA, the characteristic function of the transverse magnetization in the thermodynamic limit contains only the first term ofP + (θ):

Summary 4.5: PPA characteristic function for transverse magnetization
In the thermodynamic limit, Eq. (101) reduces tõ P PPA (θ) = 1 2Z + F (β, g, γ) k∈k + 2 (cos(2θ) cosh(β k (g, γ)) − i sin(2θ) sinh(β k (g, γ)) cos(ϑ k ) + 1) , The magnetization distribution is shown in Fig. 8 for different values of g and β for a fixed system size L = 50. The distribution P(m) vanishes for odd values of m for even L. It is naturally symmetric for g = 0 and approximately so for finite g in the high-temperature case at low magnetic fields, when it approaches a binomial distribution. The accuracy of the PPA is remarkable as a function of g and β with discrepancies being noticeable in the pure ferromagnet (g = 0) at low temperature. As the magnetic field is cranked up at constant β, the alignment of the spins is favored shifting the mean and increasing the negative skewness of the distribution in the paramagnetic phase. Figure 9 provides a systematic characterization of the first two cumulants as a function of the inverse temperature for different values of g. In contrast with the kink-number distribution, in the ferromagnetic phase the variance always exceeds the mean, and thus the magnetization distribution remains super-Poissonian. In the paramagnetic phase, at any fixed value of g the variance decreases with temperature, while the converse is true for the mean magnetization. As a result, the character of the distribution changes from super-Poissonian to sub-Poissonian as the the temperature is lowered. The behavior of P(m) is shown to be robust as a function of the system size, with discrepancies between the exact results and the PPA being restricted to the critical point. The relative error of the PPA remains below 10% as a function of g and β as shown in Fig. 10.
As in the case of kink number distribution, we close with a characterization of the magnetization distribution in the limits of infinite and vanishing inverse temperature β.  (106) Exact infinite-temperature characteristic function of transverse magnetization: The behavior of the ground-state magnetization distribution is the reverse of the kink-number distribution in the sense that it becomes approximately symmetric in the ferromagnetic phase and sharply peaked at m = L in the paramagnetic phase. Using formulas (106) and (63), one can find the first cumulants of the ground-state distribution explicitly. In particular, the first few cumulants read The second cumulant turns out to have a particularly simple form due to its close relation to the ground-state fidelity susceptibility [74,94] and reads Figure 10: Relative error for the first two cumulants of the magnetization distribution as a function of magnetic field g. Using the full characteristic function in Eq. (101) and the corresponding PPA Eq. (105), the relative error is evaluated as a function of the magnetic field for a system size L = 12 and different temperatures.
By contrast, in the infinite-temperature case, in which the PPA is exact, the distribution is symmetric, centered at m = 0 and independent of the magnetic field. The magnetization distribution describes in this case the sum of L independent discrete random variables with outcomes ±1 with equal probability 1/2. As a result κ 1 = 0, κ 2 = L/4. In the infinite temperature limit, P(m) is equal to that of a classical Ising chain and can be written explicitly: As a result, in the normal approximation P β→0 (n = 2l) is given by Gaussian distribution with zero mean and variance κ 2 = L. Fig. 11 shows this Gaussian distribution as a black envelope, accurately approximating the exact results.

Conclusion
We have provided an exact treatment of the thermal equilibrium properties for a class of integrable spin chains that admit a description in terms of free fermions. Instances of this family are the onedimensional transverse-field Ising, XY and Kitaev models, among other examples. Whenever the system Hamiltonian commutes with parity operator, the complete Hilbert spaces is the direct sum of the corresponding even and odd parity subspaces. For an exact treatment of thermal equilibrium, we have detailed an algebraic approach in the complete Hilbert spaces and provided the exact expression for the partition function. We have identified the limitations of the approximate description of thermal equilibrium in terms of the positive-parity sector, frequently adopted in the literature. This approximate approach fails in what can be considered the most interesting regime: the neighborhood of a quantum critical point at low temperatures. In particular, we have shown that the discrepancies between the exact and approximate results can lead to significant errors in this regime. Making use of the exact algebraic framework, we have computed as well the eigenvalue probability distribution of different observables. As an application, we have characterized in detail the distribution of the number of kinks as well as the transverse magnetization, covering all regimes from zero temperature (ground-state behavior) to infinite temperature. Our results are of direct relevance to the study of thermal equilibrium properties of integrable spin chains as well as the study of the nonequilibrium dynamics generated by driving a thermal state out of equilibrium. They are thus expected to find applications in the description of quantum simulation experiments, quantum annealing and quantum thermodynamics of spin systems. As a prospect, it is interesting to extend our results to the generalized Gibbs state whenever the relaxing dynamics of an initial state preserves a set of integrals of motion.

A Proof Proposition 2: Identities for Traces
First, the formulas given by Eq. (34) are true for n = 1. We assume that they are true for some n ≥ 1 and we compute tr and an inductive step is completed.