Exact Thermal Properties of Integrable Spin Chains

Michał Białończyk1, Fernando J. Gómez-Ruiz2,3, Adolfo del Campo4,2,5,6,*

1 Institute of Physics, Jagiellonian University, Łojasiewicza 11, 30-348
Kraków, Poland

2 Donostia International Physics Center, E-20018 San Sebastián, Spain

3 Departamento de Física, Universidad de Los Andes, A.A. 4976, Bogotá,
Colombia

4 Department of Physics and Materials Science, University of Luxembourg,
L-1511 Luxembourg, G. D. Luxembourg

5 IKERBASQUE, Basque Foundation for Science, E-48013 Bilbao, Spain

6 Department of Physics, University of Massachusetts, Boston, MA 02125, USA

*<EMAIL_ADDRESS>
## Abstract

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
an arbitrary observable at thermal equilibrium and characterize in detail the
thermal distribution of the kink number and transverse magnetization in the
transverse-field quantum Ising chain.

###### Contents

  1. 1 Introduction
  2. 2 Full Diagonalization of Spin-$\frac{1}{2}$ XY Model
    1. 2.1 Mathematical tools for the complete Hilbert space
  3. 3 The Canonical Partition Function
  4. 4 Full Counting Statistics in Integrable Spin Chains
    1. 4.1 Probability distribution of the number of kinks at thermal equilibrium
    2. 4.2 Probability distribution for the transverse magnetization at thermal equilibrium
  5. 5 Conclusion
  6. A Proof Proposition 2: Identities for Traces

##  1 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]. 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 [4, 5]. These family of
quasi-free fermion models include as well the Kitaev spin model in one
dimension and on a honeycomb lattice [6], among other examples [1, 2, 3].

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 [7, 8, 9]. 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 [10,
11], as reported in a variety of experiments [12, 13, 14, 15]. These results
have also been extended to nonlinear quenches [16, 17] and inhomogeneous
systems [18, 19, 20, 21, 22], while their breakdown has been characterized in
open systems [23, 24, 25]. 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 [26, 27, 28, 15,
29, 30]. 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 [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41], 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 [42, 43, 44,
45] and as a working substance in a quantum thermodynamic cycle [46].

Quasi-free fermion models provided an effective description of a variety of
condensed-matter systems, where they can be realized with high accuracy in
[47]. They are further amenable to quantum simulation with trapped ions [48,
49, 50, 51, 52], ultracold gases in optical lattices [53] and superconducting
qubits [54]. Digital quantum simulation provides yet another avenue for their
study in the laboratory [55, 56, 57, 14].

In many applications, it is generally desirable to consider a thermal state
and analyze the finite-temperature 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 positive-parity subspace, while disregarding the rest of the
spectrum [58, 1, 2, 3, 59, 60, 61, 62, 43]. 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]. However, an exact treatment requires taking into
account parity properly and at finite temperature both subspaces are
populated. Kapitonov and Il’inskii provided a derivation of the closed form
expression of the exact partition function using functional integrals over
Grassmann variables [64]. More recently, Fei and Quan [44] used group theory
methods to calculate the exact partition function and quantum work
distribution.

In this manuscript, we present an elementary derivation of the exact partition
function and generalize the results from [64], giving a precise and clear
prescription to characterize the eigenvalue distribution 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 [65, 66]. 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.

##  2 Full Diagonalization of Spin-$\frac{1}{2}$ XY Model

We consider the anisotropic one-dimensional XY Hamiltonian for spins $1/2$ in
a transverse magnetic field $g$. The Hamiltonian reads:

$\hat{\mathcal{H}}\left(g,\gamma\right)=-J\left[\sum_{n=1}^{L}\left(\frac{1+\gamma}{2}\right)\hat{X}_{n}\hat{X}_{n+1}+\left(\frac{1-\gamma}{2}\right)\hat{Y}_{n}\hat{Y}_{n+1}+g\hat{Z}_{n}\right].$
(1)

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 $\gamma>0$ and $L$ is the number of sites in the chain. For $\gamma=1$, the
Hamiltonian (1) corresponds to the Ising model in a transverse magnetic field,
which possesses a $\mathbb{Z}_{2}$ symmetry. The limit $\gamma=0$ describes
the isotropic XY model. For the anisotropic case $0<\gamma\leq 1$ the model
belongs to the Ising universality class, and its phase diagram is determined
by the ratio $\nu=g/J$. When $\nu>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\leq\nu<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 $\nu=1$. Finally, for the isotropic case $\gamma=0$, a QPT is
observed between gapless ($\nu<1$) and ferromagnetic ($\nu>1$) phases.
The operators $\hat{X}_{n}$, $\hat{Y}_{n}$, and $\hat{Z}_{n}$ are matrices of
order $2^{L}$ defined by the relations

$\begin{split}\hat{X}_{n}&=\hat{\mathbb{I}}_{1}\otimes\ldots\otimes\hat{\mathbb{I}}_{n-1}\otimes\hat{\sigma}_{n}^{x}\otimes\hat{\mathbb{I}}_{n+1}\otimes\ldots\otimes\hat{\mathbb{I}}_{L},\\\
\hat{Y}_{n}&=\hat{\mathbb{I}}_{1}\otimes\ldots\otimes\hat{\mathbb{I}}_{n-1}\otimes\hat{\sigma}_{n}^{y}\otimes\hat{\mathbb{I}}_{n+1}\otimes\ldots\otimes\hat{\mathbb{I}}_{L},\\\
\hat{Z}_{n}&=\hat{\mathbb{I}}_{1}\otimes\ldots\otimes\hat{\mathbb{I}}_{n-1}\otimes\hat{\sigma}_{n}^{z}\otimes\hat{\mathbb{I}}_{n+1}\otimes\ldots\otimes\hat{\mathbb{I}}_{L}.\end{split}$
(2)

Here, $\hat{\sigma}_{n}^{\alpha}$ denotes the Pauli operator at site $n$ along
the axis $\alpha=x,y,z$, $\hat{\mathbb{I}}_{n}$ is the identity matrix of
order $2$ at the site $n$, and periodic boundary conditions are assumed,
$\hat{\sigma}_{L+1}^{\alpha}=\hat{\sigma}_{1}^{\alpha}$. A standard way to
diagonalize the Hamiltonian in Eq. (1) relies on introducing a new set of
Fermionic operators given by

$\begin{split}\hat{\sigma}_{n}^{x}&=\left(\hat{c}_{n}^{\dagger}+\hat{c}_{n}\right)\prod_{m<n}\left(\hat{\mathbb{I}}-2\hat{c}_{m}^{\dagger}\hat{c}_{m}\right),\\\
\hat{\sigma}_{n}^{y}&=-i\left(\hat{c}_{n}^{\dagger}-\hat{c}_{n}\right)\prod_{m<n}\left(\hat{\mathbb{I}}-2\hat{c}_{m}^{\dagger}\hat{c}_{m}\right),\\\
\hat{\sigma}_{n}^{z}&=\hat{\mathbb{I}}-2\hat{c}_{n}^{\dagger}\hat{c}_{n}.\end{split}$
(3)

These expressions represent the well-known Jordan-Wigner transformation [67].
Here, $\hat{c}_{n}$ and $\hat{c}_{n}^{\dagger}$ are ladder Fermionic operators
at site $n$, which satisfy anti-commutation relations
$\left\\{\hat{c}_{i},\hat{c}_{j}^{\dagger}\right\\}=\delta_{i,j}$ and
$\left\\{\hat{c}_{i},\hat{c}_{j}\right\\}=\left\\{\hat{c}_{i}^{\dagger},\hat{c}_{j}^{\dagger}\right\\}=0$.
This is in contrast to the Pauli matrices, which satisfy commutation relations
$\left[\hat{\sigma}_{n}^{\dagger},\hat{\sigma}_{m}^{-}\right]=\delta_{n,m}\hat{\sigma}_{n}^{z}$
and $\left[\hat{\sigma}_{n}^{z},\hat{\sigma}_{m}^{\pm}\right]=\pm
2\delta_{n,m}\hat{\sigma}_{n}^{\pm}$ with
$\hat{\sigma}_{n}^{\pm}=\hat{\sigma}_{n}^{x}\pm i\hat{\sigma}_{n}^{y}$. With
periodic boundary conditions in the spin representation, the Fermionic
operators $\hat{c}_{n}$ and $\hat{c}_{n}^{\dagger}$ satisfy nontrivial
boundary conditions

$\displaystyle\hat{c}_{L+1}^{\dagger}$
$\displaystyle=\left(-1\right)^{\hat{{\rm N}}}\hat{c}_{1}^{\dagger},$
$\displaystyle\hat{c}_{L+1}$ $\displaystyle=\left(-1\right)^{\hat{{\rm
N}}}\hat{c}_{1},$ (4)

where $\hat{{\rm N}}=\sum_{n=1}^{L}\hat{c}_{n}^{\dagger}\hat{c}_{n}$ is the
Fermionic number operator. By direct substitution of Eq. (3) into Eq. (1), the
Hamiltonian can be written as a quadratic form

$\begin{split}\hat{H}\left(g,\gamma\right)&=-\sum_{n=1}^{L-1}\left[\hat{c}_{n}^{\dagger}\hat{c}_{n+1}+\hat{c}_{n+1}^{\dagger}\hat{c}_{n}+\gamma\left(\hat{c}_{n}^{\dagger}\hat{c}_{n+1}^{\dagger}+\hat{c}_{n+1}\hat{c}_{n}\right)\right]\\\
&\quad+\hat{\Pi}\left[\hat{c}_{L+1}^{\dagger}\hat{c}_{1}+\hat{c}_{1}^{\dagger}\hat{c}_{L+1}+\gamma\left(\hat{c}_{L+1}\hat{c}_{1}+\hat{c}_{1}\hat{c}_{L+1}\right)\right]-g\sum_{n=1}^{L}\left(1-2\hat{c}_{n}^{\dagger}\hat{c}_{n}\right).\end{split}$
(5)

