Uhlmann fidelity and fidelity susceptibility for integrable spin chains at finite temperature: exact results

Quantum phase transitions (QPTs), unlike classical phase transitions, are induced by quantum fluctuations related to Heisenberg uncertainty and can occur at zero temperature. They are usually triggered by a change of a macroscopic parameter, like an external magnetic field, when two competing parts of the system Hamiltonian (represented by noncommuting operator terms) exchange magnitude. The crossing of a phase transition is signaled by a profound change in ground-state properties that are reflected in the order parameter (e.g. the magnetization in a ferromagnet) and correlation functions [1–3]. Quantum information theory provides a plethora of theoretical tools to analyze these properties and acquire more insights into the mechanism of QPTs [4]. Here, we focus on the fidelity that for pure states is defined in terms of the overlap between the ground states corresponding to close values of the external parameter g that drives the phase transition:

$$\begin{equation}{F}_{0}\left({g}_{1}\vert {g}_{2}\right)=\vert \langle \psi ({g}_{1})\vert \psi ({g}_{2})\rangle \vert ,\end{equation} \tag{ 1 }$$

where |ψ(g)⟩ denotes the ground state of parameter-dependent Hamiltonian H(g).

The ground-state fidelity quantifies how sensitive the ground-state is to changes in the control parameter g. The potential of the ground-state fidelity as an indicator of phase transitions became apparent in the work of Quan et al [5], where it was observed that the decay of Loschmidt echo is enhanced by criticality [6]. Shortly after, Zanardi and Paunković [7] proposed that static fidelity can be a good indicator of the phase transition and demonstrated this in the one-dimensional transverse field XY model and the Dicke model. These results initiated a vivid research activity focused on the use of fidelity to characterize quantum critical systems [8]. The success of this approach lies in the ability of the fidelity to capture abrupt changes in the ground state caused by a small variation of the driving parameter, with no prior knowledge of the order parameter or the location of the critical point. This observation suggests identifying the leading nontrivial term in a Taylor expansion in the shift of the control parameter

$$\begin{equation}{F}_{0}\left(g\vert g+\delta \right)=\vert \langle \psi (g)\vert \psi (g+\delta )\rangle \vert =1-\frac{{\chi }_{0}(g)}{2}{\delta }^{2}+\dots .\end{equation} \tag{ 2 }$$

The importance of the quantity χ(g), known as the fidelity susceptibility, was first pointed out by You et al [9]. More generally, one can consider the Riemann metric on the ground-state manifold spanned by varying all driving parameters, with the diagonal elements of the Riemann tensor corresponding to fidelity susceptibility, and any singularity in them indicating a phase transition [10].

The role of fidelity has been investigated in many models, including integrable spin chains like the Ising model in a transverse field and the XY models (and more general integrable spin chains), the Lipkin–Meshkov–Glick model [11], and the Bose–Hubbard model [12], to name some representative examples. For generic models of correlated fermions, bosons, and spin systems, the use of Monte Carlo methods has been proposed [13, 14]. Likewise, the use of tensor networks has also been advanced to compute the fidelity susceptibility, in the study of many-body quantum metrology [15]. The family of integrable spin chains is particularly useful in this context given its representation in terms of non-interacting fermions which allows for efficient computations and plenty of analytical results. This family is a good testbed for theoretical concepts and it is believed to capture properties of quantum many-body systems with long-range interactions, that are realized in laboratory [16, 17]. A phase transition in these systems is induced by varying the external transverse magnetic field (say, in the in z direction) represented by the Hamiltonian term $-g{\sum }_{i=1}^{N}{\sigma }_{i}^{z}$, where N is the system size.

Efforts to characterize the ground-state fidelity have been extended to the fidelity between thermal states [18, 19]. The latter introduces the inverse temperature β = 1/kT as an external control parameter, in addition to the driving parameter g, and requires generalizing equation (1) to mixed states using the Uhlmann fidelity [20–22]

$$\begin{equation}F\left({\beta }_{1},{g}_{1}\vert {\beta }_{2},{g}_{2}\right)=F({\rho }_{1},{\sigma }_{2})=\mathrm{Tr}\left(\sqrt{\sqrt{{\rho }_{1}}{\sigma }_{2}\sqrt{{\rho }_{1}}}\right),\end{equation} \tag{ 3 }$$

where ${\rho }_{1}\left({\beta }_{1},{g}_{1}\right)={\mathrm{e}}^{-{\beta }_{1}H({g}_{1})}/\mathrm{Tr}({\mathrm{e}}^{-{\beta }_{1}H({g}_{1})})$ and ${\sigma }_{2}\left({\beta }_{2},{g}_{2}\right)={\mathrm{e}}^{-{\beta }_{2}H({g}_{2})}/\mathrm{Tr}({\mathrm{e}}^{-{\beta }_{2}H({g}_{2})})$. One can define two kinds of susceptibilities, one with respect to a change of temperature with fixed value of the driving parameter, and the second one, with respect to a change of the parameter g at constant temperature:

$$\begin{equation}\xi (\beta ,g)={\left.-\frac{{\mathrm{d}}^{2}}{\mathrm{d}{\delta }^{2}}F\left(\beta ,g\vert \beta +\delta ,g\right)\right\vert }_{\delta =0},\end{equation} \tag{ 4 }$$

$$\begin{equation}\chi (\beta ,g)={\left.-\frac{{\mathrm{d}}^{2}}{\mathrm{d}{\delta }^{2}}F\left(\beta ,g\vert \beta ,g+\delta \right)\right\vert }_{\delta =0}.\end{equation} \tag{ 5 }$$

It was shown [19] that the thermal fidelity susceptibility ξ(β, g) is proportional to the specific heat. Further, for a sufficiently large temperature, the fidelity susceptibility related to the field χ(β, g) is proportional to magnetic susceptibility, as in this limit density matrices approximately commute and the phase transition can be treated as approximately classical.

On the other hand, one can consider the limit β → ∞, when one expects convergence to the results for ground-state fidelity. This limit justifies the expectation for the fidelity to remain a good indicator of phase transitions at finite temperature.

In this work, we focus on the exact characterization of the fidelity and fidelity susceptibility of finite-size integrable spin systems, paying particular attention to the intermediate temperature regime. Specifically, we consider that the thermal energy is comparable with the energy gap (or elementary excitation) of the system. At the critical point (and its vicinity), the fidelity susceptibility exhibits a sharp peak. The height of this peak measures how different are the states in the two phases. We provide exact formulas for the fidelity between thermal states of integrable spin chains and analyze the dependence of the height of the peak on temperature. Many spin chains widely considered in the literature share a common property, namely their Hamiltonian commutes with the parity operator

$$\begin{equation}P=\prod\limits _{i=1}^{N}{\sigma }_{i}^{z}.\end{equation} \tag{ 6 }$$

The latter has eigenvalues ±1 and therefore, the energy spectrum splits into two parts of positive and negative parity. For instance, in the Ising model for an even number of spins the ground state is in the positive parity subspace, but dealing with thermal states requires careful treatment of eigenstates of both parities. The necessity of proper handling of parities in finite-temperature spin chains was first stated by Katsura [23]. Later, Kapitonov and Il'inskii [24] provided a derivation of full partition function using integrals over Grassmann variables. In our recent work [25], we provide an elementary method to derive exact expressions for partition function and characteristic function of a wide class of observables using only the structure of Hilbert space. The big discrepancy, even in the thermodynamic limit, between the fidelity susceptibility of lowest energy states in positive and negative parity subspaces, was first reported by Damski and Rams [26]. A thorough analysis of the scaling of the ground state fidelity susceptibility in the XY model was conducted in [27], while closed-form expressions for the fidelity susceptibility in the quantum Ising model in a transverse field were derived in [28]. However, to our best knowledge, an exact treatment of the fidelity between thermal states of integrable spin chains in the complete Hilbert spaces seems to be lacking at the time of writing. We aim at filling this gap by using methods developed in [25]. We provide exact expressions for fidelity between two arbitrary thermal states in the quantum XY models and analyze the temperature dependence of fidelity susceptibility in a paradigmatic test-bed for quantum critical phenomena, the quantum Ising model in a transverse field.

In the remaining part of the paper, we first review the diagonalization of the XY model with particular stress on the exact treatment of the parity subspaces in section 2. In section 3 we recall the methods developed in [25] and use them to derive the full expression for Uhlmann fidelity between arbitrary thermal states in section 4. We then focus on fidelity susceptibility computed at the critical point, analyze its temperature dependence and compare it with the simplified expression obtained considering only the positive parity contribution; see section 5. Although such simplification is often justified in the thermodynamic limit, we show that this is not the case for fidelity susceptibility. In particular, we show that the discrepancy between exact and simplified expressions persists when increasing the system size. We also introduce an accurate low-temperature approximation using two lowest-lying energy eigenstates, that we use to characterize the thermal susceptibility and specific heat. Our results are thus of direct relevance to the use of the fidelity and fidelity susceptibility to characterize quantum critical phenomena.

In this section, we outline the basic steps to diagonalize the XY model in one spatial dimension. The method we use can be applied as well in extended XY models [29], with special care of the parity of the ground state. The Hamiltonian of an XY chain is

$$\begin{equation}H(\gamma ,g)=-\sum\limits _{i=1}^{N}\left(\frac{1+\gamma }{2}{\sigma }_{i}^{x}{\sigma }_{i+1}^{x}+\frac{1-\gamma }{2}{\sigma }_{i}^{y}{\sigma }_{i+1}^{y}\right)-g\sum\limits _{i=1}^{N}{\sigma }_{i}^{z}.\end{equation} \tag{ 7 }$$

