Locality of temperature in spin chains

In traditional thermodynamics, temperature is a local quantity: a subsystem of a large thermal system is in a thermal state at the same temperature as the original system. For strongly interacting systems, however, the locality of temperature breaks down. We study the possibility of associating an effective thermal state to subsystems of infinite chains of interacting spin particles of arbitrary finite dimension. We study the effect of correlations and criticality in the definition of this effective thermal state and discuss the possible implications for the classical simulation of thermal quantum systems.

Abstract. In traditional thermodynamics, temperature is a local quantity: a subsystem of a large thermal system is in a thermal state at the same temperature as the original system. For strongly interacting systems, however, the locality of temperature breaks down. We study the possibility of associating an effective thermal state to subsystems of infinite chains of interacting spin particles of arbitrary finite dimension. We study the effect of correlations and criticality in the definition of this effective thermal state and discuss the possible implications for the classical simulation of thermal quantum systems.

Introduction
The question whether the standard notions of thermodynamics are still applicable in the quantum regime has experienced a renewed interest in the recent years. This refreshed motivation can be explained as the consequence of two successes. On the one hand, the spectacular progress of the experiments encompassed in the so called quantum simulators already allows for a direct observation of thermodynamic phenomena in many different quantum systems, such as ultra cold atoms in optical lattices, ion traps, superconductor qubits, etc [1,2,3,4]. On the other hand, the inflow of ideas from quantum information theory provided significant insight into the thermodynamics of quantum systems [5,6,7,8].
Specifically, qualitative improvements have been made in understanding how the methods of statistical mechanics can be justified from quantum mechanics as its underlying theory [2,6,9,10,11,12]. One of the fundamental postulates of thermodynamics is the so called Zeroth Law: two bodies, each in thermodynamic equilibrium with a third system, are in equilibrium with each other [13,14]. This is the law that stands behind the notion of temperature [13,14]. In fact, the above formulation of the Zeroth law consists of three parts: (i) there exists a thermal equilibrium state which is characterized by a single parameter called temperature, and isolated systems tend to this state [2,6,9,10]; (ii) the temperature is local, namely, each part of the whole is in a thermal state [13,14]; and (iii) the temperature is an intensive quantity: if the whole is in equilibrium, all the parts have the same temperature [13,14,15,16,17,18]. The last two points are usually derived from statistical mechanics under the assumption of weakly interacting systems. Nevertheless, when the interactions present in the system are non-negligible, the points (ii) and (iii) need to be revised. Following the direction given by Refs. [15,16], in this work we concentrate on the clarification and generalization of the aforementioned aspects of the Zeroth Law of the thermodynamics for spin chains with strong, short range interactions. The general setting of the problem is as follows. The system with Hamiltonian H is in thermal equilibrium, described by a canonical state at inverse temperature β, and we seek to understand the thermal properties of a finite part of the system. Obviously, in the presence of strong interactions, the reduced density matrix of a subsystem of the global system (especially in the quantum regime [18]) will not generally have the same form as (1). In lattice systems, where the Hamiltonian is a sum of local terms interacting according to some underlying graph, it is unclear how one can locally assign temperature to a subsystem. More precisely, the reduced density matrix of the subsystem A (see Fig. 1) of a global thermal state is described by which will not be thermal unless the particles in A do not interact with its environmentĀ. Hence, given only a subsystem state ρ A and its Hamiltonian H A , it is not possible to assign a temperature to it, since this would totally depend on the features of the environment and the interactions that couple the subsystem to it.
In the context of quantum information, a first step to circumvent the problem of assigning temperature to a subsystem was made in Ref. [15]. There, for harmonic lattices it was shown that it is sufficient to extend the subsystem A by a boundary region B that, when traced out, disregards the correlations and the boundary effects (see Fig. 1). If the size of such a boundary region is independent of the total system size, temperature can still be said to be local. More explicitly, given a lattice Hamiltonian H with a subsystem A, a shell region around it B and its environment C = (A ∪ B) c , see Fig. 1 where ρ A = Tr B ω(H AB ) is the state of A for the chain truncated to AB, and f ( B ) is expected to be a monotonically decreasing function in B . The width of the boundary region B is defined as the graph-distance between the sets of vertices (regions) A and C. Surely, the differences (3) fully characterize the distance of ρ A from ρ A . Indeed, the trace distance, a distinguishability measure for quantum states, D(ρ A , ρ A ) = 1 2 ||ρ A − ρ A || 1 , has the following representation [19,20]: where I is the identity operator in the Hilbert space of A. In Ref. [17], it is proven that the correlations responsible for the distinguishability between the truncated and non-truncated thermal states are quantified by a generalized covariance. For any two operators O and O , full-rank quantum state ρ, and parameter τ ∈ [0, 1], the generalized covariance is defined as and the average distinguishability of the two states measurements of some observable O can provide reads as where H I is the corresponding Hamiltonian term that couples B and C, H 0 = H − H I is the truncated Hamiltonian (see Fig. 1) and ω s = ω(H(s)) is the thermal state of the interpolating Hamiltonian H(s) := H − (1 − s) H I . Hence, the generalized covariance is the quantity that measures the response in a local operator of perturbing a thermal state and ultimately at what length scales temperature can be defined.
Temperature is known to be a local quantity in a high temperature regime. More specifically, in Ref. [17], it is shown that for any local Hamiltonian there is a threshold temperature (that only depends on the connectivity of the underlying graph) above which the generalized covariance decays exponentially. Nevertheless, it is far from clear what occurs below the threshold, and, especially, at low temperatures (β 1). Note that, in that case, the right hand side of the truncation formula (6) could be significantly different from zero since the integration runs up until β/2, while the covariance is expected to decay only algebraically for critical systems. In this work we show that, for one dimensional translation-invariant systems, temperature is local for any β. Away from criticality, we rigorously bound the truncation formula (6) by mapping the generalized covariance to the contraction of a tensor network and exploiting some standard results in condensed matter. At criticality, we use some results from conformal field theory [21,22]. Finally, the results in [23], where the equivalence of microcanonical and canonical ensembles is proven for translation-invariant lattices with short range interactions, render our results valid also when, instead of being canonical, (1), the global state ω(H) is, e.g., microcanonical. The latter is defined as an equiprobable mixture of all the energy eigenstates in a narrow energy window (see [23] for details). In condensed matter physics, this problem has been considered in the context of approximating the expectation values of infinite systems by finite ones, receiving the name of finite size scaling (see, e.g., [24] and references therein). Nevertheless, the finite-sizescaling methods are more focused to find the values for the critical exponents and the transition temperature by observing how measured quantities vary for different lattice sizes. The paper is organized as follows. In Sec. 2, the generalized covariance of 1D system is mapped to the contraction of a 2D tensor network. In Sec. 3, we show that temperature is local at non-zero temperature (β < ∞) by identifying in the tensor network a gapped transfer matrix which leads to a clustering of correlations and ultimately to a clustering of the generalized covariance. In Sec. 4, locality of temperature is proven at zero temperature (β → ∞) using different methods for gapped and gapless systems. While transfer matrix arguments work satisfactory for gapped systems, conformal field theory results have to be used at criticality. In Sec. 5, all our results are illustrated in detail for the Ising model, for which we study in addition the behaviour of the generalized covariance and compute explicitly the physical distinguishability between the full and truncated Hamiltonians both at and off criticality. Finally, we conclude. Let us consider a system of spins described by a short range Hamiltonian. The structure of the Hamiltonian is given by a graph G(V, E). The spins correspond to the set of vertices V and the two-body interactions to the edges E. Such a Hamiltonian can be written as where h u are the Hamiltonian terms acting non-trivially on the adjacent vertices of u. In Ref. [26,27], it is shown that, for any error ε > 0, the matrix e −βH of a local Hamiltonian can be approximated in one norm by its Trotter-Suzuki expansion, such that where m > 360β 2 |E| 2 /ε and the products over u and v in Eq. (8) are realized in the same order.
To illustrate the previous approximation, let us consider in detail the one dimensional case: a spin chain with nearest neighbour interactions. By decomposing in the standard way the Hamiltonian in its odd and even terms, the tensor networkρ T N becomes in this casẽ where H odd/even = u∈odd/even h u and H = H odd + H even . Let us think about each exp(β/mh u ) as a tensor. In this way,ρ T N can be seen as the contraction of several of such tensors, that is, a tensor network see Fig. 2 (c). Starting from a one dimensional quantum system,ρ T N can be interpreted as a tensor network spanning two dimensions, with the extra dimension of length m. We will refer to this extra dimension as the β direction, while the original dimension will be called spatial direction. In Fig. 2 (a), the diagrammatic representation ofρ T N is presented. Its tensors can be decomposed and arranged in order to form a square lattice of elementary tensors as shown in Fig. 2 . Diagramatic representation of (a) the expectation value of a one site operator and (b) the generalized covariance (a two-point correlation function) between two one site operators. In both cases, the final result is computed as the ratio between the contraction of two tensor networks.