Here, the parity operator $\hat{\Pi}$ is given by $\left(-1\right)^{\hat{{\rm
N}}}=\exp\left(i\pi\hat{{\rm N}}\right)$ and has eigenvalues $\pm 1$. The
parity operator anticommmutes with the creation $\hat{c}^{\dagger}_{n}$ and
annihilation $\hat{c}_{n}$ Fermionic operators,
$\left\\{\left(-1\right)^{\hat{{\rm
N}}},\hat{c}_{n}^{\dagger}\right\\}=\left\\{\left(-1\right)^{\hat{{\rm
N}}},\hat{c}_{n}\right\\}=0$, and therefore, it commutes with any operator
bilinear in $\hat{c}_{n}^{\dagger}$ and $\hat{c}_{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 $\mathbb{Z}_{2}$ symmetry and, thus,
the parity operator $\hat{\Pi}$ 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 $\left(+1\right)$ and negative
$\left(-1\right)$ parity. Using the projectors $\hat{\Pi}^{\pm}$,

$\hat{\Pi}^{\pm}=\frac{1}{2}\left[\hat{\mathbb{I}}\pm\left(-1\right)^{\hat{{\rm
N}}}\right],$ (6)

the Hamiltonian in Eq. (5) is represented in the form

$\hat{H}=\hat{H}^{+}\hat{\Pi}^{+}+\hat{H}^{-}\hat{\Pi}^{-},$ (7)

with the reduced Hamiltonians $\hat{H}^{\pm}$ being given by

$\hat{H}^{\pm}\left(g,\gamma\right)=-\sum_{n=1}^{L}\left[\hat{c}_{n}^{\dagger}\hat{c}_{n+1}+\hat{c}_{n+1}^{\dagger}\hat{c}_{n}+\gamma\left(\hat{c}_{n}^{\dagger}\hat{c}_{n+1}^{\dagger}+\hat{c}_{n+1}\hat{c}_{n}\right)+g\left(1-2\hat{c}_{n}^{\dagger}\hat{c}_{n}\right)\right].$
(8)

A subtle difference between $\hat{H}^{+}$ and $\hat{H}^{-}$ is found in the
boundary conditions for the Fermion operators. $\hat{H}^{+}$ obeys
antiperiodic boundary conditions ($\hat{c}_{L+1}=-\hat{c}_{1}$ and
$\hat{c}_{L+1}^{\dagger}=-\hat{c}_{1}^{\dagger}$) while $\hat{H}^{-}$
satisfies periodic boundary conditions ($\hat{c}_{L+1}=\hat{c}_{1}$ and
$\hat{c}_{L+1}^{\dagger}=\hat{c}_{1}^{\dagger}$). The Hamiltonian given by Eq.
(8) is quadratic in the Fermionic operators and is thus exactly diagonalizable
using Fourier and Bogoliubov transformations [68, 58, 69, 70]. We expand the
operator $\hat{c}_{n}$ via a Fourier transformation in momentum space,

$\displaystyle\hat{c}_{n}$
$\displaystyle=\frac{e^{-i\pi/4}}{\sqrt{L}}\sum_{k\in\mathbf{K}^{\pm}}\hat{c}_{k}\exp\left(ink\right),$
$\displaystyle\hat{c}_{n}^{\dagger}$
$\displaystyle=\frac{e^{i\pi/4}}{\sqrt{L}}\sum_{k\in\mathbf{K}^{\pm}}\hat{c}_{k}^{\dagger}\exp\left(-ink\right).$
(9)

The wavevector $k$ takes values in the positive $\left(\mathbf{K}^{+}\right)$
and negative $\left(\mathbf{K}^{-}\right)$ parity sectors

$\displaystyle\mathbf{K}^{+}$
$\displaystyle=\Bigg{\\{}k\left|\frac{\pi}{L}\left(2m-1\right)\right.,$
$\displaystyle
m=-\frac{L}{2}+1,-\frac{L}{2}+2,\ldots,\frac{L}{2}-1,\frac{L}{2}\Bigg{\\}},$
(10) $\displaystyle\mathbf{K}^{-}$
$\displaystyle=\Bigg{\\{}k\left|\frac{2\pi}{L}m\right.,$ $\displaystyle
m=-\frac{L}{2}+1,-\frac{L}{2}+2,\ldots,\frac{L}{2}-1,\frac{L}{2}\Bigg{\\}}.$
(11)

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 $\mathbf{k}=0$ and $\mathbf{k}=\pi$ are included in the
negative parity sector. For even $L$, we can rewrite conveniently the momentum
values as

$\displaystyle\mathbf{K}^{+}$
$\displaystyle=\left\\{\pm\frac{\pi}{L},\pm\frac{3\pi}{L},\pm\frac{5\pi}{L},\ldots,\pm\frac{\pi\left(L-1\right)}{L}\right\\}=\mathbf{k}^{+}\cup\\{-\mathbf{k}^{+}\\},$
$\displaystyle\mathbf{K}^{-}$
$\displaystyle=\left\\{0,\pm\frac{2\pi}{L},\pm\frac{4\pi}{L},\ldots,\pm\frac{\pi\left(L-2\right)}{L},\pi\right\\}=\mathbf{k}^{-}\cup\\{-\mathbf{k}^{-}\\}\cup\\{0,\pi\\},$

with

$\displaystyle\mathbf{k}^{+}=\left\\{\frac{\pi}{L},\frac{3\pi}{L},\ldots,\frac{\pi\left(L-1\right)}{L}\right\\},$
and
$\displaystyle\mathbf{k}^{-}=\left\\{\frac{2\pi}{L},\frac{4\pi}{L},\ldots,\frac{\pi\left(L-2\right)}{L}\right\\}.$
(12)

By direct substitution of Eq. (9) into Eq. (8), the reduced Hamiltonians
$\hat{H}^{+}$ and $\hat{H}^{-}$ are expressed in terms of $\hat{c}_{k}$ and
$\hat{c}^{\dagger}_{k}$ as

$\begin{split}\hat{H}^{+}\left(g,\gamma\right)&=\sum_{k\in\mathbf{k}^{+}}\hat{H}_{k}\left(g,\gamma\right),\\\
\hat{H}^{-}\left(g,\gamma\right)&=\sum_{k\in\mathbf{k}^{-}}\hat{H}_{k}\left(g,\gamma\right)+\hat{H}_{0}\left(g\right)+\hat{H}_{\pi}\left(g\right),\end{split}$
(13)

where

$\begin{split}\hat{H}_{k}\left(g,\gamma\right)&=2\left[\left(g-\cos\left(k\right)\right)\left(\hat{c}_{k}^{\dagger}\hat{c}_{k}-\hat{c}_{-k}\hat{c}_{-k}^{\dagger}\right)\right.+\left.\gamma\sin\left(k\right)\left(\hat{c}_{k}^{\dagger}\hat{c}_{-k}^{\dagger}-\hat{c}_{-k}\hat{c}_{k}\right)\right],\\\
\hat{H}_{0}\left(g\right)&=\left(g-1\right)\left(\hat{c}^{\dagger}_{0}\hat{c}_{0}-\hat{c}_{0}\hat{c}_{0}^{\dagger}\right),\\\
\hat{H}_{\pi}\left(g\right)&=\left(g+1\right)\left(\hat{c}^{\dagger}_{\pi}\hat{c}_{\pi}-\hat{c}_{\pi}\hat{c}_{\pi}^{\dagger}\right).\end{split}$
(14)

We next make use of a Bogoliubov transformation, and define a new set of
fermion operators $\hat{\gamma}_{k}$ and $\hat{\gamma}_{k}^{\dagger}$ given by

$\displaystyle\hat{\gamma}_{k}$
$\displaystyle=u_{k}\hat{c}_{k}-iv_{k}\hat{c}_{-k}^{\dagger},$
$\displaystyle\hat{\gamma}_{k}^{\dagger}=u_{k}\hat{c}_{k}^{\dagger}+iv_{k}\hat{c}_{-k},$
(15)

where the real numbers $u_{k}$ and $v_{k}$ satisfy $u_{k}=u_{-k}$,
$v_{k}=-v_{-k}$ and $\left|u_{k}\right|^{2}+\left|v_{k}\right|^{2}=1$. The
canonical anti-commutation relations for the operators $\hat{c}_{k}$ and
$\hat{c}_{k}^{\dagger}$ imply that the same relations are also satisfied by
$\hat{\gamma}_{k}$ and $\hat{\gamma}_{k}^{\dagger}$, that is,
$\left\\{\hat{\gamma}_{k},\hat{\gamma}_{k^{\prime}}^{\dagger}\right\\}=\delta_{k,k^{\prime}}$,
and
$\left\\{\hat{\gamma}_{k}^{\dagger},\hat{\gamma}_{k^{\prime}}^{\dagger}\right\\}=\left\\{\hat{\gamma}_{k},\hat{\gamma}_{k^{\prime}}\right\\}=0$.
By direct substitution of the Bogoliubov transformations into Eq. (13), after
a some algebra, we obtain

$\begin{split}\hat{H}_{k}\left(g,\gamma\right)&=2\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{k}\left[u_{k}^{2}\left(\cos\left(k\right)-g\right)+\gamma\sin\left(k\right)u_{k}v_{k}\right]\\\
&\qquad+2\hat{\gamma}_{k}\hat{\gamma}_{k}^{\dagger}\left[\left(\cos\left(k\right)-g\right)v_{k}^{2}-\gamma\sin\left(ak\right)u_{k}v_{k}\right]\\\
&\qquad-i\hat{\gamma}_{k}\hat{\gamma}_{-k}\left[\gamma\sin\left(k\right)\left(u_{k}^{2}-v_{k}^{2}\right)+2\left(\cos\left(k\right)-g\right)u_{k}v_{k}\right]\\\
&\qquad-i\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{-k}^{\dagger}\left[\gamma\sin\left(k\right)\left(u_{k}^{2}-v_{k}^{2}\right)+2\left(\cos\left(k\right)-g\right)u_{k}v_{k}\right]+g.\end{split}$
(16)

The terms proportional to $\gamma^{\dagger}_{k}\gamma^{\dagger}_{-k}$ and
$\gamma_{k}\gamma_{-k}$ should vanish for the Hamiltonian to acquire a
diagonal form. Writing $u_{k}=\cos\left(\vartheta_{k}/2\right)$ and
$v_{k}=\sin\left(\vartheta_{k}/2\right)$, the Bogoliubov angles satisfy

$\tan\left(\vartheta_{k}\right)=\frac{\gamma\sin\left(k\right)}{g-\cos\left(k\right)}.$
(17)

For numerical simulations, the last condition can be rewritten as
$\gamma\sin\left(k\right)\left\\{u_{k}^{2}-v_{k}^{2}\right\\}+2\left(\cos\left(k\right)-g\right)u_{k}v_{k}=0$.
Finally, the Hamiltonian (13) can be rewritten as a sum of noninteracting
terms

$\begin{split}\hat{H}^{+}\left(g,\gamma\right)&=\sum_{k\in\mathbf{k}^{+}}\epsilon_{k}\left(g,\gamma\right)\left(\hat{n}_{k}+\hat{n}_{-k}-1\right),\\\
\hat{H}^{-}\left(g,\gamma\right)&=\sum_{k\in\mathbf{k}^{-}}\epsilon_{k}\left(g,\gamma\right)\left(\hat{n}_{k}+\hat{n}_{-k}-1\right)+\left(g-1\right)\left(2\hat{n}_{0}-1\right)+\left(g+1\right)\left(2\hat{n}_{\pi}-1\right),\end{split}$
(18)

with $\hat{n}_{k}=\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{k}$ denoting the
fermion number operator and
$\epsilon_{k}\left(g,\gamma\right)=2\sqrt{\left(g-\cos
k\right)^{2}+\gamma^{2}\sin^{2}k}$ being the quasiparticle energy of mode
$\mathbf{k}\neq 0,\pi$ per particle.

###  2.1 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
$\varmathbb{H}_{k}$.

  * •

Operators are denoted by letters with a hat, such as $\hat{O}_{k}$ and
$\hat{h}_{k_{i}}$.

  * •

Operations on tensor products of Hilbert spaces are denoted with calligraphic
letters $\mathcal{P}$ and $\mathcal{N}$.

To begin with, we note that the positive-parity Hilbert subspace
$\varmathbb{H}^{+}$ can be written as the tensor product of subspaces
corresponding to each _pair of momenta_ ($k$ _and_ $-k$)

$\varmathbb{H}^{+}=\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}.$ (19)

Each subspace $\varmathbb{H}_{k}$ is the linear span of the vacuum and states
involving one and two Fermionic excitations with a given momentum

$\begin{split}\varmathbb{H}_{k}&=\mathrm{span}\\{\ket{0}_{k},\hat{c}_{k}^{\dagger}\hat{c}_{-k}^{\dagger}\ket{0}_{k},\hat{c}_{k}^{\dagger}\ket{0}_{k},\hat{c}_{-k}^{\dagger}\ket{0}_{k}\\}\\\
&=\\{\ket{00}_{k},\ket{11}_{k},\ket{10}_{k},\ket{01}_{k}\\},\quad\forall\quad
k\in\mathbf{k}^{+}.\end{split}$ (20)

Here, $\ket{0}_{k}$ is the vector annihilated by both $\hat{c}_{k}$ and
$\hat{c}_{-k}$. Each of the subspaces can be divided into the sectors with
even $\varmathbb{H}_{k}^{(p)}$ and odd $\varmathbb{H}_{k}^{(n)}$ number of
excitations

$\begin{split}\varmathbb{H}_{k}^{(p)}&=\mathrm{span}\\{\ket{0}_{k},\hat{c}_{k}^{\dagger}\hat{c}_{-k}^{\dagger}\ket{0}_{k}\\}=\\{\ket{00}_{k},\ket{11}_{k}\\},\\\
\varmathbb{H}_{k}^{(n)}&=\mathrm{span}\\{\hat{c}_{-k}^{\dagger}\ket{0}_{k},\hat{c}_{k}^{\dagger}\ket{0}_{k}\\}=\\{\ket{01}_{k},\ket{10}_{k}\\}.\end{split}$
(21)

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 $\hat{\Pi}$ has eigenvalue $+1$.
Thus, the subspace is only spanned by vectors associated with an even number
of quasiparticles. We denote this subspace by
$\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right)$

$\mathcal{P}=\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right)=\mathrm{span}\left\\{\bigotimes_{k\in\mathbf{k}^{+}}\ket{i_{k}j_{k}}\colon\;i_{k},j_{k}\in\\{0,1\\},\;\sum_{k\in\mathbf{k}^{+}}(i_{k}+j_{k})\,\,\text{is
even}\right\\}.$ (22)

Similarly, we define the subspace spanned by odd number of quasi-particle
excitations and denote it by
$\mathcal{N}=\mathcal{N}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right)$.
It is easy to see that both spaces
$\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right)$ and
$\mathcal{N}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right)$
have dimension $2^{L-1}$ and satisfy

$\varmathbb{H}^{+}=\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right)\oplus\mathcal{N}\left(\bigotimes_{k\in\mathbf{k}^{+}}\varmathbb{H}_{k}\right).$
(23)

For the positive-parity subspace only $\mathcal{P}$ is relevant; vectors in
$\mathcal{N}$ have no physical meaning for the system described by the
Hamiltonian $\hat{H}^{+}$. However, the spaces $\mathcal{P}$ and $\mathcal{N}$
(defined for proper momenta) exchange their roles for $\hat{H}^{-}$; 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 $\mathcal{P}$. This can be done for states
and observables. Before dealing with observables, we introduce an alternative
recursive definition of the spaces $\mathcal{P}$ and $\mathcal{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),

$\mathcal{P}\left(\varmathbb{H}_{k_{1}}\right)=\varmathbb{H}_{k_{1}}^{(p)},\quad\mathcal{N}\left(\varmathbb{H}_{k_{1}}\right)=\varmathbb{H}_{k_{1}}^{(n)}.$
(24)

Next, we specify how to construct spaces $\mathcal{P}$ and $\mathcal{N}$ when
a mode with momentum $k_{n+1}$ is added:

$\begin{split}\mathcal{P}\left(\bigotimes_{i=1}^{n+1}\varmathbb{H}_{k_{i}}\right)&=\mathcal{P}\left(\bigotimes_{i=1}^{n}\varmathbb{H}_{k_{i}}\right)\otimes\varmathbb{H}_{k_{n+1}}^{(p)}\oplus\mathcal{N}\left(\bigotimes_{i=1}^{n}\varmathbb{H}_{k_{i}}\right)\otimes\varmathbb{H}_{k_{n+1}}^{(n)},\,\,\,n\geq
1,\\\
\mathcal{N}\left(\bigotimes_{i=1}^{n+1}\varmathbb{H}_{k_{i}}\right)&=\mathcal{N}\left(\bigotimes_{i=1}^{n}\varmathbb{H}_{k_{i}}\right)\otimes\varmathbb{H}_{k_{n+1}}^{(p)}\oplus\mathcal{P}\left(\bigotimes_{i=1}^{n}\varmathbb{H}_{k_{i}}\right)\otimes\varmathbb{H}_{k_{n+1}}^{(n)},\,\,\,n\geq
1.\end{split}$ (25)

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 $\hat{O}_{k}$ act independently on each subspace
$\varmathbb{H}_{k}$ and each $\hat{O}_{k}$ can be written as a sum of an even
part $\hat{O}_{k}^{(p)}$ and an odd part $\hat{O}_{k}^{(n)}$ as