When γ = 1, the Hamiltonian (7) corresponds to the Ising model in a transverse magnetic field (TFQIM), while the limit γ = 0 describes the isotropic XY model. For the anisotropic case 0 < γ ⩽ 1 the model belongs to the Ising universality class, and its phase diagram is determined by the value of g. When g > 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. On the other hand, when 0 ⩽ g < 1 the ground state of the system corresponds to a ferromagnetic configuration with polarization along the xy plane. These phases are separated by a second-order QPT at the critical point g = 1. Finally, for the isotropic case γ = 0, a QPT is observed between the gapless phase (g < 1) and the ferromagnetic phase (g > 1).

In the following, we assume periodic boundary conditions, i.e. ${\sigma }_{N+1}^{\alpha }={\sigma }_{1}^{\alpha }$. The Hamiltonian (7) can be diagonalized with mapping onto noninteracting fermions using the Jordan–Wigner transformation [30]

$$\begin{equation}{\sigma }_{i}^{z}=1-2{c}_{i}^{{\dagger}}{c}_{i},\hspace{25.0pt}{\sigma }_{i}^{x}=({c}_{i}+{c}_{i}^{{\dagger}})\prod\limits _{j{< }i}(1-2{c}_{j}^{{\dagger}}{c}_{j}),\hspace{25.0pt}{\sigma }_{i}^{y}=i({c}_{i}-{c}_{i}^{{\dagger}})\prod\limits _{j{< }i}(1-2{c}_{j}^{{\dagger}}{c}_{j}),\end{equation} \tag{ 8 }$$

where the fermionic operators satisfy $\left\{{c}_{i},{c}_{j}\right\}=0$ and $\left\{{c}_{i},{c}_{j}^{{\dagger}}\right\}={\delta }_{ij}$. Using the fermionic representation in real space, the Hamiltonian takes the form

$$\begin{align}H\left(\gamma ,g\right)& =-\sum\limits _{i=1}^{N-1}\left[{c}_{i}^{{\dagger}}{c}_{i+1}+{c}_{i+1}^{{\dagger}}{c}_{i}+\gamma \left({c}_{i}^{{\dagger}}{c}_{i+1}^{{\dagger}}+{c}_{i+1}{c}_{i}\right)\right]\\ & \quad +P\left[{c}_{N}^{{\dagger}}{c}_{1}+{c}_{1}^{{\dagger}}{c}_{N}+\gamma \left({c}_{N}{c}_{1}+{c}_{1}{c}_{N}\right)\right]\\ & \quad -g\sum\limits _{i=1}^{N}\left(1-2{c}_{i}^{{\dagger}}{c}_{i}\right),\end{align} \tag{ 9 }$$

where $P={\prod }_{i=1}^{N}{\sigma }_{i}^{z}$ is the parity operator already mentioned in the introduction. The next step is to make use of the Fourier transform

$$\begin{equation}{c}_{j}=\frac{\mathrm{exp}\left(-i\pi /4\right)}{\sqrt{N}}\sum\limits _{k\in \mathbf{K}}{c}_{k}\enspace \mathrm{exp}\left(ikj\right).\end{equation} \tag{ 10 }$$

Because of the presence of the parity operator, Hamiltonian (9) is not in the form of noninteracting fermions yet. However, the operator P has eigenvalues ±1 and commutes with the Hamiltonian. Therefore, the Hamiltonian can be diagonalized separately in two sectors of the total Hilbert space, with P = 1 (even number of quasiparticles) and P = −1 (odd number of quasiparticles). The difficulty lies in the fact that the boundary conditions obeyed by the fermionic operators depend on the sector of the Hilbert space. This changes the set of momenta K relevant for the Fourier transform. In the following, we write down the most important results in both subspaces. We assume an even number of particles N; for details see for example [26, 31].

Positive parity subspace. The condition P = 1 implies c_N+1 = −c₁ and the set of momenta is the following:

$$\begin{equation}{\mathbf{K}}^{+}=\left\{{\pm}\frac{\pi }{N},{\pm}\frac{3\pi }{N},\dots ,{\pm}\left(\pi -\frac{\pi }{N}\right)\right\}={\mathbf{k}}^{+}\cup -{\mathbf{k}}^{+},\end{equation} \tag{ 11 }$$

where k⁺ denotes the set of only positive momenta: ${\mathbf{k}}^{+}=\left\{\frac{\pi }{N},\frac{3\pi }{N},\dots ,\left(\pi -\frac{\pi }{N}\right)\right\}$. The Hamiltonian takes the form

$$\begin{equation}H(\gamma ,g)=2\sum\limits _{k\in {\mathbf{k}}^{+}}\left[\left({c}_{k}^{{\dagger}}{c}_{k}-{c}_{-k}{c}_{-k}^{{\dagger}}\right)\left(g-\mathrm{cos}\enspace k\right)+\left({c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}+{c}_{k}{c}_{-k}\right)\gamma \enspace \mathrm{sin}\enspace k\right].\end{equation} \tag{ 12 }$$

After a Bogoliubov transformation, it can be written as

$$\begin{equation}{H}^{+}(\gamma ,g)=\sum\limits _{k\in {\mathbf{k}}^{+}}2{{\epsilon}}_{k}\left(\gamma ,g\right)\left({\gamma }_{k}^{{\dagger}}{\gamma }_{k}+{\gamma }_{-k}^{{\dagger}}{\gamma }_{-k}-1\right),\end{equation} \tag{ 13 }$$

where, for all k ∈ k⁺:

$$\begin{equation}{{\epsilon}}_{k}\left(\gamma ,g\right)=\sqrt{{(g-\mathrm{cos}\enspace k)}^{2}+{(\gamma \enspace \mathrm{sin}\enspace k)}^{2}},\end{equation} \tag{ 14a }$$

$$\begin{equation}{\gamma }_{k}=\mathrm{cos}\left(\frac{{\vartheta }_{k}}{2}\right){c}_{k}+\mathrm{sin}\left(\frac{{\vartheta }_{k}}{2}\right){c}_{-k}^{{\dagger}},\end{equation} \tag{ 14b }$$

$$\begin{equation}(\mathrm{sin}\enspace {\vartheta }_{k},\mathrm{cos}\enspace {\vartheta }_{k})=\left(\frac{\gamma \enspace \mathrm{sin}\enspace k}{{{\epsilon}}_{k}\left(\gamma ,g\right)},\frac{g-\mathrm{cos}\enspace k}{{{\epsilon}}_{k}\left(\gamma ,g\right)}\right).\end{equation} \tag{ 14c }$$

The lowest energy state in this subspace and its energy are respectively given by

$$\begin{equation}\vert \gamma ,{g\rangle }^{+}=\prod\limits _{k\in {\mathbf{k}}^{+}}\left[\mathrm{cos}\left(\frac{{\vartheta }_{k}}{2}\right)-\mathrm{sin}\left(\frac{{\vartheta }_{k}}{2}\right){c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}\right]\left\vert \mathrm{v}\mathrm{a}\mathrm{c}\right\rangle .\end{equation} \tag{ 15a }$$

$$\begin{equation}{E}_{0}^{+}\left(\gamma ,g\right)=-2\sum\limits _{k\in {\mathbf{k}}^{+}}{{\epsilon}}_{k}\left(\gamma ,g\right),\end{equation} \tag{ 15b }$$

where $\left\vert \mathrm{v}\mathrm{a}\mathrm{c}\right\rangle$ is a state annihilated by fermionic operators c_k for k ∈ K⁺.

Negative parity subspace. Similarly, the condition P = −1 implies c_N+1 = c₁ and the set of allowed momenta

$$\begin{equation}{\mathbf{K}}^{-}=\left\{0,{\pm}\frac{2\pi }{N},{\pm}\frac{4\pi }{N},\dots ,{\pm}\left(\pi -\frac{2\pi }{N}\right),\pi \right\}={\mathbf{k}}^{-}\cup -{\mathbf{k}}^{-}\cup \left\{0,\pi \right\},\end{equation} \tag{ 16 }$$

where the set of positive momenta different then 0, π is ${\mathbf{k}}^{-}=\left\{\frac{2\pi }{N},\frac{4\pi }{N},\dots \left(\pi -\frac{2\pi }{N}\right)\right\}$. All steps of diagonalization are the same as for the positive-parity subspace, except for the fact that the modes with 0 and π momenta require careful treatment. We thus repeat the steps for k ∈ k⁻ and treat 0, π momenta separately.

The Hamiltonian takes the form

$$\begin{align}{H}^{-}(\gamma ,g)& =2\sum\limits _{k\in {\mathbf{k}}^{-}}\left[\left({c}_{k}^{{\dagger}}{c}_{k}-{c}_{-k}{c}_{-k}^{{\dagger}}\right)\left(g-\mathrm{cos}\enspace k\right)\right.\\ & \quad +\left.\left({c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}+{c}_{k}{c}_{-k}\right)\gamma \enspace \mathrm{sin}\enspace k\right]\\ & \quad +(g-1)({c}_{0}^{{\dagger}}{c}_{0}-{c}_{0}{c}_{0}^{{\dagger}})+(g+1)({c}_{\pi }^{{\dagger}}{c}_{\pi }-{c}_{\pi }{c}_{\pi }^{{\dagger}}).\end{align} \tag{ 17 }$$

After a Bogoliubov transformation, it can be written as