Generalized covariance as the contraction of tensor networks
The expectation value of a local operator is given by By using Eq. (9), the fact that e −βH 1 = Z and some elementary algebra, the expectation value of a local operator A can be approximated by the ratio between the contraction of two tensor networks This is represented diagramatically in Fig. 3(a). The generalized covariance can be rewritten as whereÕ = O − Tr(Oω(H)) for any operator O. Hence, in a similar way as it has been made for the expectation values, the generalized covariance can be also approximated as the ratio between two tensor network contractions as shown in Fig. 3 From this perspective, the generalized covariance can be seen as a two point correlation function on a 2 dimensional lattice in which τ m is the separation in the β direction and the distance between the non-trivial supports of O and O is the separation in the spatial direction (see Fig. 3). This construction can be generalized to approximate expectation values of local operators and n-point correlation functions of a D dimensional quantum model by the ratio of the contraction of two D + 1 dimensional tensor networks.

Transfer matrices
It is also very useful to define two extra objects: the transfer matrices along the spatial T and β directions T β . The first is obtained by contracting a column of the elementary tensors of the network, while the second is obtained by contracting several rows of elementary tensors, see Fig. 4. The number of rows that need to be contracted in order to obtain the transfer matrix in the β direction, T β , is chosen such that its spectral gap between the largest and second-largest eigenvalues is independent of both β and m. This can be achieved by contracting m/β rows, leading to a transfer matrix with two largest eigenvalues µ 1 and µ 2 where ∆ is the gap of the Hamiltonian.