$\hat{O}_{k}=\hat{O}_{k}^{(p)}+\hat{O}_{k}^{(n)},\quad\hat{O}_{k}^{(p)}\bigg{|}_{\varmathbb{H}_{k}^{(n)}}=0,\quad\hat{O}_{k}^{(n)}\bigg{|}_{\varmathbb{H}_{k}^{(p)}}=0.$
(26)

The operators $\hat{O}_{k}^{(p)}$ and $\hat{O}_{k}^{(n)}$ act on the total
space $\varmathbb{H}_{k}$, but have a $2\times 2$ zero block $0_{2}$ in the
respective subspace. The proper restrictions of the tensor product of
operators $\hat{O}_{k}$ can be defined in a similar way as in Eqs. (24) and
(25) for $\mathcal{P}\left(\hat{O}_{k_{1}}\right)=\hat{O}_{k_{1}}^{(p)}$ and
$\mathcal{N}\left(\hat{O}_{k_{1}}\right)=\hat{O}_{k_{1}}^{(n)}$, and are given
by

$\begin{split}\mathcal{P}\left(\bigotimes_{i=1}^{n+1}\hat{O}_{k_{i}}\right)&=\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\otimes\hat{O}_{k_{n+1}}^{(p)}+\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\otimes\hat{O}_{k_{n+1}}^{(n)},\,\,\,n\geq
1,\\\
\mathcal{N}\left(\bigotimes_{i=1}^{n+1}\hat{O}_{k_{i}}\right)&=\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\otimes\hat{O}_{k_{n+1}}^{(p)}+\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\otimes\hat{O}_{k_{n+1}}^{(n)},\,\,\,n\geq
1.\end{split}$ (27)

Example 2.1:  Even and odd parity parts of the Hamiltonian For $\hat{H}_{k}$
given by Eq. (14), note that for a each mode $k_{n}$ the Hamiltonian can be
rewritten as
$\hat{H}_{k}=\mathcal{P}\left(\hat{\mathbb{I}}_{k_{1}}\otimes\hat{\mathbb{I}}_{k_{2}}\otimes\ldots\otimes\hat{h}_{k_{n}}\otimes\ldots\otimes\hat{\mathbb{I}}_{k_{L/2}}\right),$
where, in the basis
$\\{\ket{00}_{k},\ket{11}_{k},\ket{01}_{k},\ket{10}_{k}\\}$,
$\hat{h}_{k_{n}}=2\left(\begin{array}[]{cccc}\cos\left(k_{n}\right)-g&\gamma\sin\left(k_{n}\right)&0&0\\\
\gamma\sin\left(k_{n}\right)&g-\cos\left(k_{n}\right)&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ \end{array}\right).$ Here, $\hat{h}_{k_{n}}^{(n)}$ is $4\times 4$
zero matrix (with no odd part), and $\hat{h}_{k_{n}}^{(p)}=\hat{h}_{k_{n}}$.

As the odd part of Hamiltonian is zero, the description using ordinary tensor
products instead of over $\mathcal{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$:
$\hat{\rho}_{k}=\exp\left(-\beta\hat{h}_{k}\right).$ (28) Using the expression
for $\hat{h}_{k}$ in the the basis
$\\{\ket{00}_{k},\ket{11}_{k},\ket{01}_{k},\ket{10}_{k}\\}$,
$\hat{\rho}_{k}=\exp\left[-2\beta\left(\begin{array}[]{cc}\cos(k)-g&\gamma\sin(k)\\\
\gamma\sin(k)&g-\cos(k)\\\ \end{array}\right)\right]\oplus\mathbb{I}_{2}.$
(29) Therefore, the even and odd parts read:
$\hat{\rho}_{k}^{(p)}=\exp\left[-2\beta\left(\begin{array}[]{cc}\cos(k)-g&\gamma\sin(k)\\\
\gamma\sin(k)&g-\cos(k)\\\ \end{array}\right)\right]\oplus
0_{2},\quad\hat{\rho}_{k}^{(n)}=0_{2}\oplus\mathbb{I}_{2}.$ (30) Using the
fact that $\hat{h}_{k}$ has eigenvalues $\pm\epsilon_{k}$, we have:
$\mathrm{Tr}\left(\hat{\rho}_{k}^{(p)}\right)=2\cosh\left(\beta\epsilon_{k}\left(g,\gamma\right)\right),\quad\mathrm{Tr}\left(\rho_{k}^{(n)}\right)=2.$
(31)

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
$\hat{O}_{k}$ and $\hat{R}_{k}$ acting independently on each subspace
$\varmathbb{H}_{k}$. Then, the following identities are true for operator
multiplication
$\begin{split}\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\,\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{R}_{k_{i}}\right)&=\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\,\hat{R}_{k_{i}}\right),\\\
\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\,\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{R}_{k_{i}}\right)&=\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\,\hat{R}_{k_{i}}\right).\end{split}$
(32)

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 $\varmathbb{H}_{k}$, the following
identities for exponents of operators hold:
$\begin{split}\exp\left[\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\right]&=\mathcal{P}\left(\bigotimes_{i=1}^{n}\exp\left(\hat{O}_{k_{i}}\right)\right),\\\
\exp\left[\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\right]&=\mathcal{N}\left(\bigotimes_{i=1}^{n}\exp\left(\hat{O}_{k_{i}}\right)\right).\end{split}$
(33)

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 $\hat{O}_{k}$ that act
independently on each subspace $\varmathbb{H}_{k}$. Then, the traces of the
restricted tensor products can be expressed as follows,
$\begin{split}\mathrm{tr}\left[\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\right]&=\frac{1}{2}\left(\prod_{i=1}^{n}\mathrm{tr}\left(\hat{O}_{k_{i}}\right)+\prod_{i=1}^{n}\left(\mathrm{tr}\left(\hat{O}_{k_{i}}^{(p)}\right)-\mathrm{tr}\left(\hat{O}_{k_{i}}^{(n)}\right)\right)\right),\\\
\mathrm{tr}\left[\mathcal{N}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\right]&=\frac{1}{2}\left(\prod_{i=1}^{n}\mathrm{tr}\left(\hat{O}_{k_{i}}\right)-\prod_{i=1}^{n}\left(\mathrm{tr}\left(\hat{O}_{k_{i}}^{(p)}\right)-\mathrm{tr}\left(\hat{O}_{k_{i}}^{(n)}\right)\right)\right).\end{split}$
(34) We present a proof ot Eq. (34) in the Appendix A.

_Negative-parity subspace._ In the negative-parity subspace, all formulas
derived for the positive-parity subspace remain valid. In particular, for all
momenta $k\neq 0,\pi$ expressions from examples 2.1, 2.1 apply. The only
difference is that one has to treat carefully the parts of the Hilbert space
associated with momenta $0$ and $\pi$. They are spanned by the following
bases:

$\begin{split}\varmathbb{H}_{0}&=\mathrm{span}\\{\ket{0}_{0}\,,\,\hat{c}_{0}^{\dagger}\ket{0}_{0}\\},\\\
\varmathbb{H}_{\pi}&=\mathrm{span}\\{\ket{0}_{\pi}\,,\,\hat{c}_{\pi}^{\dagger}\ket{0}_{\pi}\\}.\end{split}$
(35)

As a result, matrices describing the Hamiltonian and Gibbs state are $2\times
2$ instead of $4\times 4$. In the following example we give formulas for the
even- and odd-parity parts of the Gibbs state in modes $k=0,\pi$:

Example 2.6:  Even- and odd-parity parts of the exact Gibbs state for $0,\pi$
momenta Using equation (16), the explicit form of the Gibbs state of the modes
with momenta $0,\pi$, in the bases
$\\{\ket{0}_{0}\,,\,\hat{c}_{0}^{\dagger}\ket{0}_{0}\\}$,
$\\{\ket{0}_{\pi}\,,\,\hat{c}_{\pi}^{\dagger}\ket{0}_{\pi}\\}$, are
respectively given by
$\hat{\rho}_{0}=\left(\begin{array}[]{cc}e^{-\beta(g-1)}&0\\\
0&e^{\beta(g-1)}\\\
\end{array}\right),\quad\hat{\rho}_{\pi}=\left(\begin{array}[]{cc}e^{-\beta(g+1)}&0\\\
0&e^{\beta(g+1)}\\\ \end{array}\right).$ (36) Thus, the corresponding even-
and odd-parity parts read $\displaystyle\hat{\rho}_{0}^{(p)}$
$\displaystyle=\left(\begin{array}[]{cc}e^{-\beta\left(g-1\right)}&0\\\ 0&0\\\
\end{array}\right),$
$\displaystyle\hat{\rho}_{\pi}^{(p)}=\left(\begin{array}[]{cc}e^{-\beta\left(g+1\right)}&0\\\
0&0\\\ \end{array}\right),$ (37e) $\displaystyle\hat{\rho}_{0}^{(n)}$
$\displaystyle=\left(\begin{array}[]{cc}0&0\\\ 0&e^{\beta\left(g-1\right)}\\\
\end{array}\right),$
$\displaystyle\hat{\rho}_{\pi}^{(n)}=\left(\begin{array}[]{cc}0&0\\\
0&e^{\beta\left(g+1\right)}\\\ \end{array}\right).$ (37j)

In closing this section, we point out that when $L$ is odd, the momenta $0$
and $\pi$ appear in the positive-parity subspace; the general formulas (24)
and (26) are always valid.

##  3 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 [71].

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 $\beta$ and characterized by the initial density operator

$\hat{\rho}_{{\rm
Gibbs}}\left(\beta,g,\gamma\right)=\frac{\exp\left(-\beta\hat{H}\left(g,\gamma\right)\right)}{{Z\left(\beta,g,\gamma\right)}},$
(38)

where $Z\left(\beta,g,\gamma\right)$ is the canonical partition function given
by

$Z\left(\beta,g,\gamma\right)={\rm
tr}\left[\exp\left(-\beta\hat{H}\left(g,\gamma\right)\right)\right].$ (39)

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 $\hat{\rho}=\exp\left(-\beta\hat{H}\right)$, where
$\hat{H}$ 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
$\hat{\rho}$ is given by

$\hat{\rho}=\exp\left[-\beta\left(\hat{H}^{+}\hat{\Pi}^{+}+\hat{H}^{-}\hat{\Pi}^{-}\right)\right]=\hat{\rho}^{+}\oplus\hat{\rho}^{-},$
(40)

where

$\displaystyle\hat{\rho}^{+}=\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\hat{\rho}_{k}\right),\qquad\hat{\rho}^{-}=\mathcal{N}\left(\bigotimes_{k\in\mathbf{k}^{-}}\hat{\rho}_{k}\otimes\hat{\rho}_{0}\otimes\hat{\rho}_{\pi}\right),$
(41)

and $\hat{\rho}_{k}$ are defined in Examples 2.1, with the sets
$\mathbf{k}^{+}$ and $\mathbf{k}^{-}$ given in Eq. (12). For these operators
the corresponding reduced partition functions are

$\displaystyle
Z^{+}(\beta,g,\gamma)=\mathrm{tr}\left[\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\hat{\rho}_{k}\right)\right],$
and $\displaystyle
Z^{-}(\beta,g,\gamma)=\mathrm{tr}\left[\mathcal{N}\left(\bigotimes_{k\in\mathbf{k}^{-}}\hat{\rho}_{k}\otimes\hat{\rho}_{0}\otimes\hat{\rho}_{\pi}\right)\right].$
(42)

For simplicity, we calculate $Z^{+}$ and $Z^{-}$ separately, and focus on
$Z^{+}$ first. Using the formulas from Example 2.1, one finds

$\begin{split}\mathrm{tr}\left(\hat{\rho}_{k}\right)&=2\cosh\left(\beta\epsilon_{k}\right)+2=4\cosh^{2}\left(\frac{\beta\epsilon_{k}}{2}\right),\\\
\mathrm{tr}\left(\hat{\rho}_{k}^{(p)}\right)-\mathrm{tr}\left(\hat{\rho}_{k}^{(n)}\right)&=2\cosh\left(\beta\epsilon_{k}\right)-2=4\sinh^{2}\left(\frac{\beta\epsilon_{k}}{2}\right).\end{split}$
(43)

Making use of the first identity in (34), we obtain an expression for
canonical partition function in the positive-parity sector

$Z^{+}\left(\beta,g,\gamma\right)=\frac{1}{2}\left(\prod_{k\in\mathbf{k}^{+}}2^{2}\cosh^{2}\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)+\prod_{k\in\mathbf{k}^{+}}2^{2}\sinh^{2}\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)\right).$
(44)

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

$\begin{split}Z^{-}\left(\beta,g,\gamma\right)&=\frac{1}{2}\Bigg{(}2^{2}\cosh\left(\beta\left(g+1\right)\right)\cosh\left(\beta\left(g-1\right)\right)\prod_{k\in\mathbf{k}^{-}}2^{2}\cosh^{2}\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)\\\
&\;\;-2^{2}\sinh\left(\beta\left(g+1\right)\right)\sinh\left(\beta\left(g-1\right)\right)\prod_{k\in\mathbf{k}^{-}}2^{2}\sinh^{2}\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)\Bigg{)}.\end{split}$
(45)