$$\begin{align}{H}^{-}(\gamma ,g)& =\sum\limits _{k\in {\mathbf{k}}^{-}}2{{\epsilon}}_{k}\left(\gamma ,g\right)\left({\gamma }_{k}^{{\dagger}}{\gamma }_{k}+{\gamma }_{-k}^{{\dagger}}{\gamma }_{-k}-1\right)\\ & \quad +(g-1)({c}_{0}^{{\dagger}}{c}_{0}-{c}_{0}{c}_{0}^{{\dagger}})+(g+1)({c}_{\pi }^{{\dagger}}{c}_{\pi }-{c}_{\pi }{c}_{\pi }^{{\dagger}}),\end{align} \tag{ 18 }$$

where, for all k ∈ k⁻ the equations (14a), (14b) and (14c) take the same forms. Additionally, for 0, π momenta

$$\begin{equation}{{\epsilon}}_{0}=g-1,\hspace{25.0pt}{{\epsilon}}_{\pi }=g+1,\end{equation} \tag{ 19 }$$

$$\begin{equation}{\vartheta }_{0}=0,\hspace{25.0pt}{\vartheta }_{\pi }=0.\end{equation} \tag{ 20 }$$

The ground state in this subspace and its energy read

$$\begin{equation}\vert \gamma ,{g\rangle }^{-}={c}_{0}^{{\dagger}}\prod\limits _{k\in {\mathbf{k}}^{-}}\left[\mathrm{cos}\left(\frac{{\vartheta }_{k}}{2}\right)-\mathrm{sin}\left(\frac{{\vartheta }_{k}}{2}\right){c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}\right]\left\vert \mathrm{v}\mathrm{a}\mathrm{c}\right\rangle .\end{equation} \tag{ 21a }$$

$$\begin{equation}{E}_{0}^{-}\left(\gamma ,g\right)=-2\sum\limits _{k\in {\mathbf{k}}^{-}}{{\epsilon}}_{k}\left(\gamma ,g\right)-2.\end{equation} \tag{ 21b }$$

In the limit case γ = 1, the ground state of the total Hamiltonian, for an even number of spins, always lies in the positive-parity subspace [26]. In the case of the general XY model, this is not always true [31]. Later on, we will analyze in detail the Ising model, in which the energy gap is well defined

$$\begin{equation}{\Delta}(g)={\Delta}(\gamma =1,g)={E}_{0}^{-}(g)-{E}_{0}^{+}(g){\geqslant}0.\end{equation} \tag{ 22 }$$

In the following, the explicit expression for this gap at the critical point (g = 1) will be useful [26]

$$\begin{equation}{\Delta}(g=1)=2\enspace \mathrm{tan}\left(\frac{\pi }{4N}\right)\approx \frac{\pi }{2N}.\end{equation} \tag{ 23 }$$

In the ferromagnetic phase the gap Δ(g) vanishes exponentially with the system size and we refer to it as a 'symmetry-breaking gap' to distinguish it from the 'dynamical gap', which is defined as the lowest excitation within the positive-parity subspace [32]:

$$\begin{equation}\mathrm{g}\mathrm{a}\mathrm{p}(g)=\mathrm{g}\mathrm{a}\mathrm{p}(\gamma =1,g)=\sqrt{{g}^{2}-2g\enspace \mathrm{cos}(\pi /N)+1}.\end{equation} \tag{ 24 }$$

At the critical point, this gap behaves as

$$\begin{equation}\mathrm{g}\mathrm{a}\mathrm{p}(g=1)\approx \frac{4\pi }{N}.\end{equation} \tag{ 25 }$$

Therefore, at the critical point, the dynamical gap is approximately eight times bigger than the symmetry-breaking gap. By contrast, in the ferromagnetic phase (g < 1) the symmetry-breaking gap is negligible in comparison with the dynamical gap. As we shall see, this fact is important in approximating the Gibbs state. The importance of the symmetry-breaking gap and its dependence on the boundary conditions was first highlighted in [33].

Calculations in the positive-parity subspace are sufficient at zero temperature. Moreover, the positive parity is preserved during dynamics. However, at finite-temperature both positive and negative parity subspaces play a role, given their contribution to the unnormalized Gibbs state

$$\begin{equation}{\tilde{\rho}}_{\text{Gibbs}}(\beta ,\gamma ,g)=\mathrm{exp}\left[-\beta H(\gamma ,g)\right].\end{equation} \tag{ 26 }$$

Here and elsewhere, we use a tilde to denote unnormalized density matrices.

In this section, we briefly recall the methods and results from [25]. To begin, we note that from the transformation to fermionic quasiparticles, the total Hilbert space $\mathcal{H}$ can be written as a tensor product of four-dimensional Hilbert subspaces corresponding to each pair of momenta

$$\begin{equation}\mathcal{H}=\underset{k\in {\mathbf{k}}^{+}}{\bigotimes}{\mathcal{H}}_{k},\quad {\mathcal{H}}_{k}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert {0\rangle }_{k},{c}_{k}^{{\dagger}}\vert {0\rangle }_{k},{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k},{c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k}\right\},\end{equation} \tag{ 27 }$$

where |0⟩_k is annihilated by c_k and c_−k . Note that dimensions match: there are N/2 momenta, which gives the total 2^N dimension. For negative subspace one has 0 and π momenta with corresponding two-dimensional Hilbert spaces span by $\left\{\vert {0\rangle }_{0},{c}_{0}^{{\dagger}}\vert {0\rangle }_{0}\right\}$ and $\left\{\vert {0\rangle }_{\pi },{c}_{\pi }^{{\dagger}}\vert {0\rangle }_{\pi }\right\}$, together with (N/2) − 1 momenta in k⁻. On the other hand, the subspaces of the given parity have dimension two times smaller. Therefore, it is clear that usual tensor product is not adapted for manipulations involving a definite parity. In order to handle this, we introduced the operations of 'positive' and 'negative' tensor products which pick only vectors of correct parity. They can be defined recursively:

$$\begin{equation}\mathcal{P}\left({\mathcal{H}}_{{k}_{1}}\right)={\mathcal{H}}_{{k}_{1}}^{(p)},\hspace{25.0pt}\mathcal{N}\left({\mathcal{H}}_{{k}_{1}}\right)={\mathcal{H}}_{{k}_{1}}^{(n)},\end{equation} \tag{ 28 }$$

$$\begin{align}\mathcal{P}\left(\underset{i=1}{\overset{n+1}{\bigotimes}}{\mathcal{H}}_{{k}_{i}}\right)=\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{\mathcal{H}}_{{k}_{i}}\right)\bigotimes{\mathcal{H}}_{{k}_{n+1}}^{(p)}\oplus \mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{\mathcal{H}}_{{k}_{i}}\right)\bigotimes{\mathcal{H}}_{{k}_{n+1}}^{(n)},\quad n{\geqslant}1\\ \mathcal{N}\left(\underset{i=1}{\overset{n+1}{\bigotimes}}{\mathcal{H}}_{{k}_{i}}\right)=\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{\mathcal{H}}_{{k}_{i}}\right)\bigotimes{\mathcal{H}}_{{k}_{n+1}}^{(n)}\oplus \mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{\mathcal{H}}_{{k}_{i}}\right)\bigotimes{\mathcal{H}}_{{k}_{n+1}}^{(p)},\quad n{\geqslant}1,\end{align}$$

where ${\mathcal{H}}_{k}^{(p)}$ and ${\mathcal{H}}_{k}^{(n)}$ are subspaces of ${\mathcal{H}}_{k}$ spanned by all positive parity an negative parity vectors:

$$\begin{align*}{\mathcal{H}}_{k}^{(p)}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert {0\rangle }_{k},{c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k}\right\},\hspace{25.0pt}{\mathcal{H}}_{k}^{(n)}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{{c}_{k}^{{\dagger}}\vert {0\rangle }_{k},{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k}\right\},\\ {\mathcal{H}}_{0}^{(p)}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert {0\rangle }_{0}\right\},\hspace{25.0pt}{\mathcal{H}}_{0}^{(n)}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{{c}_{0}^{{\dagger}}\vert {0\rangle }_{0}\right\},\\ {\mathcal{H}}_{\pi }^{(p)}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert {0\rangle }_{\pi }\right\},\hspace{25.0pt}{\mathcal{H}}_{\pi }^{(n)}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{{c}_{\pi }^{{\dagger}}\vert {0\rangle }_{\pi }\right\}.\end{align*}$$

Similarly, one can define analogous relations for positive and negative parity part of a tensor product of operators. Let us assume that each O_k has zero matrix elements between vectors of different parity (such as ⟨0_k |c_k O_k |0⟩_k or $\langle {0}_{k}\vert {c}_{k}{O}_{k}{c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k}$ for instance). For such operators we define a positive and negative part in an intuitive way through the following relations:

$$\begin{align*}\left.{O}_{k}^{(p)}\right\vert {\mathcal{H}}_{k}^{(p)}=\left.{O}_{k}\right\vert {\mathcal{H}}_{k}^{(p)},\hspace{25.0pt}\left.{O}_{k}^{(p)}\right\vert {\mathcal{H}}_{k}^{(n)}=0,\\ \left.{O}_{k}^{(n)}\right\vert {\mathcal{H}}_{k}^{(n)}=\left.{O}_{k}\right\vert {\mathcal{H}}_{k}^{(n)},\hspace{25.0pt}\left.{O}_{k}^{(n)}\right\vert {\mathcal{H}}_{k}^{(p)}=0,\end{align*}$$

with recursive relations taking the form