Locality of temperature at non-zero temperature
Let us consider now the case in which β is of order one. The physical distinguishability in A between the full and the truncated Hamiltonians can be bounded by where H I = h L +h R , i ∈ {L, R}, and we have used the linearity of the generalized covariance with respect of its operators. Without loss of generality let us assume that the term of H I that maximizes the generalized covariance is the one of the left, h L . Hence, the quantity to bound is cov τ ωs (h L , O). In order to do so, let us rewrite it as where Z s := 1 L |T s T 2 B +1 T s |1 R is the partition function, T is the transfer matrix in the spatial direction, see Fig. 4 (left), and |1 L/R is the left/right dominant eigenvector of T i. e. the eigenvector associated to its largest eigenvalue. The matrix T s is the transfer matrix corresponding to the boundaries BC where the elementary tensors of the network are different from the rest for s < 1 and T 1 = T . The matrix Y corresponds to the slice of the region A where the operator O is supported, and the matrix X s is the transfer matrix T s with the insertion of the operators h L located at a distance τ β from O in the transverse direction. The diagrammatic representations of the matrices X s , Y and T s are shown in Fig. 5. To simplify the calculations, the transfer matrix can always be normalized such that its dominant eigenvalue is λ 1 = 1.
To bound the generalized covariance (17) it is useful to rewrite it in terms of 2-point correlation functions of the uniform system (s = 1) and we have used that In short range one dimensional systems, the absence of phase transitions at non-zero temperature [25] implies that the transfer matrix T is gapped, with a gap related to the spatial correlation length as where λ 2 is the second largest eigenvalue of the transfer matrix T . For gapped transfer matrices, the 2-point correlation function (20) can be proven to be upperbounded by The complete proof of the previous statement can be found in Lemma 1 of the Appendix A. Furthermore, Lemma 1 allows us to bound all the terms in Eq. (18), and, as it is shown in Appendix B, the following inequality holds for the left hand side of Eq. (16): where The quantity c is a constant of order one that depends on the model considered. Hence, the temperature is proven to be intensive for any one dimensional translationally invariant model at non-zero β.