Using (40), the exact partition is the sum of contributions of positive and
negative parity:
$Z(\beta,g,\lambda)=Z^{+}(\beta,g,\gamma)+Z^{-}(\beta,g,\gamma)$. To sum up,
one can rewrite exact partition function in closed-form.

Summary 3.1:  Exact partition function for spin-$\frac{1}{2}$ XY model
$\begin{split}Z\left(\beta,g,\gamma\right)&=\frac{1}{2}\Bigg{(}\prod_{k\in\mathbf{K}^{+}}2\cosh\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)+\prod_{k\in\mathbf{K}^{+}}2\sinh\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)\\\
&\quad+\prod_{k\in\mathbf{K}^{-}}2\cosh\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)-\prod_{k\in\mathbf{K}^{-}}2\sinh\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right)\Bigg{)},\end{split}$
(46) where
$\displaystyle\epsilon_{k}\left(g,\gamma\right)=2\sqrt{\left(g-\cos\left(k\right)\right)^{2}+\left(\gamma\sin\left(k\right)\right)^{2}},\quad\epsilon_{k=0}=2(g-1),\quad\epsilon_{k=\pi}=2(g+1).$
(47)

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,

$\begin{split}Z\left(\beta,g,\gamma\right)&=\frac{1}{2}\left[Z_{F}^{+}\left(\beta,g,\gamma\right)+Z_{F}^{-}\left(\beta,g,\gamma\right)+Z_{B}^{+}\left(\beta,g,\gamma\right)-Z_{B}^{-}\left(\beta,g,\gamma\right)\right]\end{split}$
(48)

where
$Z_{F}^{\pm}\left(\beta,g,\gamma\right)=\prod_{k\in\mathbf{K}^{\pm}}2\cosh\left(\beta\epsilon_{k}\left(g,\gamma\right)/2\right)$
and
$Z_{B}^{\pm}\left(\beta,g,\gamma\right)=\prod_{k\in\mathbf{K}^{\pm}}2\sinh\left(\beta\epsilon_{k}\left(g,\gamma\right)/2\right)$
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\rightarrow\infty$ [1, 59, 60, 43, 58, 62]

Summary 3.2:  PPA partition function $Z_{\rm
PPA}(\beta,g,\gamma)=Z_{F}^{+}(\beta,g,\gamma)=\prod_{k\in\mathbf{K}^{+}}2\cosh\left(\frac{\beta}{2}\epsilon_{k}\left(g,\gamma\right)\right).$
(49)

Figure 1:  Comparison of the exact and PPA canonical partition functions. The
ratio between the total partition function 3 and the PPA Eq. (49) is shown in
the $\beta$-$g$ plane for finite system size $L=50,\,100,\,5000,\,10000$,
increasing from left to right (anisotropic parameter $\gamma=1$). Significant
differences appear close to the critical point $g=g_{c}=1$, with the magnitude
of $Z_{F}^{+}(\beta,g,\gamma)$ 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.

The complete expression for the partition function 3 was already derived using
an alternative method based on Grassmann variables, although without a
numerical characterization [64]; see as well [44]. It is thus natural to
analyze the extent to which the PPA $Z_{F}^{+}(\beta,g,\gamma)$ provides a
valid approximation to the exact partition function.
Fig. 1 shows the difference between the ratio
$Z_{F}^{+}\left(\beta\right)/Z\left(\beta\right)$ 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}^{+}(\beta,g,\gamma)\approx 0.5\,Z\left(\beta,g,\gamma\right)$.

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}^{\alpha}$ and $\ket{g^{\alpha}}$ the lowest energy level and the
corresponding eigenstate in the subspace of parity $\alpha=\pm$. The
diagonalization procedure of the Ising model yields explicit formulas for
these eigenvalues. For even number of spins [72]

$\begin{split}E_{g}^{+}&=-\sum_{k\in\mathbf{k}^{+}}\epsilon_{k},\\\
E_{g}^{-}&=-\sum_{k\in\mathbf{k}^{-}}\epsilon_{k}-2.\end{split}$ (50)

The corresponding eigenstates read

$\begin{split}\ket{g^{+}}=\prod_{k\in\mathbf{k}^{+}}(\cos(\vartheta_{k}/2)-\sin(\vartheta_{k}/2)\hat{c}_{k}^{\dagger}\hat{c}_{-k}^{\dagger})\ket{\mathrm{vac}},\\\
\ket{g^{-}}=c_{0}^{\dagger}\prod_{k\in\mathbf{k}^{-}}(\cos(\vartheta_{k}/2)-\sin(\vartheta_{k}/2)\hat{c}_{k}^{\dagger}\hat{c}_{-k}^{\dagger})\ket{\mathrm{vac}},\end{split}$
(51)

where $\ket{\mathrm{vac}}$ is annihilated by all $\hat{c}_{k}$ for
$k\in\mathbf{K}^{+}\cup\mathbf{K}^{-}$ (including $0$ and $\pi$ modes). In
what follows, we restrict ourselves to the TFQIM ($\gamma=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
[73]). The energy gap $\delta(g)$ between these two lowest energy states plays
a crucial role. We recall its asymptotic behavior [72]

$\begin{split}\delta(0<g<1)&=\mathcal{O}\left[\sim\mathrm{exp}\left(-L/\xi(g)\right)\right],\\\
\delta(g=1)&=2\mathrm{tan}\left[\frac{\pi}{4L}\right]\approx\frac{\pi}{2L},\\\
\delta(g>1)&=2g-2+\mathcal{O}\left(g^{-L}\right),\end{split}$ (52)

where $\xi(g)$ denotes the correlation length. In the low temperature regime,
the Gibbs state is effectively spanned by the two lowest energy states,
$\ket{g}^{+}$ and $\ket{g}^{-}$. In this truncation, the partition function
and Gibbs state read

$Z_{\rm approx}(\beta,g)=e^{-\beta E_{g}^{+}}+e^{-\beta E_{g}^{-}},$ (53)
$\rho_{\rm Gibbs}(\beta,g)\approx\frac{1}{Z_{\rm
approx}(\beta,g)}\left(e^{-\beta E_{g}^{+}}\ket{g^{+}}\bra{g^{+}}+e^{-\beta
E_{g}^{-}}\ket{g^{-}}\bra{g^{-}}\right).$ (54)

This low-temperature two-level approximation relies on (51) and disregards the
contribution from higher excited states, that are energetically separated from
$\ket{g^{+}}$ and $\ket{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

$\Delta(g)=4\sqrt{g^{2}-2g\cos\left(\frac{\pi}{L}\right)+1}=4|g-1|+\mathcal{O}\left(\frac{1}{L^{2}}\right),\quad
g>0,\,g\neq 1,$ (55)

while at the critical point, this gap behaves as

$\Delta(g=1)\approx\frac{4\pi}{L}.$ (56)

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 $\beta\Delta(g)\gg 1$.
Similarly, in 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,

$\frac{\Delta(g=1)}{\delta(g=1)}\xrightarrow[L\to\infty]{}\frac{1}{8}.$ (57)

However, for large $\beta$ the error is very small. The accuracy of the the
two-level approximation for different phases is shown in Fig. 2.

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
$\delta(g)$ and $\Delta(g)$ in the neighbourhood of $g_{c}=1$ are comparable,
a lower temperature is required to obtain a desired level of accuracy. For
given $\beta$, the accuracy decreases with increasing system size.

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
$\ket{g^{+}}$ and can be approximated by $e^{-\beta E_{g}^{+}}$, while the
complete partition function is approximately

$Z_{\rm approx}(\beta,g)\approx e^{-\beta E_{g}^{+}}+e^{-\beta
E_{g}^{-}}\approx 2e^{-\beta E_{g}^{+}}.$ (58)

This explains the observed error of about $50\%$ between the exact and PPA
partition functions.

##  4 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 [74]. 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 [75, 76,
77]. Distributions of other observables such as the total energy play a key
role in quantum chaos [78] and the statistics of related positive-operator
valued measures (POVMs, such as work) are at the core of fluctuation theorems
in quantum thermodynamics [79]. In the context of spin chains, the
distribution of the order parameter has long been recognized as a probe for
criticality and turbulence [80, 81, 82, 83, 84, 85, 86, 87]. 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 [26, 27, 28, 15, 29].

The full counting statistics is characterized by the probability
$P\left(\omega\right)$ to obtain the eigenvalue $\omega$ of a general operator
$\hat{W}$. It is defined as the expectation value

$P\left(\omega\right)=\bigg{\langle}\delta\left(\hat{W}-\omega\right)\bigg{\rangle},$
(59)

where the $\delta$ function is to be interpreted as a Kronecker or Dirac delta
function, depending on whether the spectrum of $\hat{W}$ is point-wise or
continuous. The angular bracket denotes the quantum expectation value with
respect to a general state characterized by a density matrix $\hat{\rho}$. We
introduce the Fourier transform representation