$$\begin{equation}\mathcal{P}\left({O}_{{k}_{1}}\right)={O}_{{k}_{1}}^{(p)},\hspace{25.0pt}\mathcal{N}\left({O}_{{k}_{1}}\right)={\mathcal{O}}_{{k}_{1}}^{(n)},\end{equation} \tag{ 29 }$$

$$\begin{equation*}\mathcal{P}\left(\underset{i=1}{\overset{n+1}{\bigotimes}}{O}_{{k}_{i}}\right)=\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)\bigotimes{O}_{{k}_{n+1}}^{(p)}+\mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)\bigotimes{O}_{{k}_{n+1}}^{(n)},\quad n{\geqslant}1,\end{equation*}$$

$$\begin{equation*}\mathcal{N}\left(\underset{i=1}{\overset{n+1}{\bigotimes}}{\mathcal{O}}_{{k}_{i}}\right)=\mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)\bigotimes{O}_{{k}_{n+1}}^{(p)}+\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)\bigotimes{O}_{{k}_{n+1}}^{(n)},\quad n{\geqslant}1.\end{equation*}$$

These relations are true provided that momenta k₁, ..., k_n are all relevant for one subspace (positive or negative parity). At this point we emphasize that the condition of vanishing 'mixing' matrix elements is an important restriction that does not apply to many important properties, such as the longitudinal magnetizations ${\sigma }_{i}^{x}$ or ${\sigma }_{i}^{y}$. Dealing with such operators is particularly difficult because subspaces with different momenta mix. For some methods and results involving ground states, see for example [34]. We shall make use of the following propositions:

Proposition 1. For every operators O_k and R_k :

$$\begin{equation}\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)\enspace \mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{R}_{{k}_{i}}\right)=\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\enspace {R}_{{k}_{i}}\right),\end{equation} \tag{ 30 }$$

$$\begin{equation}\mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)\enspace \mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{R}_{{k}_{i}}\right)=\mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\enspace {R}_{{k}_{i}}\right).\end{equation} \tag{ 31 }$$

Proposition 2. The exponentials of restricted tensor products obey:

$$\begin{equation}\mathrm{exp}\enspace \mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)=\mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}\mathrm{exp}\enspace {O}_{{k}_{i}}\right),\end{equation} \tag{ 32 }$$

$$\begin{equation}\mathrm{exp}\enspace \mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)=\mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}\mathrm{exp}\enspace {O}_{{k}_{i}}\right).\end{equation} \tag{ 33 }$$

Proposition 3. Traces of restricted tensor products read:

$$\begin{equation}\mathrm{Tr}\enspace \mathcal{P}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)=\frac{1}{2}\left(\prod\limits _{i=1}^{n}\mathrm{Tr}\enspace {O}_{{k}_{i}}+\prod\limits _{i=1}^{n}(\mathrm{Tr}\enspace {O}_{{k}_{i}}^{(p)}-\mathrm{Tr}\enspace {O}_{{k}_{i}}^{(n)})\right),\end{equation} \tag{ 34 }$$

$$\begin{equation}\mathrm{Tr}\enspace \mathcal{N}\left(\underset{i=1}{\overset{n}{\bigotimes}}{O}_{{k}_{i}}\right)=\frac{1}{2}\left(\prod\limits _{i=1}^{n}\mathrm{Tr}\enspace {O}_{{k}_{i}}-\prod\limits _{i=1}^{n}(\mathrm{Tr}\enspace {O}_{{k}_{i}}^{(p)}-\mathrm{Tr}\enspace {O}_{{k}_{i}}^{(n)})\right).\end{equation} \tag{ 35 }$$

These propositions suffice to derive the formula for fidelity between arbitrary thermal states. We finish this section with an example:

Example 1. Canonical Gibbs state

The complete Hamiltonian can be written in the form:

$$\begin{equation}H(\gamma ,g)={H}^{+}(\gamma ,g)\oplus {H}^{-}(\gamma ,g),\end{equation} \tag{ 36 }$$

where

$$\begin{equation}{H}^{+}(\gamma ,g)=\mathcal{P}\left(\underset{k\in {\mathbf{k}}^{+}}{\bigotimes}{H}_{k}(\gamma ,g)\right),\hspace{25.0pt}{H}^{-}(\gamma ,g)=\mathcal{N}\left(\underset{k\in {\mathbf{k}}^{-}\cup \left\{0,\pi \right\}}{\bigotimes}{H}_{k}(\gamma ,g)\right).\end{equation} \tag{ 37 }$$

Operators H_k have the following matrix representations

$$\begin{equation}{H}_{k}(\gamma ,g)=2\left(\begin{matrix}\hfill \mathrm{cos}\enspace k-g\hfill & \hfill \gamma \enspace \mathrm{sin}\enspace k\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \gamma \enspace \mathrm{sin}\enspace k\hfill & \hfill g-\mathrm{cos}\enspace k\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 38 }$$

$$\begin{equation}{H}_{0}(g)=\left(\begin{matrix}\hfill 1-g\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill g-1\hfill \end{matrix}\right),\end{equation} \tag{ 39 }$$

$$\begin{equation}{H}_{\pi }(g)=\left(\begin{matrix}\hfill -1-g\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1+g\hfill \end{matrix}\right),\end{equation} \tag{ 40 }$$

in the basis $\left\{\vert {0\rangle }_{k},{c}_{k}^{{\dagger}}{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k},{c}_{k}^{{\dagger}}\vert {0\rangle }_{k},{c}_{-k}^{{\dagger}}\vert {0\rangle }_{k}\right\}$ and $\left\{\vert {0\rangle }_{0,\pi },{c}_{0,\pi }^{{\dagger}}\vert {0\rangle }_{0,\pi }\right\}$, respectively. Using proposition 2, we can write the unnormalized thermal Gibbs state as:

$$\begin{equation}{\tilde{\rho}}_{\text{Gibbs}}(\beta ,\gamma ,g)=\mathrm{exp}\left(-\beta H(\gamma ,g)\right)=\mathcal{P}\left(\underset{k\in {\mathbf{k}}^{+}}{\bigotimes}{{\tilde{\rho }}_{k}}\right)\oplus \mathcal{N}\left(\underset{k\in {\mathbf{k}}^{-}\cup \left\{0,\pi \right\}}{\bigotimes}{{\tilde{\rho }}_{k}}\right),\end{equation} \tag{ 41 }$$

where

$$\begin{equation}{{\tilde{\rho }}_{k}}=\mathrm{exp}\left[-2\beta \left(\begin{matrix}\hfill \mathrm{cos}\enspace k-g\hfill & \hfill \gamma \enspace \mathrm{sin}\enspace k\hfill \\ \hfill \gamma \enspace \mathrm{sin}\enspace k\hfill & \hfill g-\mathrm{cos}\enspace k\hfill \end{matrix}\right)\right]\oplus {\mathbb{I}}_{2},\quad \end{equation} \tag{ 42 }$$

$$\begin{equation}{\tilde{\rho}}_{0}=\left(\begin{matrix}\hfill \mathrm{exp}\left[\beta (g-1)\right]\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathrm{exp}\left[\beta (1-g)\right]\hfill \end{matrix}\right),\end{equation} \tag{ 43 }$$

$$\begin{equation}{\tilde{\rho}}_{\pi }=\left(\begin{matrix}\hfill \mathrm{exp}\left[-\beta (g+1)\right]\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathrm{exp}\left[\beta (g+1)\right]\hfill \end{matrix}\right).\end{equation} \tag{ 44 }$$

The even and odd parity parts of ρ_k are

$$\begin{equation}{\tilde{\rho}}_{k}^{(p)}=\mathrm{exp}\left[-2\beta \left(\begin{matrix}\hfill \mathrm{cos}\enspace k-g\hfill & \hfill \gamma \enspace \mathrm{sin}\enspace k\hfill \\ \hfill \gamma \enspace \mathrm{sin}\enspace k\hfill & \hfill g-\mathrm{cos}\enspace k\hfill \end{matrix}\right)\right],\end{equation} \tag{ 45 }$$

$$\begin{equation}{\tilde{\rho}}_{k}^{(n)}={\mathbb{I}}_{2},\hspace{25.0pt}{\tilde{\rho}}_{0}^{(p)}=\mathrm{exp}\left[\beta (g-1)\right],\hspace{25.0pt}{\tilde{\rho}}_{0}^{(n)}=\mathrm{exp}\left[\beta (1-g)\right].\end{equation} \tag{ 46 }$$

Using equation (14a) one has $\mathrm{Tr}\left({\tilde{\rho}}_{k}\right)=2\enspace \mathrm{cosh}\left(2\beta {{\epsilon}}_{k}\right)$. Thanks to proposition 3, we can easily write down the formula for the partition function, that is the trace of the unnormalized Gibbs state:

$$\begin{align}Z\left(\beta ,\gamma ,g\right)& =\mathrm{Tr}\left({\tilde{\rho}}_{\text{Gibbs}}(\beta ,\gamma ,g)\right)\\ & =\frac{1}{2}\left[\prod\limits _{k\in {\mathbf{K}}^{+}}2\enspace \mathrm{cosh}\left(\beta {{\epsilon}}_{k}\left(\gamma ,g\right)\right)+\prod\limits _{k\in {\mathbf{K}}^{+}}2\enspace \mathrm{sinh}\left(\beta {{\epsilon}}_{k}\left(\gamma ,g\right)\right)\right.\\ & \quad +\left.\prod\limits _{k\in {\mathbf{K}}^{-}}2\enspace \mathrm{cosh}\left(\beta {{\epsilon}}_{k}\left(\gamma ,g\right)\right)-\prod\limits _{k\in {\mathbf{K}}^{-}}2\enspace \mathrm{sinh}\left(\beta {{\epsilon}}_{k}\left(\gamma ,g\right)\right)\right].\end{align} \tag{ 47 }$$