Gapped systems
Given a Hamiltonian with gap ∆, here we study the regime in which β −1 ∆. This implies that the lattice in its β direction is much larger than the correlation length β ξ β , with In the limit of temperature tending to zero, the 2D network that represents the partition function becomes infinite in the β direction (see Fig. 5).
In order to see that the temperature is also local in this case, let us decompose the integral over t of the generalized covariance into two pieces where L is a cut-off that will be chosen afterwards to minimize a bound on the right hand side, and β will be made to tend to infinity. Concerning the integral over 0 ≤ t ≤ L, we will exploit the fact that the system is gapped, and hence its ground state is known to have a finite correlation length ξ in the spatial direction and to be represented by a Matrix Product State of bond dimenson D, with D ∝ poly(ξ) [40,41,42]. As argued in the previous section, a finite correlation length guarantees a gap in the transfer matrix in the corresponding direction. By performing an analogous calculation to the one described in the previous section, one obtains The second integral over t > L can be bounded by taking the transfer matrix in the β direction which is also gapped for gapped Hamiltonians. More specifically, the generalized covariance can be written as where we have identified T β as the transfer matrix for which the ground state of the Hamiltonian |GS is its dominant eigenvector. As previously, we make use of Lemma 1 in the Appendix A and obtain The integration is then bounded by Putting the previous bounds together, and after an optimization over L, we get showing that temperature can be locally assigned to subsystems for arbitrarily large β and gapped Hamiltonians.

Criticality
A system at zero temperature is said to be critical when the gap between the energy ground state (space) and the first excited state closes to zero in the thermodynamic limit. The critical exponents zν control how the spectral gap ∆ tends to zero where N is the system's size and ν is the critical exponent that controls the divergence of the correlation length The previous divergences are a signature of the scale invariance that the system experiences at criticality. If the critical exponent z = 1, there is a further symmetry enhancement and the system becomes conformal invariant.
up to higher order terms, where y is the scaling dimension of the operator H I [38,39]. If H I is a standard Hamiltonian term, in the sense that the system is homogeneous, its leading scaling dimension is y = 2.
Once more, we see that by increasing the buffer region temperature can be arbitrarily well assigned.

A case study: The Ising chain
Now we illustrate our results for the quantum Ising chain, which is described by the Hamiltonian where σ i x and σ i z correspond to the Pauli matrices, h characterizes the strength of the magnetic field and N is the number of spins. Notice that the interactions in the above Hamiltonian are of finite range, a crucial assumption in our derivations, see (7). This model has a quantum phase transition at h = 1, so it well exemplifies the different regimes discussed above: criticality (only at zero temperature) and away from it (for zero and non-zero temperatures).