$P\left(\omega\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\theta\tilde{P}\left(\theta\right)\exp\left(-i\theta\omega\right),$
(60)

where $\tilde{P}\left(\theta\right)$ is the characteristic function given by

$\tilde{P}\left(\theta\right)=\mathrm{tr}\left[\hat{\rho}\exp\left(i\theta\hat{W}\right)\right].$
(61)

In cases such as the kink number and the transverse magnetization, the
eigenvalues are integers $\omega\in\mathbb{Z}$ and the range of the integral
can be restricted from $-\pi$ to $\pi$. 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 $\hat{W}$ according
to

$\langle\hat{W}^{m}\rangle=\frac{1}{i^{m}}{\left.\kern-1.2pt\frac{d^{m}}{d\theta^{m}}\tilde{P}(\theta)\vphantom{\big{|}}\right|_{\theta=0}}.\\\
$ (62)

Further, its logarithm is the cumulant generating function used to derive the
cumulants of the distribution through the identity

$\kappa_{m}=(-i)^{m}{\left.\kern-1.2pt\frac{d^{m}}{d\theta^{m}}\ln\tilde{P}(\theta)\vphantom{\big{|}}\right|_{\theta=0}}.$
(63)

The first cumulant $\kappa_{1}$ is just the mean value, $\kappa_{2}$ is the
variance, and $\kappa_{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 $\kappa_{1},\kappa_{2}\neq 0$ and
vanishing $\kappa_{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
$\mathcal{W}$ of observables. This class is defined by the property that any
operator $\hat{W}\in\mathcal{W}$, in each parity subspace, can be written in
the form

$\hat{W}=\sum_{k}\hat{W}_{k},$ (64)

where

$\hat{W_{k}}=\hat{\Psi}_{k}^{\dagger}\hat{w}_{k}\hat{\Psi}_{k},\quad\hat{\Psi}^{\dagger}=\left(\hat{c}_{-k},\>\hat{c}_{k}^{\dagger},\>\hat{c}_{k},\>\hat{c}_{-k}^{\dagger}\right)$
(65)

and the matrix $\hat{w}_{k}$ has the block-diagonal form

$\hat{w}_{k}=\left(\begin{array}[]{cc}\hat{w}_{k}^{(1)}&0\\\
0&\hat{w}_{k}^{(2)}\\\ \end{array}\right).$ (66)

Here, $\hat{w}_{k}^{(1)}$ and $\hat{w}_{k}^{(2)}$ are $2\times 2$ are matrices
for momenta different from $0,\pi$ and $1\times 1$ matrices for $0,\pi$
momenta. We point out that the notation in equations (65, 66) is compatible
with matrix expressions from Examples 2.1, 2.1, written in the basis
$\\{\ket{00}_{k},\ket{11}_{k},\ket{01}_{k},\ket{10}_{k}\\}$. When off-diagonal
blocks vanish, the operator $\hat{W}_{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 $\mathcal{W}$ as
these observables mix the subspaces with different parities. The treatment of
such operators is beyond the scope of this paper. Some examples of the
intricacies involved in characterizing the longitudinal magnetization in spin
systems can be found in [61, 88, 89, 87].

In the following we present the detailed procedure for computing
characteristic function $\tilde{P}\left(\theta\right)$ of a given observable
$\hat{W}$ in the class $\mathcal{W}$.

  1. 1.

First, we fix the state $\hat{\rho}$ to be the thermal-equilibrium Gibbs
state, $\hat{\rho}=\hat{\rho}_{{\rm Gibbs}}$ given by equation (38). Then,
using formulas from Section 2 we can diagonalize the even-parity part of
$\hat{\rho}_{k}$

$\exp\left[-2\beta\left(\begin{array}[]{cc}\cos(k)-g&\gamma\sin(k)\\\
\gamma\sin(k)&g-\cos(k)\\\
\end{array}\right)\right]=\hat{S}_{k}^{\dagger}\,\mathrm{diag}\left(e^{-\beta\epsilon_{k}(g,\gamma)},\,e^{\beta\epsilon_{k}(g,\gamma)}\right)\hat{S}_{k},$
(67)

where

$\hat{S}_{k}=\left(\begin{array}[]{cc}\cos\left(\frac{\vartheta_{k}}{2}\right)&\sin\left(\frac{\vartheta_{k}}{2}\right)\\\
\sin\left(\frac{\vartheta_{k}}{2}\right)&-\cos\left(\frac{\vartheta_{k}}{2}\right)\\\
\end{array}\right)$ (68)

and the angle $\vartheta_{k}$ satisfies

$\cos\left(\vartheta_{k}\right)=\frac{2(\cos(k)-g)}{\epsilon_{k}(g,\gamma)},\quad\sin\left(\vartheta_{k}\right)=\frac{2\gamma\sin(k)}{\epsilon_{k}(g,\gamma)}.$
(69)

  2. 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:

$\tilde{P}\left(\theta\right)=\frac{1}{Z(\beta,g,\gamma)}\left(\tilde{P}^{+}\left(\theta\right)+\tilde{P}^{-}\left(\theta\right)\right).$
(70)

Using Propositions 2.1 and 2.1, we aim at calculating

$\displaystyle\tilde{P}^{+}\left(\theta\right)=\mathrm{tr}\left[\mathcal{P}\left(\bigotimes_{k\in\mathbf{k}^{+}}\hat{\rho}_{k}\exp\left(i\theta\hat{w}_{k}\right)\right)\right],$
$\displaystyle\tilde{P}^{-}\left(\theta\right)=\mathrm{tr}\left[\mathcal{N}\left(\bigotimes_{k\in\mathbf{k}^{-}}\hat{\rho}_{k}\exp\left(i\theta\hat{w}_{k}\right)\right)\right].$
(71)

Next, we define the matrix

$\hat{\sigma}_{k}=\hat{S}_{k}\exp\left(i\theta\hat{w}_{k}^{(1)}\right)\hat{S}_{k}^{\dagger}.$
(72)

Denoting the eigenvalues of $\hat{w}_{k}^{(2)}$ by $\mu_{k}$ and $\lambda_{k}$
we find

$\begin{split}\mathrm{tr}\left[\hat{\rho}_{k}\exp\left(i\theta\hat{w}_{k}\right)\right]=\hat{\sigma}_{k}^{11}\,e^{-\beta\epsilon_{k}(g,\gamma)}+\hat{\sigma}_{k}^{22}\,e^{\beta\epsilon_{k}(g,\gamma)}+e^{i\theta\mu_{k}}+e^{i\theta\lambda_{k}},\\\
\mathrm{tr}\left[\hat{\rho}_{k}^{(p)}\exp\left(i\theta\hat{w}_{k}^{(p)}\right)\right]-\mathrm{tr}\left[\hat{\rho}_{k}^{(n)}\exp\left(i\theta\hat{w}_{k}^{(n)}\right)\right]=\hat{\sigma}_{k}^{11}\,e^{-\beta\epsilon_{k}(g,\gamma)}+\hat{\sigma}_{k}^{22}\,e^{\beta\epsilon_{k}(g,\gamma)}-e^{i\theta\mu_{k}}-e^{i\theta\lambda_{k}}.\end{split}$
(73)

Using Proposition 2.1 we obtain

$\begin{split}2\tilde{P}^{+}\left(\theta\right)=&\prod_{k\in\mathbf{k}^{+}}\left(\hat{\sigma}_{k}^{11}\,e^{-\beta\epsilon_{k}(g,\gamma)}+\hat{\sigma}_{k}^{22}\,e^{\beta\epsilon_{k}(g,\gamma)}+e^{i\theta\mu}+e^{i\theta\lambda}\right)\\\
&+\prod_{k\in\mathbf{k}^{+}}\left(\hat{\sigma}_{k}^{11}\,e^{-\beta\epsilon_{k}(g,\gamma)}+\hat{\sigma}_{k}^{22}\,e^{\beta\epsilon_{k}(g,\gamma)}-e^{i\theta\mu}-e^{i\theta\lambda}\right).\end{split}$
(74)

  3. 3.

To determine $\tilde{P}^{-}(\theta)$ it remains to compute the contributions
corresponding to $0,\pi$ momenta. Denoting

$\hat{w}_{0}=\mathrm{diag}\left(w_{0}^{1},\,w_{0}^{2}\right),\quad\hat{w}_{\pi}=\mathrm{diag}\left(w_{\pi}^{1},\,w_{\pi}^{2}\right),$
(75)

one finds

$\begin{split}\hat{\rho}_{0}\exp\left(i\theta\hat{w}_{0}\right)=\mathrm{diag}\left(e^{\beta(g-1)+i\theta
w_{0}^{1}},\,e^{-\beta(g-1)+i\theta w_{0}^{2}}\right),\\\
\hat{\rho}_{\pi}\exp\left(i\theta\hat{w}_{\pi}\right)=\mathrm{diag}\left(e^{\beta(g+1)+i\theta
w_{\pi}^{1}},\,e^{-\beta(g+1)+i\theta w_{\pi}^{2}}\right).\end{split}$ (76)

Therefore, the negative-parity part of the characteristic function is

$\begin{split}2\tilde{P}^{-}\left(\theta\right)=&\tilde{P}^{F}(\theta)\prod_{k\in\mathbf{k}^{-}}\left(\hat{\sigma}_{k}^{11}\,e^{-\beta\epsilon_{k}(g,\gamma)}+\hat{\sigma}_{k}^{22}\,e^{\beta\epsilon_{k}(g,\gamma)}+e^{i\theta\mu}+e^{i\theta\lambda}\right)\\\
&-\tilde{P}^{B}(\theta)\prod_{k\in\mathbf{k}^{-}}\left(\hat{\sigma}_{k}^{11}\,e^{-\beta\epsilon_{k}(g,\gamma)}+\hat{\sigma}_{k}^{22}\,e^{\beta\epsilon_{k}(g,\gamma)}-e^{i\theta\mu}-e^{i\theta\lambda}\right),\end{split}$
(77)

where

$\begin{split}\tilde{P}^{F}(\theta)=\left(e^{\beta(g-1)+i\theta
w_{0}^{1}}+e^{-\beta(g-1)+i\theta w_{0}^{2}}\right)\left(e^{\beta(g+1)+i\theta
w_{\pi}^{1}}+e^{-\beta(g+1)+i\theta w_{\pi}^{2}}\right),\\\
\tilde{P}^{B}(\theta)=\left(e^{\beta(g-1)+i\theta
w_{0}^{1}}-e^{-\beta(g-1)+i\theta w_{0}^{2}}\right)\left(e^{\beta(g+1)+i\theta
w_{\pi}^{1}}-e^{-\beta(g+1)+i\theta w_{\pi}^{2}}\right).\end{split}$ (78)

Note that this is not the only way to calculate the characteristic function:
instead of diagonalizing $\hat{\rho}_{k}$, one could diagonalize an observable
$\hat{w}_{k}$. However, in our approach the role of the Boltzmann factor set
by $\beta\epsilon_{k}(g,\gamma)$, 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.

###  4.1 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
[11, 26, 27, 28, 15, 29]. 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 $\gamma=1$ for simplicity. The
explicit form of kink-number operator reads

$\hat{N}=\frac{1}{2}\sum_{n=1}^{L}\left(1-\hat{X}_{n}\,\hat{X}_{n+1}\right),$
(79)

with eigenvalues $n=0,1,\dots,L$ under periodic boundary conditions.

Comparing the Ising Hamiltonian Eq. (1), with $\gamma=1$ and $g=0$, with the
Bogoliubov Hamiltonian (18) at $\gamma=1$ and $g=0$, the kink operator takes a
simple form as the sum of the number operators of quasiparticles in each
momentum [11]. 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:

$\hat{N}=\frac{L}{2}+\sum_{k}\hat{N}_{k}.$ (80)

By analogy with Eq. (65) and Eq. (66), we define a new set of operators
$\hat{n}_{k}$, $\hat{n}_{0}$, and $\hat{n}_{\pi}$; taking for any mode $k\neq
0,\pi$ the basis given by
$\\{\ket{00}_{k},\ket{11}_{k},\ket{01}_{k},\ket{10}_{k}\\}$, while selecting
for $0,\pi$ momenta the basis $\\{\ket{0}_{0},c_{0}^{\dagger}\ket{0}_{0}\\}$,
$\\{\ket{0}_{\pi},c_{\pi}^{\dagger}\ket{0}_{\pi}\\}$. Therefore, we define the
operators

$\displaystyle\hat{n}_{k}$
$\displaystyle=\left(\begin{array}[]{cccc}\cos\left(k\right)&\sin\left(k\right)&0&0\\\
\sin\left(k\right)&-\cos\left(k\right)&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{array}\right),$ $\displaystyle\hat{n}_{0}$
$\displaystyle=\left(\begin{array}[]{cc}\frac{1}{2}&0\\\ 0&-\frac{1}{2}\\\
\end{array}\right)$
$\displaystyle\hat{n}_{\pi}=\left(\begin{array}[]{cc}-\frac{1}{2}&0\\\
0&\frac{1}{2}\\\ \end{array}\right),$ (89)

and thus

$\hat{n}_{k}^{(1)}=\left(\begin{array}[]{cc}\cos\left(k\right)&\sin\left(k\right)\\\
\sin\left(k\right)&-\cos\left(k\right)\\\
\end{array}\right),\quad\hat{n}_{k}^{(2)}=0_{2}.$ (90)

Note that $\exp\left(i\theta\hat{n}_{k}^{(1)}\right)$ has the simple form

$\exp\left(i\theta\hat{n}_{k}^{(1)}\right)=\left(\begin{array}[]{cc}\cos(\theta)+i\sin(\theta)\cos(k)&i\sin(\theta)\sin(k)\\\
i\sin(\theta)\sin(k)&\cos(\theta)-i\sin(\theta)\cos(k)\\\ \end{array}\right).$
(91)

Using expressions (68) and (72), one finds

$\begin{split}\sigma_{k}^{11}=\cos(\theta)+i\sin(\theta)\cos(k-\vartheta_{k}),\\\
\sigma_{k}^{22}=\cos(\theta)-i\sin(\theta)\cos(k-\vartheta_{k}).\end{split}$
(92)

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 reads
$\tilde{P}\left(\theta\right)=\frac{1}{Z(\beta,g,\gamma)}\left[\tilde{P}^{+}\left(\theta\right)+\tilde{P}^{-}\left(\theta\right)\right].$
(93) Positive part of characteristic function:
$\begin{split}\tilde{P}^{+}\left(\theta\right)=&\frac{\exp\left(iL\theta/2\right)}{2}\Bigg{[}\prod_{k\in\mathbf{k}^{+}}2\left(\cos(\theta)\cosh[\beta\epsilon_{k}(g,\gamma)]-i\sin(\theta)\sinh[\beta\epsilon_{k}(g,\gamma)]\cos(k-\vartheta_{k})+1\right)\\\
&+\prod_{k\in\mathbf{k}^{+}}2\left(\cos(\theta)\cosh[\beta\epsilon_{k}(g,\gamma)]-i\sin(\theta)\sinh[\beta\epsilon_{k}(g,\gamma)]\cos(k-\vartheta_{k})-1\right)\Bigg{]}.\end{split}$
(94) Negative part of characteristic function:
$\begin{split}\tilde{P}^{-}\left(\theta\right)=\frac{\exp\left(iL\theta/2\right)}{2}\Bigg{[}\tilde{P}^{F}(\theta)\prod_{k\in\mathbf{k}^{-}}2\left(\cos(\theta)\cosh[\beta\epsilon_{k}(g,\gamma)]-i\sin(\theta)\sinh[\beta\epsilon_{k}(g,\gamma)]\cos(k-\vartheta_{k})+1\right)\\\
-\tilde{P}^{B}(\theta)\prod_{k\in\mathbf{k}^{-}}2\left(\cos(\theta)\cosh[\beta\epsilon_{k}(g,\gamma)]-i\sin(\theta)\sinh[\beta\epsilon_{k}(g,\gamma)]\cos(k-\vartheta_{k})-1\right)\Bigg{]},\end{split}$
(95) where
$\begin{split}\tilde{P}^{F}(\theta)&=2^{2}\cosh\left(\frac{\beta\epsilon_{k=0}+i\theta}{2}\right)\cosh\left(\frac{\beta\epsilon_{k=\pi}-i\theta}{2}\right),\\\
\tilde{P}^{B}(\theta)&=2^{2}\sinh\left(\frac{\beta\epsilon_{k=0}+i\theta}{2}\right)\sinh\left(\frac{\beta\epsilon_{k=\pi}-i\theta}{2}\right).\end{split}$
(96) The exact total partition function is given by Eq. (48), with the
eigenenergies $\epsilon_{k}\left(g,\gamma\right)$ and $\epsilon_{k=0}$ given
by Eq. (47), and the Bogoliubov angles $\vartheta_{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 of
$\tilde{P}^{+}(\theta)$:

Summary 4.2:  PPA characteristic function for kink number In the thermodynamic
limit, Eq. (93) reduces to $\tilde{P}_{\rm
PPA}(\theta)=\frac{\exp\left(iL\theta/2\right)}{Z^{+}_{F}(\beta,g,\gamma)}\prod_{k\in\mathbf{k}^{+}}2\left(\cos(\theta)\cosh(\beta\epsilon_{k}(g,\gamma))-i\sin(\theta)\sinh(\beta\epsilon_{k}(g,\gamma))\cos(k-\vartheta_{k})+1\right),$
(97) where $Z_{F}^{+}(\beta,g,\gamma)$ is defined in (49).

Figure 3: Kink-number distribution at thermodynamic equilibrium. Probability
distribution of the number of kinks $P(n)$ as a function of the magnetic field
$g$ and temperature $T$ for a chain of $L=50$ spins. The exact probability
distribution Eq. (93) (red bars) is compared with the simplified expression in
Eq. (97) (blue bars). Only in the low-temperature paramagnet the PPA is
accurate. Further, the normal (Gaussian) approximation to the histograms is
also shown (dashed lines).

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\left(n\right)$ using Eq. (93). Additionally, we evaluate the PPA
probability distribution function using Eq. (97). 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 [58, 1, 59, 60, 62, 43]. 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:

$P_{\rm
PPA}\left(n=2\ell+1\right)=\frac{1}{2\pi}\int_{-\pi}^{\pi}d\theta\tilde{P}\left(\theta\right)\exp\left[-i\theta(2\ell+1)\right]\neq
0,$ (98)

but for large $g$ and $\beta$ as shown in 3, when $P_{\rm
PPA}\left(n=2\ell+1\right)\approx 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 $\prod_{n=1}^{L}\hat{X}_{n}\hat{X}_{n+1}$
which is $1$ for even kink number and $-1$ for an odd number. Using
$\hat{X}_{L+1}=\hat{X}_{1}$ and
$\left(\hat{X}_{n}\right)^{2}=\bigotimes_{n=1}^{L}\hat{\mathbb{I}}_{n}$, it
satisfies:

$\prod_{n=1}^{L}\hat{X}_{n}\hat{X}_{n+1}=1.$ (99)

The PPA characteristic function, $\tilde{P}^{+}(\theta)$ and
$\tilde{P}^{-}(\theta)$ 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\left(n\right)\approx P_{\rm PPA}\left(n\right)+\frac{1}{2}\left[P_{\rm
PPA}\left(n-1\right)+P_{\rm PPA}\left(n+1\right)\right],$ (100)

as shown in Figure 4.

Figure 4: Exact and Coarse-grained PPA kink-number probability distributions
at thermal equilibrium. The exact kink-number probability distribution
evaluated using Eq. (93) (red) is compared with the coarse-grained PPA
probability distribution Eq. (100) (blue). The numerical histograms are
compared with the Gaussian $N(\kappa_{1},\kappa_{2})$ with fitted numerical
values for $\kappa_{1}$ and $\kappa_{2}$ (dashed lines). In as much as the
exact distribution is symmetric and its left tail is negligible near the
origin, the coarse-graining of the PPA distribution in Eq. (100) reproduces
accurately the exact distribution. Deviations are manifested at low $g$ and
temperature, when the distribution is asymmetric.

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.

Figure 5: Cumulants of the kink-number distribution as a function of the
inverse of temperature $\beta$. Using the exact characteristic function given
by Eq. (93), the mean kink number $\kappa_{1}$ and the variance $\kappa_{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. (97). 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$.

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.

Figure 6: Relative error for the first two cumulants of the kink-number
distribution as a function of magnetic field $g$. Using the full
characteristic function in Eq. (93) and the PPA characteristic function Eq.
(97), the relative error is evaluated as a function of the magnetic field for
a system size $L=12$ and different temperatures.

To complete the characterization of the kink-number distribution we consider
the limiting cases of the ground-state distribution ($\beta\rightarrow\infty$)
and the infinite-temperature case ($\beta\rightarrow 0$) in an exact approach,
without using the PPA. The first can be easily described using (91), while in
the second we consider a maximally-mixed Gibbs state and apply trace formulas
2.1. For $\beta=0$, the exact result and the PPA coincide.

Summary 4.3:  Limiting cases of kink number distribution Exact ground-state
characteristic function of the kink-number distribution:
$\tilde{P}_{\beta\rightarrow\infty}(\theta)=\exp(iL\theta/2)\prod_{k\in\mathbf{k}^{+}}(\cos\theta-i\sin\theta\cos(k-\vartheta_{k})).$
(101) Exact infinite-temperature characteristic function of the kink-number
distribution: $\tilde{P}_{\beta\rightarrow
0}(\theta)=\exp(iL\theta/2)\,\left(\cos^{L}\frac{\theta}{2}+(-1)^{L/2}\sin^{L}\frac{\theta}{2}\right).$
(102)

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.
(102).

Figure 7: Limiting cases of kink number distribution. Probability distribution
of the number of kinks $P(n)$ as a function of the magnetic field $g$ and
inverse temperature $\beta$ for a chain of $L=50$ spins. The left panel shows
the kink-number distribution for different values of the magnetic field and is
obtained using the ground-state characteristic function Eq. (101). The right
panel shows the kink number distribution at infinity temperature, computed
using the characterization function given by Eq. (102). The vertical dashed
red line is located at $\kappa_{1}=L/2$, while the long-dashed black line
corresponds to the Gaussian approximation $N(L/2,L/4)$.

In fact, full probability distribution for infinite temperature can be found
by a combinatorial argument. Working in the basis of eigenstates of
$\sigma^{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\binom{L}{2l}$ ways.
Therefore, the full probability distribution has the form:

$P_{\beta\rightarrow 0}(n=2l)=\frac{1}{2^{L-1}}\,\binom{L}{2l},\quad
l=0,1,\ldots\frac{L}{2}.$ (103)

The corresponding cumulant values read

$\kappa_{1}=\frac{L}{2},\quad\kappa_{2}=\frac{L}{4},\quad\kappa_{3}=0,\quad\kappa_{4}=-\frac{L}{8},\quad\kappa_{5}=0,\quad\kappa_{6}=\frac{L}{4},\quad\dots$
(104)

By keeping the first two cumulants and setting the rest to zero,
$P_{\beta\rightarrow 0}(n=2l)$ can be approximated by a Gaussian distribution
$N(\kappa_{1},\kappa_{2})$ with mean $\kappa_{1}=L/2$ and variance
$\kappa_{2}=L/4$. As shown in Fig. 7 this approximation describes the envelope
of the distribution with great accuracy.

###  4.2 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 magnetization

$\hat{M}=\sum_{n=1}^{L}\hat{Z}_{n},$ (105)

with eigenvalues $m=-L,-L+2,\dots,L-2,L$ for even $L$. The latter has been
studied in the PPA and continuous approximations and finds broad applications
in the characterization of quantum critical behavior [80, 81, 82, 83, 84, 86,
87] and the identification of various many-body states in ultracold-atom
quantum simulators [85].

In the Fourier representation, it is the sum of two different contributions:

$\hat{M}^{+}=\sum_{k\in\mathbf{k}^{+}}2(\hat{c}_{k}\hat{c}_{k}^{\dagger}-\hat{c}_{k}^{\dagger}\hat{c}_{k}),\quad\hat{M}^{-}=\sum_{k\in\mathbf{k}^{-}}2(\hat{c}_{k}\hat{c}_{k}^{\dagger}-\hat{c}_{k}^{\dagger}\hat{c}_{k})+\hat{c}_{0}\hat{c}_{0}^{\dagger}-\hat{c}_{0}^{\dagger}\hat{c}_{0}+\hat{c}_{\pi}\hat{c}_{\pi}^{\dagger}-\hat{c}_{\pi}^{\dagger}\hat{c}_{\pi}.$
(106)

In parallel with Eq. (89), we define a new set of a single-mode operators
$\hat{m}_{k}$, $\hat{m}_{0}$, and $\hat{m}_{\pi}$,

$\displaystyle\hat{m}_{k}$
$\displaystyle=\left(\begin{array}[]{cccc}2&0&0&0\\\ 0&-2&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ \end{array}\right),$ $\displaystyle\hat{m}_{k}^{(1)}$
$\displaystyle=\left(\begin{array}[]{cc}2&0\\\ 0&-2\\\ \end{array}\right),$
$\displaystyle\hat{m}_{k}^{(2)}$ $\displaystyle=0_{2}.$ (113)

In addition, in the negative-parity sector, the matrix $\hat{m}_{k}$ has the
same form for the momenta $0,\pi$ that is given by
$\hat{m}_{0}=\hat{m}_{\pi}={\rm diag}\left(1,-1\right)$. We can easily compute
$\exp\left(i\theta\hat{m}_{k}^{(1)}\right)$ and the $\hat{\sigma}_{k}$ matrix
to obtain

$\begin{split}\hat{\sigma}_{k}^{(11)}=\cos(2\theta)+i\cos(\vartheta_{k})\sin(2\theta),\\\
\hat{\sigma}_{k}^{(22)}=\cos(2\theta)-i\cos(\vartheta_{k})\sin(2\theta).\end{split}$
(114)

Summary 4.4:  Full generating function of transverse magnetization The full
characteristic function for the transverse magnetization Eq. (105) at thermal
equilibrium reads
$\tilde{P}\left(\theta\right)=\frac{1}{Z(\beta,g,\gamma)}\left(\tilde{P}^{+}\left(\theta\right)+\tilde{P}^{-}\left(\theta\right)\right).$
(115) Positive part of characteristic function:
$\begin{split}\tilde{P}^{+}\left(\theta\right)=&\frac{1}{2}\Bigg{[}\prod_{k\in\mathbf{k}^{+}}2\left(\cos(2\theta)\cosh(\beta\epsilon_{k}(g,\gamma))-i\sin(2\theta)\sinh(\beta\epsilon_{k}(g,\gamma))\cos(\vartheta_{k})+1\right)\\\
&+\prod_{k\in\mathbf{k}^{+}}2\left(\cos(2\theta)\cosh(\beta\epsilon_{k}(g,\gamma))-i\sin(2\theta)\sinh(\beta\epsilon_{k}(g,\gamma))\cos(\vartheta_{k})-1\right)\Bigg{]}.\end{split}$
(116) Negative part of characteristic function:
$\begin{split}\tilde{P}^{-}\left(\theta\right)=\frac{1}{2}\Bigg{[}\tilde{P}^{F}(\theta)\prod_{k\in\mathbf{k}^{-}}2\left(\cos(2\theta)\cosh(\beta\epsilon_{k}(g,\gamma))-i\sin(2\theta)\sinh(\beta\epsilon_{k}(g,\gamma))\cos(\vartheta_{k})+1\right)\\\
-\tilde{P}^{B}(\theta)\prod_{k\in\mathbf{k}^{-}}2\left(\cos(2\theta)\cosh(\beta\epsilon_{k}(g,\gamma))-i\sin(2\theta)\sinh(\beta\epsilon_{k}(g,\gamma))\cos(\vartheta_{k})-1\right)\Bigg{]},\end{split}$
(117) with
$\begin{split}\tilde{P}^{F}(\theta)&=2^{2}\cosh\left(\frac{\beta\epsilon_{k=0}+2i\theta}{2}\right)\cosh\left(\frac{\beta\epsilon_{k=\pi}+2i\theta}{2}\right),\\\
\tilde{P}^{B}(\theta)&=2^{2}\sinh\left(\frac{\beta\epsilon_{k=0}+2i\theta}{2}\right)\sinh\left(\frac{\beta\epsilon_{k=\pi}+2i\theta}{2}\right).\end{split}$
(118) The exact partition function is given by Eq. (48), with the
eigenenergies $\epsilon_{k}\left(g,\gamma\right)$ and $\epsilon_{k=0}$ given
by Eq. (47), and the Bogoliubov angles $\vartheta_{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 of
$\tilde{P}^{+}(\theta)$:

Summary 4.5:  PPA characteristic function for transverse magnetization In the
thermodynamic limit, Eq. (115) reduces to $\tilde{P}_{\rm
PPA}(\theta)=\frac{1}{2Z_{F}^{+}(\beta,g,\gamma)}\prod_{k\in\mathbf{k}^{+}}2\left(\cos(2\theta)\cosh(\beta\epsilon_{k}(g,\gamma))-i\sin(2\theta)\sinh(\beta\epsilon_{k}(g,\gamma))\cos(\vartheta_{k})+1\right),$
(119) where $Z_{F}^{+}(\beta,g,\gamma)$ is defined in (49).

Figure 8: Magnetization distribution at thermodynamic equilibrium. Probability
distribution of the transverse magnetization $P(m)$ for different values of
the magnetic field $g$ and inverse temperature $\beta$ in a chain of $L=50$
spins. The exact probability distribution Eq. (115) (red bars) is compared
with the simplified expression in Eq. (119) (blue bars).

The magnetization distribution is shown in Fig. 8 for different values of $g$
and $\beta$ 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 $\beta$ with discrepancies being
noticeable in the pure ferromagnet ($g=0$) at low temperature. As the magnetic
field is cranked up at constant $\beta$, 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 $\beta$ as shown in Fig. 10.

Figure 9: Cumulants of the magnetization distribution as a function of the
inverse of temperature $\beta$. Using the full characteristic function given
by Eq. (115), the mean value of the transversal magnetization $\kappa_{1}$ and
the variance $\kappa_{2}$ are shown by red circles and blue squares,
respectively. The dashed lines correspond to the numerical results using the
simplified characteristic function (Eq. (119)). In the ferromagnetic phase the
statistics is super-Poissonian, while it changes from super- to sub-Poissonian
in the paramagnetic phase as the temperature is decreased. The magnetic field
varies from $0.0$ to $2.0$ from left to right in steps of $0.5$. The system
size is $L=12$ in the upper row and $L=100$ in the lower one.  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. (115) and the corresponding PPA Eq. (119), the relative error is evaluated
as a function of the magnetic field for a system size $L=12$ and different
temperatures.

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 $\beta$.

Summary 4.6:  Limiting cases of transverse magnetization distribution Exact
ground-state characteristic function of transverse magnetization:
$\tilde{P}_{\beta\rightarrow\infty}(\theta)=\prod_{k\in\mathbf{k}^{+}}(\cos
2\theta-i\sin 2\theta\cos\vartheta_{k}).$ (120) Exact infinite-temperature
characteristic function of transverse magnetization:
$\tilde{P}_{\beta\rightarrow 0}(\theta)=\cos^{L}\theta.$ (121)

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 (120) and (63), one can find the first
cumulants of the ground-state distribution explicitly. In particular, the
first few cumulants read

$\displaystyle\kappa_{1}$ $\displaystyle=$
$\displaystyle-\sum_{k\in\mathbf{k}^{+}}2\cos\vartheta_{k},$ (122)
$\displaystyle\kappa_{2}$ $\displaystyle=$ $\displaystyle
L-2\sum_{k\in\mathbf{k}^{+}}\cos(2\vartheta_{k}),$ (123)
$\displaystyle\kappa_{3}$ $\displaystyle=$ $\displaystyle
4\sum_{k\in\mathbf{k}^{+}}\left[\cos(\vartheta_{k})-\cos(3\vartheta_{k})\right].$
(124)

The second cumulant turns out to have a particularly simple form due to its
close relation to the ground-state fidelity susceptibility [90, 72] and reads

$\kappa_{2}=L\,\frac{1+g^{L-2}}{1+g^{L}}.$ (125)

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 $\pm 1$ with equal
probability $1/2$. As a result $\kappa_{1}=0$, $\kappa_{2}=L/4$. In the
infinite temperature limit, $P(m)$ is equal to that of a classical Ising chain
and can be written explicitly:

$P_{\beta\rightarrow
0}(m)=\frac{1}{2^{L}}\,\binom{L}{\frac{1}{2}\left(m+L\right)},\quad
m=-L,-L+2,\ldots,L-2,L.$ (126)

Odd cumulant identically vanish, while the first even ones read

$\kappa_{2}=L,\quad\kappa_{4}=-2L,\quad\kappa_{6}=16L,\quad\kappa_{8}=-272L,\quad\kappa_{10}=7936L,\quad\dots$
(127)

As a result, in the normal approximation $P_{\beta\rightarrow 0}(n=2l)$ is
given by Gaussian distribution with zero mean and variance $\kappa_{2}=L$.
Fig. 11 shows this Gaussian distribution as a black envelope, accurately
approximating the exact results.

Figure 11: Limiting cases of transverse magnetization distribution.
Probability distribution of the transverse magnetization $P(m)$ as a function
of the magnetic field $g$ and inverse temperature $\beta$ for a chain of
$L=50$ spins. The left panel shows the ground-state transverse magnetization
distribution for different values of the magnetic field, and is computed using
the characteristic function Eq. (120). The right panel shows the transverse
magnetization distribution at infinity temperature, obtained using the
characterization function given by Eq. (121). The envelope of the distribution
is reproduced by the Gaussian approximation $N(0,L)$ shown as a dashed black
line.

##  5 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 one-dimensional 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.

## Acknowledgements

#### Funding information

We acknowledge funding support from the Spanish Ministerio de Ciencia e
Innovación (PID2019-109007GA-I00). MB acknowledges support of Polish National
Science Center (NCN) grant DEC-2016/23/B/ST3/01152 and Faculty of Physics,
astronomy and Computer Science internal grant N17/MNS/000008.

##  Appendix 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\geq 1$ and we compute

$\begin{split}\mathrm{tr}\left[\mathcal{P}\left(\bigotimes_{i=1}^{n+1}\hat{O}_{k_{i}}\right)\right]&=\mathrm{tr}\left[\mathcal{P}\left(\bigotimes_{i=1}^{n}\hat{O}_{k_{i}}\right)\right]\cdot\mathrm{tr}\left[\hat{O}_{k_{n+1}}^{(p)}\right]+\mathrm{tr}\left[\mathcal{N}\left(\bigotimes_{i=1}^{n}O_{k_{i}}\right)\right]\cdot\mathrm{tr}\left[\hat{O}_{k_{n+1}}^{(n)}\right]\\\
&=\frac{1}{2}\Bigg{[}\left(\mathrm{tr}\left[\hat{O}_{k_{n+1}}^{(p)}\right]+\mathrm{tr}\left[\hat{O}_{k_{n+1}}^{(n)}\right]\right)\,\prod_{i=1}^{n}\mathrm{tr}\left[\hat{O}_{k_{i}}\right]\\\
&\;\;\quad+\left(\mathrm{tr}\left[\hat{O}_{k_{n+1}}^{(p)}\right]-\mathrm{tr}\left[\hat{O}_{k_{n+1}}^{(n)}\right]\right)\,\prod_{i=1}^{n}\left(\mathrm{tr}\left[\hat{O}_{k_{i}}^{(p)}\right]-\mathrm{tr}\left[\hat{O}_{k_{i}}^{(n)}\right]\right)\Bigg{]}\\\
&=\frac{1}{2}\left[\prod_{i=1}^{n+1}\mathrm{tr}\left[\hat{O}_{k_{i}}\right]+\prod_{i=1}^{n+1}\left(\mathrm{tr}\left[\hat{O}_{k_{i}}^{(p)}\right]-\mathrm{tr}\left[\hat{O}_{k_{i}}^{(n)}\right]\right)\right],\end{split}$

and an inductive step is completed.

## References

  * [1] S. Sachdev,  _Quantum Phase Transitions, Second Edition_ ,  Cambridge University Press, Cambridge (2011).
  * [2] B. K. Chakrabarti, A. Dutta and P. Sen,  _Quantum Ising Phases and Transitions in Transverse Ising Models_ , vol. m41,  Springer-Verlag, Berlin (1996).
  * [3] A. Dutta, G. Aeppli, B. K. Chakrabarti, U. Divakaran, T. F. Rosenbaum and D. Sen,  _Quantum Phase Transitions in Transverse Field Spin Models: From Statistical Physics to Quantum Information_ ,  Cambridge University Press,  10.1017/CBO9781107706057 (2015).
  * [4] E. Lieb, T. Schultz and D. Mattis,  _Two soluble models of an antiferromagnetic chain_ ,  Annals of Physics 16(3), 407 (1961),  https://doi.org/10.1016/0003-4916(61)90115-4.
  * [5] P. Pfeuty,  _The one-dimensional Ising model with a transverse field_ ,  Annals of Physics 57(1), 79 (1970),  https://doi.org/10.1016/0003-4916(70)90270-8.
  * [6] A. Kitaev,  _Anyons in an exactly solved model and beyond_ ,  Annals of Physics 321(1), 2 (2006),  https://doi.org/10.1016/j.aop.2005.10.005,  January Special Issue.
  * [7] S. Sachdev,  _Universal, finite-temperature, crossover functions of the quantum transition in the Ising chain in a transverse field_ ,  Nuclear Physics B 464(3), 576 (1996),  https://doi.org/10.1016/0550-3213(95)00657-5.
  * [8] A. Leclair, F. Lesage, S. Sachdev and H. Saleur,  _Finite temperature correlations in the one-dimensional quantum Ising model_ ,  Nuclear Physics B 482(3), 579 (1996),  https://doi.org/10.1016/S0550-3213(96)00456-7.
  * [9] S. Sachdev and A. P. Young,  _Low temperature relaxational dynamics of the Ising chain in a transverse field_ ,  Phys. Rev. Lett. 78, 2220 (1997),  10.1103/PhysRevLett.78.2220.
  * [10] A. Polkovnikov,  _Universal adiabatic dynamics in the vicinity of a quantum critical point_ ,  Phys. Rev. B 72, 161201 (2005),  10.1103/PhysRevB.72.161201.
  * [11] J. Dziarmaga,  _Dynamics of a quantum phase transition: Exact solution of the quantum Ising model_ ,  Phys. Rev. Lett. 95, 245701 (2005),  10.1103/PhysRevLett.95.245701.
  * [12] M. Gong, X. Wen, G. Sun, D.-W. Zhang, D. Lan, Y. Zhou, Y. Fan, Y. Liu, X. Tan, H. Yu, Y. Yu, S.-L. Zhu _et al._ ,  _Simulating the Kibble-Zurek mechanism of the Ising model with a superconducting qubit system_ ,  Scientific Reports 6(1), 22667 (2016),  10.1038/srep22667.
  * [13] J.-M. Cui, Y.-F. Huang, Z. Wang, D.-Y. Cao, J. Wang, W.-M. Lv, L. Luo, A. del Campo, Y.-J. Han, C.-F. Li and G.-C. Guo,  _Experimental trapped-ion quantum simulation of the Kibble-Zurek dynamics in momentum space_ ,  Sci. Rep. 6, 33381 (2016),  10.1038/srep33381.
  * [14] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres _et al._ ,  _Quantum Kibble-Zurek mechanism and critical dynamics on a programmable Rydberg simulator_ ,  Nature 568(7751), 207 (2019),  10.1038/s41586-019-1070-1.
  * [15] Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J. Gómez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo and H. Nishimori,  _Probing the universality of topological defect formation in a quantum annealer: Kibble-Zurek mechanism and beyond_ ,  Phys. Rev. Research 2, 033369 (2020),  10.1103/PhysRevResearch.2.033369.
  * [16] D. Sen, K. Sengupta and S. Mondal,  _Defect production in nonlinear quench across a quantum critical point_ ,  Phys. Rev. Lett. 101, 016806 (2008),  10.1103/PhysRevLett.101.016806.
  * [17] R. Barankov and A. Polkovnikov,  _Optimal nonlinear passage through a quantum critical point_ ,  Phys. Rev. Lett. 101, 076801 (2008),  10.1103/PhysRevLett.101.076801.
  * [18] A. del Campo, T. W. B. Kibble and W. H. Zurek,  _Causality and non-equilibrium second-order phase transitions in inhomogeneous systems_ ,  J. Phys.Cond. Mat. 25(40), 404210 (2013).
  * [19] J. Dziarmaga and M. M. Rams,  _Dynamics of an inhomogeneous quantum phase transition_ ,  New J. Phys. 12(5), 055007 (2010).
  * [20] J. Dziarmaga and M. M. Rams,  _Adiabatic dynamics of an inhomogeneous quantum phase transition: the case of a $z>1$ dynamical exponent_,  New J. Phys. 12(10), 103002 (2010).
  * [21] M. M. Rams, M. Mohseni and A. del Campo,  _Inhomogeneous quasi-adiabatic driving of quantum critical dynamics in weakly disordered spin chains_ ,  New J. Phys. 18, 123034 (2016),  10.1088/1367-2630/aa5079.
  * [22] F. J. Gómez-Ruiz and A. del Campo,  _Universal dynamics of inhomogeneous quantum phase transitions: Suppressing defect formation_ ,  Phys. Rev. Lett. 122, 080604 (2019),  10.1103/PhysRevLett.122.080604.
  * [23] D. Patanè, A. Silva, L. Amico, R. Fazio and G. E. Santoro,  _Adiabatic dynamics in open quantum critical many-body systems_ ,  Phys. Rev. Lett. 101, 175701 (2008),  10.1103/PhysRevLett.101.175701.
  * [24] A. Dutta, A. Rahmani and A. del Campo,  _Anti-Kibble-Zurek behavior in crossing the quantum critical point of a thermally isolated system driven by a noisy control field_ ,  Phys. Rev. Lett. 117, 080402 (2016),  10.1103/PhysRevLett.117.080402.
  * [25] M.-Z. Ai, J.-M. Cui, R. He, Z.-H. Qian, X.-X. Gao, Y.-F. Huang, C.-F. Li and G.-C. Guo,  _Experimental verification of anti–Kibble-Zurek behavior in a quantum system under a noisy control field_ ,  Phys. Rev. A 103, 012608 (2021),  10.1103/PhysRevA.103.012608.
  * [26] L. Cincio, J. Dziarmaga, M. M. Rams and W. H. Zurek,  _Entropy of entanglement and correlations induced by a quench: Dynamics of a quantum phase transition in the quantum Ising model_ ,  Phys. Rev. A 75, 052321 (2007),  10.1103/PhysRevA.75.052321.
  * [27] A. del Campo,  _Universal statistics of topological defects formed in a quantum phase transition_ ,  Phys. Rev. Lett. 121, 200601 (2018),  10.1103/PhysRevLett.121.200601.
  * [28] J.-M. Cui, F. J. Gómez-Ruiz, Y.-F. Huang, C.-F. Li, G.-C. Guo and A. del Campo,  _Experimentally testing quantum critical dynamics beyond the Kibble-Zurek mechanism_ ,  Comm. Phys. 3, 44 (2020),  10.1038/s42005-020-0306-6.
  * [29] F. J. Gómez-Ruiz, J. J. Mayo and A. del Campo,  _Full counting statistics of topological defects after crossing a phase transition_ ,  Phys. Rev. Lett. 124, 240602 (2020),  10.1103/PhysRevLett.124.240602.
  * [30] A. del Campo, F. J. Gómez-Ruiz, Z.-H. Li, C.-Y. Xia, H.-B. Zeng and H.-Q. Zhang,  _Universal statistics of vortices in a newborn holographic superconductor: Beyond the Kibble-Zurek mechanism_ (2021), 2101.02171.
  * [31] A. del Campo, M. M. Rams and W. H. Zurek,  _Assisted finite-rate adiabatic passage across a quantum critical point: Exact solution for the quantum Ising model_ ,  Phys. Rev. Lett. 109, 115703 (2012),  10.1103/PhysRevLett.109.115703.
  * [32] K. Takahashi,  _Transitionless quantum driving for spin systems_ ,  Phys. Rev. E 87, 062117 (2013),  10.1103/PhysRevE.87.062117.
  * [33] H. Saberi, T. Opatrný, K. Mølmer and A. del Campo,  _Adiabatic tracking of quantum many-body dynamics_ ,  Phys. Rev. A 90, 060301 (2014),  10.1103/PhysRevA.90.060301.
  * [34] J. D. Sau and K. Sengupta,  _Suppressing defect production during passage through a quantum critical point_ ,  Phys. Rev. B 90, 104306 (2014),  10.1103/PhysRevB.90.104306.
  * [35] N. Wu, A. Nanduri and H. Rabitz,  _Optimal suppression of defect generation during a passage across a quantum critical point_ ,  Phys. Rev. B 91, 041115 (2015),  10.1103/PhysRevB.91.041115.
  * [36] A. del Campo and K. Sengupta,  _Controlling quantum critical dynamics of isolated systems_ ,  Euro. Phys. J. Spe. Top. 224(1), 189 (2015).
  * [37] P. W. Claeys, M. Pandey, D. Sels and A. Polkovnikov,  _Floquet-engineering counterdiabatic protocols in quantum many-body systems_ ,  Phys. Rev. Lett. 123, 090602 (2019),  10.1103/PhysRevLett.123.090602.
  * [38] M. M. Rams, J. Dziarmaga and W. H. Zurek,  _Symmetry breaking bias and the dynamics of a quantum phase transition_ ,  Phys. Rev. Lett. 123, 130603 (2019),  10.1103/PhysRevLett.123.130603.
  * [39] G. Passarelli, V. Cataudella, R. Fazio and P. Lucignano,  _Counterdiabatic driving in the quantum annealing of the $p$-spin model: A variational approach_,  Phys. Rev. Research 2, 013283 (2020),  10.1103/PhysRevResearch.2.013283.
  * [40] K. Suzuki and K. Takahashi,  _Performance evaluation of adiabatic quantum computation via quantum speed limits and possible applications to many-body systems_ ,  Phys. Rev. Research 2, 032016 (2020),  10.1103/PhysRevResearch.2.032016.
  * [41] L. Prielinger, A. Hartmann, Y. Yamashiro, K. Nishimura, W. Lechner and H. Nishimori,  _Diabatic quantum annealing by counter-diabatic driving_ (2020), 2011.02691.
  * [42] A. Silva,  _Statistics of the work done on a quantum critical system by quenching a control parameter_ ,  Phys. Rev. Lett. 101, 120603 (2008),  10.1103/PhysRevLett.101.120603.
  * [43] R. Dorner, J. Goold, C. Cormick, M. Paternostro and V. Vedral,  _Emergent thermodynamics in a quenched quantum many-body system_ ,  Phys. Rev. Lett. 109, 160601 (2012),  10.1103/PhysRevLett.109.160601.
  * [44] Z. Fei and H. T. Quan,  _Group-theoretical approach to the calculation of quantum work distribution_ ,  Phys. Rev. Research 1, 033175 (2019),  10.1103/PhysRevResearch.1.033175.
  * [45] Z. Fei, N. Freitas, V. Cavina, H. T. Quan and M. Esposito,  _Work statistics across a quantum phase transition_ ,  Phys. Rev. Lett. 124, 170603 (2020),  10.1103/PhysRevLett.124.170603.
  * [46] R. B. S, V. Mukherjee, U. Divakaran and A. del Campo,  _Universal finite-time thermodynamics of many-body quantum machines from Kibble-Zurek scaling_ ,  Phys. Rev. Research 2, 043247 (2020),  10.1103/PhysRevResearch.2.043247.
  * [47] R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl and K. Kiefer,  _Quantum criticality in an Ising chain: Experimental evidence for emergent $E_{8}$ symmetry_,  Science 327(5962), 177 (2010),  10.1126/science.1180085.
  * [48] D. Porras and J. I. Cirac,  _Effective quantum spin systems with trapped ions_ ,  Phys. Rev. Lett. 92, 207901 (2004),  10.1103/PhysRevLett.92.207901.
  * [49] A. Friedenauer, H. Schmitz, J. T. Glueckert, D. Porras and T. Schaetz,  _Simulating a quantum magnet with trapped ions_ ,  Nature Physics 4(10), 757 (2008),  10.1038/nphys1032.
  * [50] K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan and C. Monroe,  _Quantum simulation of frustrated Ising spins with trapped ions_ ,  Nature 465(7298), 590 (2010),  10.1038/nature09071.
  * [51] R. Islam, E. E. Edwards, K. Kim, S. Korenblit, C. Noh, H. Carmichael, G.-D. Lin, L.-M. Duan, C.-C. Joseph Wang, J. K. Freericks and C. Monroe,  _Onset of a quantum phase transition with a trapped ion quantum simulator_ ,  Nature Communications 2(1), 377 (2011),  10.1038/ncomms1374.
  * [52] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov and C. Monroe,  _Non-local propagation of correlations in quantum systems with long-range interactions_ ,  Nature 511(7508), 198 (2014),  10.1038/nature13450.
  * [53] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss and M. Greiner,  _Quantum simulation of antiferromagnetic spin chains in an optical lattice_ ,  Nature 472(7343), 307 (2011),  10.1038/nature09994.
  * [54] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud _et al._ ,  _Quantum annealing with manufactured spins_ ,  Nature 473(7346), 194 (2011),  10.1038/nature10012.
  * [55] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp and A. Wallraff,  _Digital quantum simulation of spin models with circuit quantum electrodynamics_ ,  Phys. Rev. X 5, 021027 (2015),  10.1103/PhysRevX.5.021027.
  * [56] R. Barends, A. Shabani, L. Lamata, J. Kelly, A. Mezzacapo, U. L. Heras, R. Babbush, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro _et al._ ,  _Digitized adiabatic quantum computing with a superconducting circuit_ ,  Nature 534(7606), 222 (2016),  10.1038/nature17658.
  * [57] A. Cervera-Lierta,  _Exact Ising model simulation on a quantum computer_ ,  Quantum 2, 114 (2018),  10.22331/q-2018-12-21-114.
  * [58] E. Barouch and B. M. McCoy,  _Statistical mechanics of the $XY$ model. II. spin-correlation functions_,  Phys. Rev. A 3, 786 (1971),  10.1103/PhysRevA.3.786.
  * [59] H. T. Quan and F. M. Cucchietti,  _Quantum fidelity and thermal phase transitions_ ,  Phys. Rev. E 79, 031101 (2009),  10.1103/PhysRevE.79.031101.
  * [60] X.-M. Lu, Z. Sun, X. Wang and P. Zanardi,  _Operator fidelity susceptibility, decoherence, and quantum criticality_ ,  Phys. Rev. A 78, 032309 (2008),  10.1103/PhysRevA.78.032309.
  * [61] A. Lamacraft and P. Fendley,  _Order parameter statistics in the critical quantum Ising chain_ ,  Phys. Rev. Lett. 100, 165706 (2008),  10.1103/PhysRevLett.100.165706.
  * [62] S. Deng, G. Ortiz and L. Viola,  _Dynamical critical scaling and effective thermalization in quantum quenches: Role of the initial state_ ,  Phys. Rev. B 83, 094304 (2011),  10.1103/PhysRevB.83.094304.
  * [63] S. Katsura,  _Statistical mechanics of the anisotropic linear Heisenberg model_ ,  Phys. Rev. 127, 1508 (1962),  10.1103/PhysRev.127.1508.
  * [64] V. S. Kapitonov and K. N. Il’inskii,  _Functional representations for correlators of spin chains_ ,  Journal of Mathematical Sciences 88(2), 233 (1998),  10.1007/BF02364984.
  * [65] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong and C. Monroe,  _Observation of a many-body dynamical phase transition with a $53$-qubit quantum simulator_,  Nature 551, 601 (2017),  10.1038/nature24654.
  * [66] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres _et al._ ,  _Quantum Kibble-Zurek mechanism and critical dynamics on a programmable Rydberg simulator_ ,  Nature 568(7751), 207 (2019),  10.1038/s41586-019-1070-1.
  * [67] P. Jordan and E. Wigner,  _Über das paulische äquivalenzverbot_ ,  Zeitschrift für Physik 47(9), 631 (1928),  10.1007/BF01331938.
  * [68] E. Barouch, B. M. McCoy and M. Dresden,  _Statistical mechanics of the $\mathrm{XY}$ model. I_,  Phys. Rev. A 2, 1075 (1970),  10.1103/PhysRevA.2.1075.
  * [69] E. Barouch and B. M. McCoy,  _Statistical mechanics of the $\mathrm{XY}$ model. III_,  Phys. Rev. A 3, 2137 (1971),  10.1103/PhysRevA.3.2137.
  * [70] B. M. McCoy, E. Barouch and D. B. Abraham,  _Statistical mechanics of the $\mathrm{XY}$ model. IV. time-dependent spin-correlation functions_,  Phys. Rev. A 4, 2331 (1971),  10.1103/PhysRevA.4.2331.
  * [71] C. Itzykson and J.-M. Drouffe,  _Statistical Field Theory_ , vol. 1 of _Cambridge Monographs on Mathematical Physics_ ,  Cambridge University Press,  10.1017/CBO9780511622779 (1989).
  * [72] B. Damski and M. M. Rams,  _Exact results for fidelity susceptibility of the quantum Ising model: the interplay between parity, system size, and magnetic field_ ,  Journal of Physics A: Mathematical and Theoretical 47(2), 025303 (2013),  10.1088/1751-8113/47/2/025303.
  * [73] M. Okuyama, Y. Yamanaka, H. Nishimori and M. M. Rams,  _Anomalous behavior of the energy gap in the one-dimensional quantum $XY$ model_,  Phys. Rev. E 92, 052116 (2015),  10.1103/PhysRevE.92.052116.
  * [74] Z. Xu and A. del Campo,  _Probing the full distribution of many-body observables by single-qubit interferometry_ ,  Phys. Rev. Lett. 122, 160602 (2019),  10.1103/PhysRevLett.122.160602.
  * [75] L. S. Levitov and G. B. Lesovik,  _Charge distribution in quantum shot noise_ ,  JETP Lett. 58, 230 (1993).
  * [76] L. S. Levitov, H. Lee and G. B. Lesovik,  _Electron counting statistics and coherent states of electric current_ ,  Journal of Mathematical Physics 37(10), 4845 (1996),  10.1063/1.531672.
  * [77] I. Klich,  _An Elementary Derivation of Levitov’s Formula_ , pp. 397–402,  Springer Netherlands, Dordrecht,  ISBN 978-94-010-0089-5,  10.1007/978-94-010-0089-5_19 (2003).
  * [78] F. Haake,  _Quantum Signatures of Chaos_ ,  Springer, Berlin (2010).
  * [79] F. Binder, L. Correa, C. Gogolin, J. Anders and G. Adesso,  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_ ,  Fundamental Theories of Physics. Springer International Publishing,  ISBN 9783319990460 (2019).
  * [80] K. Binder,  _Finite size scaling analysis of Ising model block distribution functions_ ,  Z. Physik B Cond. Mat. 43(2), 119 (1981),  10.1007/BF01293604.
  * [81] A. D. Bruce and N. B. Wilding,  _Scaling fields and universality of the liquid-gas critical point_ ,  Phys. Rev. Lett. 68, 193 (1992),  10.1103/PhysRevLett.68.193.
  * [82] V. Aji and N. Goldenfeld,  _Fluctuations in finite critical and turbulent systems_ ,  Phys. Rev. Lett. 86, 1007 (2001),  10.1103/PhysRevLett.86.1007.
  * [83] S. T. Bramwell, P. C. W. Holdsworth and J.-F. Pinton,  _Universality of rare fluctuations in turbulence and critical phenomena_ ,  Nature 396, 552 (1998),  10.1038/25083.
  * [84] V. Eisler, Z. Rácz and F. van Wijland,  _Magnetization distribution in the transverse Ising chain with energy flux_ ,  Phys. Rev. E 67, 056129 (2003),  10.1103/PhysRevE.67.056129.
  * [85] R. W. Cherng and E. Demler,  _Quantum noise analysis of spin systems realized with cold atoms_ ,  New Journal of Physics 9(1), 7 (2007),  10.1088/1367-2630/9/1/007.
  * [86] S. Groha, F. H. L. Essler and P. Calabrese,  _Full counting statistics in the transverse field Ising chain_ ,  SciPost Phys. 4, 43 (2018),  10.21468/SciPostPhys.4.6.043.
  * [87] N. Wu,  _Exact one- and two-site reduced dynamics in a finite-size quantum Ising ring after a quench: a semi-analytical approach_ (2021), 2103.12509.
  * [88] M. Białończyk and B. Damski,  _Dynamics of longitudinal magnetization in transverse-field quantum Ising model: from symmetry-breaking gap to Kibble-Zurek mechanism_ ,  Journal of Statistical Mechanics: Theory and Experiment 2020(1), 013108 (2020),  10.1088/1742-5468/ab609a.
  * [89] N. Wu,  _Longitudinal magnetiation dynamics in the quantum Ising ring: A Pfaffian method based on correspondence between momentum space and real space_ ,  Phys. Rev. E 101, 042108 (2020),  10.1103/PhysRevE.101.042108.
  * [90] S.-J. GU,  _Fidelity approach to quantum phase transitions_ ,  International Journal of Modern Physics B 24(23), 4371 (2010),  10.1142/S0217979210056335.