Here, the use of the full sets of momenta K⁺ and K⁻ leads to a compact expression, but we emphasize the need to use the correct formulas for excitations in the modes with 0 and π momenta (19).

Equation (47) can be divided into four parts [24] [positive fermionic $\left({Z}_{F}^{+}(\beta ,\gamma ,g)\right)$, positive boundary $\left({Z}_{B}^{+}(\beta ,\gamma ,g)\right)$, negative fermionic $\left({Z}_{F}^{-}(\beta ,\gamma ,g)\right)$, and negative boundary $\left.\left({Z}_{F}^{+}(\beta ,\gamma ,g)\right)\right]$:

$$\begin{equation}Z(\beta ,\gamma ,g)=\frac{1}{2}\left[{Z}_{F}^{+}(\beta ,\gamma ,g)+{Z}_{B}^{+}(\beta ,\gamma ,g)+{Z}_{F}^{-}(\beta ,\gamma ,g)-{Z}_{B}^{-}(\beta ,\gamma ,g)\right].\end{equation} \tag{ 48 }$$

These four terms originate from the trace formulas in proposition 3. It is thus reasonable to assume that many physically relevant quantities (such as characteristic functions of observables or fidelity susceptibility) can be decomposed in an analogous form. It has been argued that in the thermodynamic limit only ${Z}_{F}^{+}$ is relevant, i.e. that the partition function is governed by the positive fermionic part corresponding to the first part of equation (34) [23]. While this is true in many cases, one has to be careful; see [25] for details. Explicitly, the positive fermionic part of the partition function has the form

$$\begin{equation}{Z}_{F}(\beta ,\gamma ,g)=\prod\limits _{k\in {\mathbf{K}}^{+}}2\enspace \mathrm{cosh}\left(\beta {{\epsilon}}_{k}\left(\gamma ,g\right)\right).\end{equation} \tag{ 49 }$$

We note that decompositions resembling (48) appear naturally in equations describing fidelities. In what follows, we use the acronym PPA to refer to the positive parity approximation that includes only positive fermionic contribution. As we shall see, in the study of the fidelity susceptibility the exact expression involving all contributions is significantly different from the PPA approximation frequently used in the literature [2, 3, 7, 18, 19, 35–37].

The algebraic method we have presented for the exact treatment of the partition function is natural and elementary, but not the only one. Alternative derivations are possible making use of Grassmann variables, see [24] or group theory [38].

In this section, we derive the exact expression for the Uhlmann fidelity between thermal states in the XY model with periodic boundary conditions. The method can be generalized to states that are exponentials or quadratic expressions in fermionic operators.

We begin by restating the definition of the Uhlmann fidelity between two density matrices ρ and σ [39]

$$\begin{equation}F\left(\rho ,\sigma \right)=\mathrm{Tr}\left(\sqrt{\sqrt{\rho }\sigma \sqrt{\rho }}\right).\end{equation} \tag{ 50 }$$

As ${\sqrt{\rho }}_{k}$ is the square root of ρ_k , it follows that

$$\begin{equation}{\left[\mathcal{P}\left(\underset{k}{\bigotimes}\sqrt{{\tilde{\rho}}_{k}}\right)\right]}^{2}=\mathcal{P}\left(\underset{k}{\bigotimes}\sqrt{{\tilde{\rho}}_{k}}\sqrt{{\tilde{\rho}}_{k}}\right)=\mathcal{P}\left(\underset{k}{\bigotimes}{\tilde{\rho}}_{k}\right).\end{equation} \tag{ 51 }$$

As a result, one finds

$$\begin{equation}\sqrt{\mathcal{P}\left(\underset{k}{\bigotimes}{\tilde{\rho}}_{k}\right)}=\mathcal{P}\left(\underset{k}{\bigotimes}\sqrt{{\tilde{\rho}}_{k}}\right),\end{equation} \tag{ 52 }$$

and similarly for $\sqrt{\mathcal{N}\left({\bigotimes}_{k}{\tilde{\rho}}_{k}\right)}$. Using the results from example 1 and the unnormalized thermal states ${\tilde{\sigma }}_{\text{Gibbs}}$, ${\tilde{\rho}}_{\text{Gibbs}}$,

$$\begin{equation}{\tilde{\rho}}_{\text{Gibbs}}={\tilde{\rho}}_{\text{Gibbs}}^{+}\oplus {\tilde{\rho}}_{\text{Gibbs}}^{-},\hspace{25.0pt}\sqrt{{\tilde{\rho}}_{\text{Gibbs}}}=\sqrt{{\tilde{\rho}}_{\text{Gibbs}}^{+}}\oplus \sqrt{{\tilde{\rho}}_{\text{Gibbs}}^{-}},\end{equation} \tag{ 53 }$$

which yields a decomposition of the exact fidelity into the sum of the components with positive and negative parity

$$\begin{equation}F\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)={F}^{+}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)+{F}^{-}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right).\end{equation} \tag{ 54 }$$

Using equation (52) and proposition 1 again, we find

$$\begin{equation}{F}^{+}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)=\mathrm{Tr}\left(\mathcal{P}\left(\underset{k\in {\mathbf{k}}^{+}}{\bigotimes}\sqrt{\sqrt{{\tilde{\rho}}_{k}}{\tilde{\sigma }}_{k}\sqrt{{\tilde{\rho}}_{k}}}\right)\right)\end{equation} \tag{ 55 }$$

$$\begin{equation}{F}^{-}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)=\mathrm{Tr}\left(\mathcal{N}\left(\underset{k\in {\mathbf{k}}^{-}\cup \left\{0,\pi \right\}}{\bigotimes}\sqrt{\sqrt{{\tilde{\rho}}_{k}}{\tilde{\sigma }}_{k}\sqrt{{\tilde{\rho}}_{k}}}\right)\right).\end{equation} \tag{ 56 }$$

Next we make use of proposition 3 to write down explicit expressions for general thermal states, with different temperatures and field values. We consider the fidelity between two thermal states characterized by field values g^ρ , g^σ , anisotropy coefficients γ^ρ , γ^σ and inverse temperature β^ρ , β^σ , that is:

$$\begin{equation}{\tilde{\rho}}_{k}=\mathrm{exp}\left[-2{\beta }^{\rho }\left(\begin{matrix}\hfill \mathrm{cos}\enspace k-{g}^{\rho }\hfill & \hfill {\gamma }^{\rho }\enspace \mathrm{sin}\enspace k\hfill \\ \hfill {\gamma }^{\rho }\enspace \mathrm{sin}\enspace k\hfill & \hfill {g}^{\rho }-\mathrm{cos}\enspace k\hfill \end{matrix}\right)\right]\oplus {\mathbb{I}}_{2}={\tilde{\rho}}_{k}^{(p)}\oplus {\mathbb{I}}_{2},\end{equation} \tag{ 57 }$$

$$\begin{equation}{\tilde{\sigma }}_{k}=\mathrm{exp}\left[-2{\beta }^{\sigma }\left(\begin{matrix}\hfill \mathrm{cos}\enspace k-{g}^{\sigma }\hfill & \hfill {\gamma }^{\sigma }\enspace \mathrm{sin}\enspace k\hfill \\ \hfill {\gamma }^{\sigma }\enspace \mathrm{sin}\enspace k\hfill & \hfill {g}^{\sigma }-\mathrm{cos}\enspace k\hfill \end{matrix}\right)\right]\oplus {\mathbb{I}}_{2}={\tilde{\sigma }}_{k}^{(p)}\oplus {\mathbb{I}}_{2}.\end{equation} \tag{ 58 }$$

Calculation of $\mathrm{Tr}\left(\sqrt{\sqrt{{\rho }_{k}}{\sigma }_{k}\sqrt{{\rho }_{k}}}\right)$ for thermal states boils down to calculation of fidelity between qubit states, which has a particularly simple form [20]:

$$\begin{equation}\mathrm{Tr}\left(\sqrt{\sqrt{{\tilde{\rho}}_{k}^{(p)}}{\tilde{\sigma }}_{k}^{(p)}\sqrt{{\tilde{\rho}}_{k}^{(p)}}}\right)={\left[\mathrm{Tr}\left({\tilde{\rho}}_{k}^{(p)}{\tilde{\sigma }}_{k}^{\left.(p)\right)}\right)+2\enspace \mathrm{det}\left({\tilde{\rho}}_{k}^{(p)}{\tilde{\sigma }}_{k}^{(p)}\right)\right]}^{1/2}.\end{equation} \tag{ 59 }$$

Using the Bogoliubov energy and angles for the states ρ_k and σ_k , we can derive compact expressions for fidelity of thermal states. Although the overall structure of the formulas is simple, there are many parameters and expressions become lengthy. For clarity we introduce the following notation:

$$\begin{equation}{u}_{k}\left(\rho \vert \sigma \right)=2\enspace \mathrm{cosh}\left(\frac{{\beta }^{\rho }{{\epsilon}}_{k}^{\rho }+{\beta }^{\sigma }{{\epsilon}}_{k}^{\sigma }}{2}\right)\mathrm{cos}\left(\frac{{\vartheta }_{k}^{\rho }-{\vartheta }_{k}^{\sigma }}{2}\right),\end{equation} \tag{ 60 }$$