Generalized Covariance
First of all, as in the previous sections, we split the chain in three regions, which are shown in Fig. 6. For such a splitting, and in the context of Eq. (6), we compute the generalized covariance cov    , H I ) as a function of t/β for several temperatures (β = 5, 20, 1000), and for h = 0.9, 1 (i.e., near and at criticality). We take N = 40, which already describes well the thermodynamic limit (recall that we are only interested on the local state, and that the correlations decay exponentially). The area below the curves correspond to the first integral in (6), which measures how well the local state in A can be approximated by a thermal state in AB. The results in Fig. 7 are in agreement with properties (i) and (ii) from Lemma 2 in Appendix A. The first property implies that the covariance is symmetric with respect to t = β/2, and it follows by taking l = t and n = β in (A.8). Second, property (ii) implies that it is bounded by a convex function of t with a maximum at t = 0 and t = β and with a minimum at t = β/2. Therefore, the covariance satisfies the bound (A.9). On the other hand, the covariance is not monotonic in s (see Fig. 8). This is somehow counterintuitive, as it shows that the outcomes of two observables with no overlapping support (located in A and in the intersection between B and C) do not always become more correlated as s, which quantifies the strength of the interaction between B and C, increases.  (right). The grey area below the curves corresponds to the second integral of Eq. (6). Notice that, due to the symmetry in t, the values t/β = 2/3, 1 are also considered.

Locality of temperature in the quantum Ising chain
In our analytical findings, the generalized covariance naturally appeared as a tool to solve the locality of temperature problem, see (6). This motivated the previous section, where we studied its properties in the context of the quantum Ising chain. Nevertheless, in order to obtain (6), one still needs to integrate cov t/β ωs (O, H I ) over s and τ . While this approach is useful when dealing with arbitrary generic systems, here we are dealing with a specific model that is furthermore solvable, so we can take a more direct approach. Concretely, we first compute and for different sizes N of the region AB. Secondly, we measure the distinguishability between such states via the quantum fidelity, which is advantageous for computational reasons. For two states, ρ A and ρ A , the fidelity is defined as [19] F It satisfies 0 ≤ F ≤ 1 and F [ρ A , ρ A ] = 1 if and only if ρ A = ρ A . In order to relate this approach to our previous considerations, we note the following relation between the trace distance, D[ρ A , ρ A ], and F (ρ A , ρ A ), given in [20], Therefore, the fidelity provides us with upper and lower bounds to (4). In particular, when , and in that case we say that the temperature is locally well defined. From now on, we take for A a two spin subsystem, an infinite chain as the total system, and we compute F (ρ A , ρ A ) as a function of the size of AB, with N = 2 + 2l B , and the different parameters of the Hamiltonian. In order to compute ρ A and ρ A , it is convenient to apply the Jordan-Wigner transformation to (34), which maps spin operators σ i x,y,z to fermionic operators a i , a † i (see Appendix C for details). The Hamiltonian (34) takes then the form, which is quadratic in terms of the fermionic operators. It follows that thermal states, as well as their local states, are gaussian operators. Therefore it is possible to describe them by their covariant matrix, whose size is only O(N 2 ). This allows us to compute ρ A in (37) for finite but large l B ; while in the limit N → ∞, i.e. to compute ρ A in (36), we rely on the analytical results from [36]. The explicit calculations are done in Appendix C. Figure 9 shows F (ρ A , ρ A ) as a function of β and h, for N = 4 (left) and N = 20 (right). Recall that N , with N = 2 + 2l B , defines the size of the boundary region which is used to approximate ρ A by ρ A . Even if the boundary is small, N = 4, the fidelity is close to 1 for all values of β and h, and thus the temperature is locally well-defined. As expected, F (ρ A , ρ A ) increases with N (see Fig. 10). We also observe in Fig. 9 that the fidelity becomes minimal near h = 1, which is the phase transition point. As N increases, this minimum is shifted to h = 1. At this point the spatial correlations also increase, which suggests a relation between both quantities. In order to further explore this connection, we compute the scaling of F (ρ A , ρ A ) with N , and compare it to the decay of the correlations. The behaviour of F (ρ A , ρ A ) is plotted in Fig. 10, which clearly shows that the fidelity follows an exponential law with N , given by