$$\begin{equation}{v}_{k}\left(\rho \vert \sigma \right)=2\enspace \mathrm{cosh}\left(\frac{{\beta }^{\rho }{{\epsilon}}_{k}^{\rho }-{\beta }^{\sigma }{{\epsilon}}_{k}^{\sigma }}{2}\right)\mathrm{sin}\left(\frac{{\vartheta }_{k}^{\rho }-{\vartheta }_{k}^{\sigma }}{2}\right).\end{equation} \tag{ 61 }$$

The Uhlmann fidelity between thermal states ρ_Gibbs and σ_Gibbs, respectively characterized by the parameters β^ρ , g^ρ , γ^ρ and β^σ , g^σ , γ^σ , reads

$$\begin{equation}F\left({\rho }_{\text{Gibbs}},{\sigma }_{\text{Gibbs}}\right)=\frac{{F}^{+}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)+{F}^{-}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)}{\sqrt{Z({\beta }^{\rho },{g}^{\rho },{\gamma }^{\rho })Z({\beta }^{\sigma },{g}^{\sigma },{\gamma }^{\sigma })}}.\end{equation} \tag{ 62 }$$

Here, the positive part equals

$$\begin{align}{F}^{+}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)& =\frac{1}{2}\left[\prod\limits _{k\in {\mathbf{k}}^{+}}\left(\sqrt{{u}_{k}^{2}\left(\rho \vert \sigma \right)+{v}_{k}^{2}\left(\rho \vert \sigma \right)}+2\right)\right.\\ & \quad +\left.\prod\limits _{k\in {\mathbf{k}}^{+}}\left(\sqrt{{u}_{k}^{2}\left(\rho \vert \sigma \right)+{v}_{k}^{2}\left(\rho \vert \sigma \right)}-2\right)\right],\end{align} \tag{ 63 }$$

and the negative part is given by

$$\begin{align}{F}^{-}\left({\tilde{\rho}}_{\text{Gibbs}},{\tilde{\sigma }}_{\text{Gibbs}}\right)& =\frac{1}{2}\left[4\enspace \mathrm{cosh}\left(\frac{{\beta }^{\rho }{{\epsilon}}_{0}^{\rho }+{\beta }^{\sigma }{{\epsilon}}_{0}^{\sigma }}{2}\right)\mathrm{cosh}\left(\frac{{\beta }^{\rho }{{\epsilon}}_{\pi }^{\rho }+{\beta }^{\sigma }{{\epsilon}}_{\pi }^{\sigma }}{2}\right)\right.\\ & \quad {\times}\prod\limits _{k\in {\mathbf{k}}^{-}}\left(\sqrt{{u}_{k}^{2}\left(\rho \vert \sigma \right)+{v}_{k}^{2}\left(\rho \vert \sigma \right)}+2\right)\\ & \quad -4\enspace \mathrm{sinh}\left(\frac{{\beta }^{\rho }{{\epsilon}}_{0}^{\rho }+{\beta }^{\sigma }{{\epsilon}}_{0}^{\sigma }}{2}\right)\mathrm{sinh}\left(\frac{{\beta }^{\rho }{{\epsilon}}_{\pi }^{\rho }+{\beta }^{\sigma }{{\epsilon}}_{\pi }^{\sigma }}{2}\right)\\ & \quad {\times}\left.\prod\limits _{k\in {\mathbf{k}}^{-}}\left(\sqrt{{u}_{k}^{2}\left(\rho \vert \sigma \right)+{v}_{k}^{2}\left(\rho \vert \sigma \right)}-2\right)\right],\end{align} \tag{ 64 }$$

where excitations _k , _0,π and Bogoliubov angles ϑ_k can be calculated with standard formulas (14a), (14b) and (14c) and the exact partition function Z is given by equation (47).

As in the case of partition function, we can single out the positive fermionic part

$$\begin{equation}{F}_{F}^{+}\left({\rho }_{\text{Gibbs}},{\sigma }_{\text{Gibbs}}\right)=\frac{{\prod }_{k\in {\mathbf{k}}^{+}}\left(\sqrt{{u}_{k}^{2}\left(\rho \vert \sigma \right)+{v}_{k}^{2}\left(\rho \vert \sigma \right)}+2\right)}{\sqrt{{Z}_{F}^{+}({\beta }^{\rho },{\gamma }^{\rho },{g}^{\rho }){Z}_{F}^{+}({\beta }^{\sigma },{\gamma }^{\sigma },{g}^{\sigma })}}\quad (\mathrm{P}\mathrm{P}\mathrm{A}),\end{equation} \tag{ 65 }$$

where ${Z}_{F}^{+}$ is given by the first term in decomposition (48) and explicitly reads

$$\begin{equation}{Z}_{F}^{+}\left(\beta ,\gamma ,g\right)=\prod\limits _{k\in {\mathbf{K}}^{+}}2\enspace \mathrm{cosh}(\beta {{\epsilon}}_{k}\left(\gamma ,g\right))\quad (\mathrm{P}\mathrm{P}\mathrm{A}).\end{equation} \tag{ 66 }$$

The PPA formula (65) has been commonly used in literature [18, 19]. In the following sections, we are going to show the limits of this approximation by comparing it to the exact formula for the fidelity. As it turns out, the latter includes important physics in the intermediate temperature regime, which is missed by the simplified expression (65). We will use the notation

$$\begin{equation}F\left({\rho }_{\text{Gibbs}},{\sigma }_{\text{Gibbs}}\right)=F\left({\beta }^{\rho },{\gamma }^{\rho },{g}^{\rho }\vert {\beta }^{\sigma },{\gamma }^{\sigma },{g}^{\sigma }\right).\end{equation} \tag{ 67 }$$

There are two limiting cases of these general expressions that are particularly relevant. One concerns the Uhlmann fidelity between two thermal states at equal inverse temperature but a different value of the external control parameter g. This limit is natural to quantify the role of thermal excitations above the ground state. The second case concerns the Uhlmann fidelity between two thermal states with common control parameter g but different inverse temperatures. The explicit form of the Uhlmann fidelity can be easily derived in this case given that the two density matrices commute, i.e. [ρ_Gibbs, σ_Gibbs] = 0. As a result, the Uhlmann fidelity can be written down in terms of the partition function

$$\begin{equation}F\left(\beta ,\gamma ,g\vert {\beta }^{\prime },\gamma ,g\right)=\frac{Z\left(\frac{\beta +{\beta }^{\prime }}{2},\gamma ,g\right)}{\sqrt{Z(\beta ,\gamma ,g)Z({\beta }^{\prime },\gamma ,g)}}.\end{equation} \tag{ 68 }$$

By using the full expression for the partition function (48), one obtains the exact result, without the need to resort to the PPA based on the simplified partition function in the even parity subspace, previously considered, e.g. in [2, 3, 5, 9, 18, 19, 35].

We independently verified the correctness of the formula (62) by showing that it reproduces the results obtained numerically by exact diagonalization for small system (N = 6, 8, 10). For details of the numerical simulations, see appendix A.

In this section we analyze the fidelity susceptibility with respect to external field g for fixed inverse temperature β. We also set γ = 1 (case of Ising model in transverse field). We thus focus on the quantity

$$\begin{equation}\chi (\beta ,g)=-\left.\frac{{\mathrm{d}}^{2}}{\mathrm{d}{\delta }^{2}}F\left(\beta ,\gamma =1,g\vert \beta ,\gamma =1,g+\delta \right)\right\vert \delta =0,\end{equation} \tag{ 69 }$$

which provides the leading nontrivial term in the expansion

$$\begin{equation}F\left(\beta ,\gamma =1,g\vert \beta ,\gamma =1,g+\delta \right)=1-\frac{1}{2}\chi (\beta ,g){\delta }^{2}+\dots .\end{equation} \tag{ 70 }$$

By contrast, the fidelity susceptibility obtained from the PPA reads

$$\begin{equation}{\chi }_{F}^{+}(\beta ,g)=-\left.\frac{{\mathrm{d}}^{2}}{\mathrm{d}{\delta }^{2}}{F}_{F}^{+}\left(\beta ,\gamma =1,g\vert \beta ,\gamma =1,g+\delta \right)\right\vert \delta =0\quad (\mathrm{P}\mathrm{P}\mathrm{A}).\end{equation} \tag{ 71 }$$

A representative plot of the Uhlmann fidelity is shown in figure 1 for two Gibbs states that differ in the value of magnetic field or inverse temperature. Significant discrepancies between the exact expression and the PPA are manifested at low temperatures. The choice of parameters has not been optimized to maximize the latter; discrepancies can be bigger, especially for general XY model (γ ≠ 1), where states with negative parity can appear more often then for Ising model [31]. In the following, we examine the nature of the discrepancies by a detailed analysis of fidelity susceptibility.

Zoom In Zoom Out Reset image size

In figure 2, we compare the exact and PPA results and show the dependence of fidelity susceptibility on the magnetic field for different β. Because the fidelity susceptibility is an even function of g, χ(β, g) = χ(β, − g), we present the results only for g > 0. The PPA yields qualitatively different results from the exact expression, with the difference being more pronounced in the low-temperature limit. Not only the value of the maximum differs, the PPA also leads to a shift of its location. It is clear that the discrepancy is enhanced in the vicinity of g_c = 1. Figure 3 shows the dependence of the fidelity susceptibility on β at the critical point g_c = 1. It is found that the enhancement of the exact fidelity susceptibility occurs for a specific value of β that depends on the system size, a feature we analyze next.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1. Two-level approximation of the Gibbs state

Deviations between the PPA and the exact results are particularly pronounced at intermediate and low temperatures. We next introduce an accurate approximation in this regime by truncating of the Gibbs state, taking only the two lowest energy states into account. We refer to it as the two-level approximation (TLA) which yields the following truncation for the Gibbs state:

$$\begin{equation}{\rho }_{\text{TLA}}(\beta ,g)=\frac{1}{{Z}_{\text{TLA}}(\beta ,g)}\left({\mathrm{e}}^{-\beta {E}_{0}^{+}(g)}\vert {g}^{+}\rangle \langle {g}^{+}\vert +{\mathrm{e}}^{-\beta {E}_{0}^{-}(g)}\vert {g}^{-}\rangle \langle {g}^{-}\vert \right),\end{equation} \tag{ 72 }$$

where

$$\begin{equation}{Z}_{\text{TLA}}(\beta ,g)={\mathrm{e}}^{-\beta {E}_{0}^{+}(g)}+{\mathrm{e}}^{-\beta {E}_{0}^{-}(g)}.\end{equation} \tag{ 73 }$$

Analogously, the neighboring thermal state can be written as

$$\begin{align}{\rho }_{\text{TLA}}(\beta ,g+\delta )& =\frac{1}{{Z}_{\text{TLA}}(\beta ,g+\delta )}\left({\mathrm{e}}^{-\beta {E}_{0}^{+}(g+\delta )}\vert g+{\delta }^{+}\rangle \langle g+{\delta }^{+}\vert \right.\\ & \quad +\left.{\mathrm{e}}^{-\beta {E}_{0}^{-}(g+\delta )}\vert g+{\delta }^{-}\rangle \langle g+{\delta }^{-}\vert \right).\end{align} \tag{ 74 }$$

In principle, the TLA should work well for big β, i.e. in the low-temperature regime. The accuracy of this approximation relies on the structure of energy levels and the relation between the dynamical and symmetry-breaking gap. Recall that the symmetry-breaking gap vanishes exponentially with N for g < 1 and is eight times smaller than the dynamical gap at the critical point. Because the energy ${E}_{1}^{+}$ of the first excited state in the positive-parity subspace is well separated from ${E}_{0}^{+}$ and ${E}_{0}^{-}$, one expects the contribution from ${\text{e}}^{-\beta {E}_{1}^{+}}$ and higher energy states to be small in comparison to equation (72).

The Uhlmann fidelity between thermal states in the TLA reads

$$\begin{align}{F}_{\text{TLA}}\left(\beta ,g\vert \beta ,g+\delta \right)& =\frac{1}{\sqrt{1+{\mathrm{e}}^{-\beta {\Delta}(g)}}}\frac{{F}_{0}^{+}\left(g\vert g+\delta \right)}{\sqrt{1+{\mathrm{e}}^{-\beta {\Delta}(g+\delta )}}}\\ & \quad +\frac{1}{\sqrt{1+{\mathrm{e}}^{\beta {\Delta}(g)}}}\frac{{F}_{0}^{-}\left(g\vert g+\delta \right)}{\sqrt{1+{\mathrm{e}}^{\beta {\Delta}(g+\delta )}}}\quad (\mathrm{T}\mathrm{L}\mathrm{A}),\end{align} \tag{ 75 }$$

where Δ(g) is a symmetry breaking gap defined in (22) and ${F}_{0}^{+}$ and ${F}_{0}^{-}$ denote the ground state fidelities in the even and odd parity subspace, respectively:

$$\begin{equation}{F}_{0}^{{\pm}}\left(g\vert g+\delta \right)=\vert \langle {g}^{{\pm}}\vert g+{\delta }^{{\pm}}\rangle \vert =1-\frac{{\delta }^{2}}{2}{\chi }_{0}^{{\pm}}(g)+\dots .\end{equation} \tag{ 76 }$$

The expression (75) can be used as a starting point for computation of approximated fidelity susceptibility. Expressions for ${\chi }_{0}^{{\pm}}$ were derived analytically in [28]:

$$\begin{equation}{\chi }_{0}^{+}(g)=\frac{{N}^{2}}{16{g}^{2}}\frac{{g}^{N}}{{({g}^{N}+1)}^{2}}+\frac{N}{16{g}^{2}}\frac{{g}^{N}-{g}^{2}}{({g}^{N}+1)({g}^{2}-1)},\end{equation} \tag{ 77 }$$

$$\begin{equation}{\chi }_{0}^{-}(g)=-\frac{{N}^{2}}{16{g}^{2}}\frac{{g}^{N}}{{({g}^{N}-1)}^{2}}+\frac{N}{16{g}^{2}}\frac{{g}^{N}+{g}^{2}}{({g}^{N}-1)({g}^{2}-1)}.\end{equation} \tag{ 78 }$$

To explain the temperature dependence of the fidelity susceptibility, we focus on the case g = 1. Although for a finite system the maximum of the fidelity susceptibility is not exactly at g_c = 1, it is very close to 1 and in practice it is convenient to consider χ(β) = χ(β, g = 1). The ground-state positive and negative susceptibilities are given by

$$\begin{equation}{\chi }_{0}^{+}=\frac{1}{32}N(N-1),\hspace{25.0pt}{\chi }_{0}^{-}=\frac{1}{96}({N}^{2}-3N+2).\end{equation} \tag{ 79 }$$

At g = 1, we can use formulas (75) and (76) to find χ(β) in the TLA

$$\begin{equation}{\chi }_{\text{TLA}}(\beta )=\frac{{\chi }_{0}^{+}}{1+{\text{e}}^{-\frac{\beta \pi }{2N}}}+\frac{{\chi }_{0}^{-}}{1+{\text{e}}^{\frac{\beta \pi }{2N}}}+R(\beta ),\end{equation} \tag{ 80 }$$

where

$$\begin{equation}R(\beta )=-\left.\frac{{\mathrm{d}}^{2}}{\mathrm{d}{\delta }^{2}}\left(\frac{1}{\sqrt{1+{\text{e}}^{-\frac{\beta \pi }{2N}}}\sqrt{1+{\mathit{\text{e}}}^{-\beta {\Delta}(1+\delta )}}}+\frac{1}{\sqrt{1+{\text{e}}^{\frac{\beta \pi }{2N}}}\sqrt{1+{\text{e}}^{\beta {\Delta}(1+\delta )}}}\right)\right\vert \delta =0.\end{equation} \tag{ 81 }$$

Here, we used the asymptotic behavior of the symmetry breaking gap at g_c = 1; see equations (22) and (23). Explicit evaluation yields

$$\begin{align}R(\beta )& =\frac{{\beta }^{2}{\text{e}}^{\frac{\beta \pi }{2N}}{{\Delta}}_{\text{c}}^{\prime 2}}{2{\left({\text{e}}^{\frac{\beta \pi }{2N}}+1\right)}^{2}}-\frac{3{\beta }^{2}{\text{e}}^{\frac{\beta \pi }{2N}}{{\Delta}}_{c}^{\prime 2}}{4{\left({\text{e}}^{\frac{\beta \pi }{2N}}+1\right)}^{3}}\\ & \quad +\frac{{\beta }^{2}{\text{e}}^{-\frac{\beta \pi }{2N}}{{\Delta}}_{\text{c}}^{\prime 2}}{2{\left({\text{e}}^{-\frac{\beta \pi }{2N}}+1\right)}^{2}}-\frac{3{\beta }^{2}{\text{e}}^{-\frac{\beta \pi }{2N}}{{\Delta}}_{\text{c}}^{\prime 2}}{4{\left({\text{e}}^{-\frac{\beta \pi }{2N}}+1\right)}^{3}},\end{align} \tag{ 82 }$$

where the first derivative of the energy gap with respect to g evaluated at g_c = 1 is denoted by

$$\begin{equation}{{\Delta}}_{\text{c}}^{\prime }=\left.\frac{\mathrm{d}}{\mathrm{d}\delta }{\Delta}(1+\delta )\right\vert \delta =0=\left.\frac{\mathrm{d}}{\mathrm{d}g}{\Delta}(g)\right\vert g=1.\end{equation} \tag{ 83 }$$

and we note that terms depending on the second derivative cancel out. The explicit expression for Δ_c' readily follows using (14a)

$$\begin{equation}{{\Delta}}_{\text{c}}^{\prime }=\frac{1}{2}{\Delta}+1=\mathrm{tan}\left(\frac{\pi }{4N}\right)+1\approx \frac{\pi }{4N}+1.\end{equation} \tag{ 84 }$$

Substituting

$$\begin{equation}x(\beta )=\frac{\pi }{4N}\beta \end{equation} \tag{ 85 }$$

and simplifying the expression for R(β), one finds

$$\begin{equation}R(x(\beta ))=\frac{{x}^{2}(\beta )}{16\enspace {\mathrm{cosh}}^{2}\left[x(\beta )\right]}\left(\frac{16}{{\pi }^{2}}{N}^{2}+\frac{8}{\pi }N+1\right).\end{equation} \tag{ 86 }$$

From this expression it is clear that the maximum of the susceptibility scales quadratically with the system size. Moreover, the inverse temperature for which the maximum is achieved scales linearly with the system size. Further, the presence of the maximum is due to the particular behavior of the first derivative of the symmetry breaking gap,

$$\begin{equation}\underset{N\to \infty }{\mathrm{lim}}\enspace {{\Delta}}_{\text{c}}^{\prime }=1{ >}0.\end{equation} \tag{ 87 }$$