Non-zero temperatures
where ξ S is a parameter that characterizes the slope of the function. On the other hand, the correlations between a local observable in A, σ i z , and one in the intersection of B and C, σ i+d z , can be obtained through the two-spin correlation function σ i z σ i+d z in [36]. Their asymptotic behaviour is also exponential with d, where ξ is the correlation length. Now, identifying d, the distance between particles, with N/2, which is roughly the size of B, we obtain from the numerical results in Fig. 11 the following simple relation, Roughly speaking, the quality of the approximation ρ A is directly related to the strength of the correlations in the system. This relation is in good agreement with previous considerations in [43], where the correlation length is related to the error of the cluster approximation [43,44].
In summary, temperature can be assigned to the local system for all h and non-zero β by taking a small boundary region (with N ≥ 4, and thus l B ≥ 2). We have shown that this is directly connected to the exponential decay of the correlations with the distance, which makes the local state of a thermal state only be sensible to its closest boundary.

Absolute zero temperature
The same conclusions apply at zero temperature, as the fidelity is also close to 1 for all h and N ≥ 4. It also has a minimum near the critical point. Nonetheless, the scaling of the fidelity (or more precisely 1 − F ) with N can differ from the previous case. While the scaling is generally exponential at zero temperature, it becomes a power law at the phase transition point (see Fig. 12), This type of decay is also obtained for the correlations as a function of the distance, which again shows a direct connection between the quality of the approximation (quantified by  F (ρ A , ρ A )) and the strength of the correlations.

Conclusions
In this work we studied the locality aspect of the Zeroth law of thermodynamics for quantum spin chains with strong but finite range interactions. Upon noting that in the presence of strong interactions the marginal states of a global thermal state do not take the canonical form themselves, we go on defining an effective thermal state for a subsystem. The latter being the reduced density matrix of the subsystem considered as a part of a slightly bigger, enveloping thermal system (see Fig. 1). Borrowing concepts from quantum information theory and employing methods from quantum statistical mechanics, we relate the accuracy with which the effective thermal state describes the actual state of the subsystem to the correlations present in the whole system (see Eqs. (4,5,6) and the discussion around them). We further utilize a Trotter approximation formula [26,27] to build a tensor network representation of the corresponding states of the subsystem to provide upper bounds on the aforementioned accuracy, depending on the size of the enveloping thermal system, and such physical quantities as the spectral gap of the global hamiltonian and the temperature of the parent chain. At the quantum critical point, we use already existing asymptotical formulas from the conformal field theory. Lastly, we exemplify our analytical findings by analyzing the quantum Ising chain. The latter is complex enough to have a quantum phase transition point, but simple enough to allow for an exact diagonalization by standard tools of statistical mechanics, thereby serving as a perfect testbed for our analytical upper bounds. In particular, we find that, e.g., away from criticality, the envelope which is bigger than the system only by one layer of spins, is enough to approximate the actual state with a rather high precision (see, e.g., Fig. 9).
Our results for one dimensional systems with finite range interactions suggest that investigating the properties of the effective thermal states in higher dimensions and, possibly, harbouring long range of interactions, is an interesting direction for further research, which can have far-reaching implications in efficient simulation of the subsystems of large and strongly interacting quantum systems. Another interesting open question beyond the scope of this work is whether these results can be generalized to other other types of equilibrium states, e. g. the so called Generalized Gibbs Ensemble and steady states of local Liouvillians. In a more practical vein, another field where our findings may find implications is quantum thermometry with non-negligible interactions [45].
Let us now consider second factor separately. By inserting a resolution of the identity, a straight forward calculation leads to where we have used that |λ 2 | is an upper-bound for all the |λ k | with k ≥ 2 and the Parseval inequality. Finally, we put everything together and get where the 2nd largest eigenvalue |λ 2 | has been written in terms of the correlation length ξ. By introducing the transfer matrix in its spectral representation, f ( ) can be written as where c k = ( 1|O|k k|O |1 ) and d kk = ( k|O|k k |O |k ). Note now that where the correlation length ξ k is defined as Note that as the eigenvalues of the transfer matrix are ordered, a larger k implies a shorter correlation length ξ k .
In a similar way, we can also simplify the terms in the last sum in Eq. (A.12). Note that where the length ξ kk has been defined as ξ −1 kk = ξ −1 k − ξ −1 k . Puting the previous steps together, we get The rest of terms in Eq. (18) can be analogously bounded. Note that they will only contribute to the second order. Putting everything together in Eq. (16), the physical distinguishability on the region A between the truncated and untruncated thermal states is upper-bounded by where c = 1 + 1 0 dsσ L (s) is a constant of order one, depending on the model.