Figure 4 shows the range of β where the approximation works. Obviously, the approximation fails in the classical limit associated with large temperatures. For β → 0, both ρ_Gibbs(0, g) and ρ_Gibbs(0, g + δ) are maximally mixed states, proportional to the identity operator. As a result, the fidelity F(0, g|0, g + δ) = 1 for any h and δ. Consequently, the fidelity susceptibility vanishes by definition,

$$\begin{equation}\underset{\beta \to 0}{\mathrm{lim}}\enspace {\chi }_{\text{Exact}}(\beta )=0.\end{equation} \tag{ 88 }$$

In this limit, the TLA naturally fails and predicts incorrectly a finite value

$$\begin{equation}\underset{\beta \to 0}{\mathrm{lim}}\enspace \;{\chi }_{\text{TLA}}(\beta )=\frac{1}{96}\left(2{N}^{2}-3N+1\right),\end{equation} \tag{ 89 }$$

which scales quadratically with the system size and diverges in the thermodynamic limit.

Zoom In Zoom Out Reset image size

The appearance of χ⁻ in (80) indicates the importance of the negative parity sector, which was omitted in previous studies [18, 19].

5.2. Thermal fidelity susceptibility and specific heat

When the thermal states commute with each other, the Uhlmann fidelity is given in terms of the partition function, as we have discussed. The Taylor series expansion of equation (68) yields the thermal fidelity susceptibility, which is itself proportional to the specific heat of the system at constant magnetic field. Specifically, taking β → β − δβ/2 and β' → β + δβ/2, and expanding in δβ, one find the leading term [9]

$$\begin{equation}\xi (\beta ,g)=-2\enspace \underset{\delta \beta \to 0}{\mathrm{lim}}\frac{\mathrm{ln}\enspace F(\beta -\delta \beta /2,g\vert \beta +\delta \beta /2,g)}{{(\delta \beta )}^{2}}=\frac{{C}_{v}(\beta ,g)}{4{\beta }^{2}}.\end{equation} \tag{ 90 }$$

Next, we aim at computing the thermal fidelity susceptibility (and, thus, the specific heat) using the TLA (72). Exact expressions, in this case, can be calculated with help of equation (68) and knowledge of partition function (47). Simplified expressions, used in literature, can be easily obtained just by substitution of the PPA partition function (49) into equation (68). To derive the TLA thermal fidelity susceptibility we use the approximated formula for fidelity between equilibrium states with temperatures differing by δ,

$$\begin{align}{F}_{\text{TLA}}\left(\beta ,g\vert \beta +\delta ,g\right)& =F\left(\beta ,\gamma =1,g\vert \beta +\delta ,\gamma =1,g\right)\\ & =\frac{{Z}_{\text{TLA}}(\beta +\delta /2)}{\sqrt{{Z}_{\text{TLA}}(\beta ,g)\enspace {Z}_{\text{TLA}}(\beta +\delta ,g)}},\end{align} \tag{ 91 }$$

where Z_TLA is given by formula (73). The explicit expression reads

$$\begin{equation}{F}_{\text{TLA}}\left(\beta ,g\vert \beta +\delta ,g\right)=\frac{1+{\text{e}}^{(-\beta +\delta /2){\Delta}(g)}}{\sqrt{1+{\text{e}}^{-\beta {\Delta}(g)}}\sqrt{1+{\text{e}}^{(\beta +\delta ){\Delta}(g)}}},\end{equation} \tag{ 92 }$$

where Δ(g) is symmetry-breaking gap defined in the formula (22). Using it, the calculation of specific heat yields a remarkably simple expression

$$\begin{align}{C}_{v}\left(\beta ,g\right)=4{\beta }^{2}\enspace \xi (\beta ,g)& =-4{\beta }^{2}\enspace \frac{{\partial }^{2}}{\partial {\delta }^{2}}F(\beta ,g\vert \beta +\delta ,g){\vert }_{\delta =0}\\ & =\frac{{\beta }^{2}{{\Delta}}^{2}(g)}{4}\enspace {\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}}^{2}\left[\frac{\beta {\Delta}(g)}{2}\right].\quad (\mathrm{T}\mathrm{L}\mathrm{A}).\end{align} \tag{ 93 }$$

Above formula gives TLA approximation for arbitrary values of β and g, one expects that approximation works better for larger β. In order to channel the discussion, we restrict ourselves to the case g = 1. Then, substituting $x(\beta )=\frac{\beta \pi }{4N}$, one finds

$$\begin{equation}{C}_{v}\left(\beta \right)={C}_{v}(\beta ,g=1)={x}^{2}(\beta ){\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}}^{2}\left[x(\beta )\right].\quad (\mathrm{T}\mathrm{L}\mathrm{A}).\end{equation} \tag{ 94 }$$

The temperature dependence of the specific heat at the critical point is shown in the figure 5, where the exact results are compared with the PPA and TLA. A characteristic feature of the temperature dependence of the specific heat is the big contrast between the high- and low-temperature regimes. In the high-temperature regime, the specific heat exhibits a sharp peak, which scales linearly with the system size. On the other hand, at low temperatures, the TLA formula (94) provides an accurate prediction, with a maximum ${C}_{v}^{\mathrm{max}}\approx 0.439\enspace 229$ independent of the system size. By contrast, the PPA partition function gives results very close to the exact one for high temperatures, as in the case of the susceptibility related to the field g (compare with the previous section). In short, the exact results are well reproduced by the TLA in the low-temperature regime, and by the PPA in the high-temperature regime.

Zoom In Zoom Out Reset image size

An important family of paradigmatic spin models-including the one-dimensional XY and Ising models-are integrable and can be expressed in terms of free fermions. The parity operator commutes with the Hamiltonian and it is often convenient to simplify their description by considering only the positive-parity subspace. This yields simple approximate formulae for relevant quantities such as the partition function which are ubiquitous in the literature [2, 3].

We have shown that such an approach fails in the characterization of quantum critical phenomena using the ground-state fidelity and fidelity susceptibility. Using an algebraic approach to exactly account for the complete Hilbert space [25], we have shown that the exact result for the fidelity susceptibility between thermal states qualitatively differs from the conventional approximation, away from the classical high-temperature regime. The discrepancy is pronounced in the quantum, low-temperature limit when the accurate treatment of the low-lying energy states is crucial. Furthermore, the discrepancy is robust against variations of the system size and is manifested even in the thermodynamic limit.

Our results show the limitation of disregarding the odd parity subspace in the characterization of integrable spin chains at finite temperature and are potentially relevant to applications ranging from quantum thermodynamics to parameter estimation, among other quantum technologies using quantum critical spin chains. The exact expressions we have provided for the fidelity between thermal states should be of broad interest in the quantum-information characterization of these systems at finite temperature, with applications ranging from quantum thermodynamics to many-body and criticality-enhanced quantum metrology [15, 40, 41]. Likewise, the exact treatment of the fidelity susceptibility is required for the analysis of their critical phenomena and could be used for benchmarking quantum simulators, annealing devices, and the performance of quantum algorithms using quantum spin chains as a testbed.

Our results apply to the canonical Gibbs state at thermal equilibrium resulting from thermalization dynamics without other conserved quantities than the energy in chains of fixed size. An interesting outlook concerns the characterization of the fidelity and fidelity susceptibility in the generalized Gibbs ensemble, i.e. in the presence of additional invariants of motion.

MB would like to thank Bogdan Damski and Andrzej Syrwid for insightful discussions. We further thank Marek M Rams for feedback on the manuscript. MB also acknowledges the support of the Polish National Science Center under scholarship ETIUDA (2020/36/T/ST3/00332) and research Grant DEC-2016/23/B/ST3/01152. F J G-R thanks the University of Luxembourg for hospitality during the completion of this work. This work is further supported by the Spanish Ministerio de Ciencia e Innovación (PID2019-109007GA-I00).

The data that support the findings of this study are available upon reasonable request from the authors.

Numerical simulations were performed with Wolfram Mathematica, version 12. Fidelity susceptibility from the formula (47) is calculated first by symbolic differentiation to obtain analytical formula and then insertion of the parameters β and g—facilitating a very high precision. In benchmarking the analytical expressions with the numerically-exact diagonalization, it is required to approximate the second derivative. We use central finite difference. In order to achieve satisfactory precision we have to use the fourth order according to which the second derivative of arbitrary function f(x) is given by [42]

$$\begin{align}\frac{{\mathrm{d}}^{2}}{\mathrm{d}{x}^{2}}f(x)& \approx \frac{1}{{\Delta}{x}^{2}}\left(-\frac{1}{12}f(x-2{\Delta}x)+\frac{4}{3}f(x-{\Delta}x)-\frac{5}{2}f(x)\right.\\ & \quad +\left.\frac{4}{3}f(x+{\Delta}x)-\frac{1}{12}f(x+2{\Delta}x)\right).\end{align} \tag{ A.1 }$$

This formula can be obtained by interpolating the function f with a polynomial of degree four and can be proven to approximate the second derivative in Δx⁴ order [42]. The fidelity in exact diagonalization is computed directly from the definition with help of built-in Mathematica functions MatrixExp and MatrixPower.

Figure A1 shows that the numerically-exact results match the exact analytical expressions. Moreover, plots for different system sizes N indicate that the position of the maximum approaches the thermodynamic value at g_c = 1 already for N = 10 (see also figure 2). The height of the maximum grows with the system size like N², as shown in section 5.

Zoom In Zoom Out Reset image size