Appendix C. Solving quantum Ising model
In this appendix we find the states (36) and (37) using formalism of covariance matrices.

Jordan-Wigner transformation
Let us first apply the Jordan-Wigner transformation, and where a i and a † i denote annihilation and creation operators, respectively. From this form of the Hamiltonian, we notice it is quadratic, and thus the thermal state (and their marginal states) are gaussian states. Therefore we can deal with them using the covariance matrix formalism.

The correlation matrix
In this formalism, we define the global correlation matrix, Γ, as where ... N ×N refers to a N × N matrix. Given Γ, we can obtain the correlation matrix corresponding to a reduced state by just selecting the corresponding matrix elements of Γ. For example, the correlation matrix of the fermions k, k + 1 is given by, Since the Jordan Wigner transformation is local, in the sense that it maps the kth fermion to the kth spin in the chain, this correlation matrix also corresponds to the two-spin subsystem at sites k and k + 1. This subsystem is precisely the region of interest A in section 5.2, and thus (C.3) corresponds to the correlation matrix of ρ A in (36). Given the reduced correlation matrix, the explicit form of ρ A can be easily obtained. As the reduced state of a thermal state is gaussian, there is a one to one connection between (C.3) and ρ A . Indeed, for any gaussian state, with with M a coefficient matrix, (C. 4) it is straightforward to prove that, provided that M is diagonalizable, (C. 5) or, equivalently, that

Explicit computation
Now we explicitly compute (C.3) for a finite and an infinite chain, in order to obtain ρ A and ρ A , respectively, using relation (C.6).

• Finite chain
For the case of a finite chain, we need to obtain the correlation matrix (C.10) corresponding to the global state. It is then useful to first diagonalize the Hamiltonian (C.1) by applying the Bogoliubov transformation where φ and ψ are real matrices and verify N k=1 φ 2 jk = N k=1 ψ 2 jk = 1. The Hamiltonian can then take the form, where ξ k are the fermionic excitation energies and b k and b † k denote annihilation and creation operators, respectively. The excitation energies, ξ k , and the matrices φ and ψ are obtained by solving the equation where D is a diagonal matrix whose entries correspond to the excitation energies, ξ k .
Once the Hamiltonian is diagonalized, it is easy to compute the correlation matrix of a thermal state at inverse temperature β in the diagonalized basis, obtaining where the non-zero matrices are diagonal.
From that expression we can obtain the correlation matrix in the original basis, Γ(X), via where T is the transformation matrix defined by the Bogoliubov transformation (C.7). That is, Y = T X, with T = γ µ µ * γ * . (C.12) and γ = 1 2 (φ + ψ) and µ = − 1 2 (φ − ψ). (C.13) • Infinite chain (N → ∞) In the case of an infinite chain, (C.3) can be obtained relying on the analytical results from [36]. The partial state of a two-spin subsystem is l=x,y,z σ k l σ k+1 l σ k l ⊗ σ k+1 l , (C.14) where the average σ k z and the two-spin correlation functions { σ k l σ k+1 l } l={x,y,z} are given by [36]. In order to express the state in the fermionic basis, we can compute the reduced correlation matrix (C.3) from this state, 15) with α = σ k z , β = σ k x σ k+1 x and γ = σ k y σ k+1 y .