Multipartite Entanglement at Finite Temperature

The interplay of quantum and thermal fluctuations in the vicinity of a quantum critical point characterizes the physics of strongly correlated systems. Here we investigate this interplay from a quantum information perspective presenting the universal phase diagram of the quantum Fisher information at a quantum phase transition. Different regions in the diagram are identified by characteristic scaling laws of the quantum Fisher information with respect to temperature. This feature has immediate consequences on the thermal robustness of quantum coherence and multipartite entanglement. We support the theoretical predictions with the analysis of paradigmatic spin systems showing symmetry-breaking quantum phase transitions and free-fermion models characterized by topological phases. In particular we show that topological systems are characterized by the survival of large multipartite entanglement, reaching the Heisenberg limit at finite temperature.

The interplay of quantum and thermal fluctuations in the vicinity of a quantum critical point characterizes the physics of strongly correlated systems. Here we investigate this interplay from a quantum information perspective presenting the universal phase diagram of the quantum Fisher information at a quantum phase transition. Different regions in the diagram are identified by characteristic scaling laws of the quantum Fisher information with respect to temperature. This feature has immediate consequences on the thermal robustness of quantum coherence and multipartite entanglement. We support the theoretical predictions with the analysis of paradigmatic spin systems showing symmetry-breaking quantum phase transitions and free-fermion models characterized by topological phases. In particular we show that topological systems are characterized by the survival of large multipartite entanglement, reaching the Heisenberg limit at finite temperature.
Current studies on entanglement in strongly-correlated systems 1-3 have mainly focused on bipartite and pairwise entanglement 27 . This is, however, clearly unsuited to capture the richness of multiparticle correlations and hardly accessible experimentally in systems of a large number of particles 28 that are the natural targets of quantum simulators 29,30 . Much less attention has been devoted to witnessing multipartite entanglement [31][32][33][34][35][36][37] and this has been mainly limited to spin models. While only few witnesses are known in the literature 38 , multipartite entanglement up to hundreds/thousands of spins has been successfully detected experimentally in atomic ensembles 39 . Among these witnesses, the quantum Fisher information (QFI) has proved to be especially suitable [39][40][41][42][43] and it is currently attracting considerable interest 37,[44][45][46][47][48][49][50][51] . The QFI has an appealing operational meaning in terms of statistical speed of quantum states under external parametric transformations 44,45 , it extends the class of states detectable by popular methods such as the spin squeezing 40,44,[52][53][54] , and it can witness entanglement in spin systems 37,47,55 as well as in free-fermion topological models 48,49 . Furthermore, the QFI can be extracted experimentally using a statistical distance method 44,45 , or by a weighted integral of the dynamic susceptibility across the full spectrum 37 . Measurable lower bounds to the QFI have been extracted experimentally 44,53,54 and proposed theoretically [56][57][58] .
Q plays a central role in the theory of quantum coherence [59][60][61][62][63] : it quantifies the coherent extent of a generic state ρ over the eigenstates of the operator Ô , vanishing if and only if ˆρ = O [ , ] 0. Multipartite entanglement is witnessed when the QFI overcomes certain finite bounds: as discussed below 41,42 , Q is only achievable if ρ contains (κ + 1)-partite entanglement among N parties and Ô is a local operator.
In this manuscript we show that the QFI of a many-body system at thermal equilibrium in the vicinity of a quantum critical point λ c has the universal behavior shown in Fig. 1 Here, ρ T is the thermal state at temperature T (here and in the following the Boltzmann constant is set to 1), μ and ν indicate the degeneracy of the ground state of energy E gs and first excited state of energy E ex , respectively, and Δ = E ex − E gs is the first energy gap in the many-body spectrum. Equation (1) is valid for  Δ T and shows that, regardless on the microscopical details of the system, the lower bound to ˆρ F O [ , ] Q T factorizes in a thermal and a quantum contribution. The thermal decaying function on the right side of the inequality (1) only depends on the structure of the low-energy spectrum, i.e. the energy gap and the degeneracy of the energy eigenstates. The bound is tight for Q 0 is the zero-temperature limit of the QFI and depends whether the ground state is degenerate or not.
If the ground state is nondegenerate (μ = 1), given by the pure state |ψ 0 〉, a Taylor expansion of the right-hand side of Eq. (1) givesˆˆˆ Q T Q T T 0 / / 2 and shows that the QFI is bounded from below by a constant for  T T cross , where T cross ≈ Δ/log (3 + ν). This defines a quantum plateau (QP) where the zero-temperature QFI, , is insensitive to thermal fluctuations, being protected by the finite energy gap Δ. In particular, if |ψ 0 〉 hosts multipartite entanglement witnessed by the QFI, such multipartite entanglement is robust against temperature for  T T cross . In the following we provide examples of systems characterized by large multipartite entanglement in the ground state (even approaching the Heisenberg scaling at finite T, see Sec. IV) that is insensitive to small temperatures.
Whenever the ground state is degenerate (μ > 1), in the limit T → 0 the QFI is given by , where ρ 0 is the incoherent mixture of the μ degenerate ground states, see Sec. II. According to Eq. (1),ˆ Also in this case, the lower bound remains constant for  T T cross , where T cross ≈ Δ/log(3 + ν/μ). If the ground state becomes degenerate only in the thermodynamic limit, this constant value defines a thermal plateau (TP) where thermal fluctuations strongly affect the QFI of the (pure) ground state |ψ 0 〉 outside the thermodynamic limit, but not the QFI of the incoherent mixture ρ 0 . In other words, the QFI of the ground state Q 0 may be very high -|ψ 0 〉 being given for instance by a maximally entangled state -but it exponentially decays with temperature to a much smaller value ˆρ Q 0 that remains constant up to T cross . In Fig. 1 we schematically plot the case of a typical symmetry-breaking model, where the TP (matching the ordered phase) and the QP (matching the disordered phase) are found on different sides of the critical point. Examples of symmetry-breaking models Control parameter λ versus temperature T for the QFI of a critical many-body system. We distinguish four regions depending on the scaling exponent β = d log F Q /d log T of the QFI with respect to temperature: a quantum plateau (QP), a thermal plateau (TP), a critical plateau (CP) and a maximum entropy plateau (MEP). QP and TP are defined from the lower bound Eq. (1), showing that the QFI remains at least constant (β ≥ 0) up to a crossover temperature T cross (white solid line) of the order of the first nonvanishing gap Δ in the energy spectrum (dashed line). The characteristic feature of the TP region is the degeneracy of the ground state: in the thermodynamic limit, the QFI suddenly decreases from its value at T = 0 to the plateau value. In the CP, the QFI follows a scaling law controlled by critical exponents of the model, β = −Δ Q /z, according to Eq. (4). For temperatures larger than T max (dotted line) -approximatively equal to the maximum energy of the spectrumthe QFI enters the MEP where β = −2. In the crossover grey regions the thermal decay is non-universal.
will be discussed in more details in Sec. III. In the absence of ground-state degeneracy, the TP is absent and the QP is found on both sides of the critical point. This behavior is found for topological closed chains, as shown in Sec. IV. At finite temperature and for values of λ around the critical point λ c , a scaling hypothesis for the dynamical susceptibility 37 predictsˆρ Here, N is the total number of parties in the system (e.g. the total number of spins), Δ Q is the exponent 37 that characterizes the finite-size scaling of the QFI with respect to N at Q , and z is the dynamical critical exponent. We thus identify a region of parameters in the vicinity of the critical point (T > 0) that we call critical plateau (CP) where the QFI follows the scaling behavior Eq. (4) as a function of temperature. In general, we expect that the CP extends for λ λ | − | ν  T z c , where ν is the correlation-length critical exponent. This region matches a quantum critical regime 4,16 where the scaling behavior of a quantum coherence measure, the QFI, at finite temperature is controlled by critical exponents of the transition. In Fig. 1 the CP is schematically represented as a triangular region. The CP is separated from the TP and QP by a model-dependent smooth decay for ≈ T T cross . Finally, for temperatures of the order of the interaction energy scale of the system, no multipartite entanglement is witnessed by the QFI. Moreover, for temperatures larger than the maximum energy of the spectrum, the QFI decays asˆρ This defines a fourth plateau that we identify as maximum entropy plateau (MEP). In this regime, all eigenstates are approximatively equally populated.
It is worth clarifying that the operator Ô in Eqs (1)(2)(3) and (5) is arbitrary, while Eq. (4) holds for the order parameter of the quantum phase transition.
The manuscript is organized as follows: in Sec. II, we provide a detailed derivation of the equations discussed above. In the remaining sections, we draw the finite-temperature phase diagram of the QFI in hallmark systems, recovering the schematic behavior shown in Fig. 1. In Sec. III we study symmetry-breaking QPTs, focusing on the Ising model and the bosonic Josephson junction, while in Sec. IV we consider topological QPTs, in particular the Kitaev chain also with variable range pairing. Finally, discussions and conclusions are reported in Sec. V.

Methods and Results
Quantum Fisher information, multipartite entanglement and quantum coherence. The QFI quantifies the distinguishability between nearby quantum states ρ and ρ φ related by an arbitrary transformation depending on the parameter φ. The Uhlmann fidelity 64 between ρ and ρ φ is ρ ρ provided that p k + p k′ ≠ 0. The QFI has key mathematical properties [65][66][67][68] , that allow the derivation of relevant bounds, see Fig. 2: i) Convexity. The QFI is nonnegative and convex in the state:∑  ii) Additivity. The QFI is additive under tensor product:ˆˆρ iii) Monotonicity. The QFI always decreases under arbitrary parameter-independent completely positive trace-preserving map Λ:ˆρ Q Q with equality for φ-independent unitary transformations.
In the following we will restrict to unitary transformations, ˆˆρ where Ô is a generic Hermitian operator that we will specify below. The unitary transformation only evolves the eigenstates of ρ and leave its eigenvalues unchanged. For unitary transformations, the QFI has the following further properties: iv) The QFI satisfiesˆˆˆˆρ [ , ] 4( ) , with equality for pure states. v) The QFI vanishes if and only if ρ and Ô can be diagonalized simultaneously:ˆˆρ [ , ] 0 (10) Q QFI and quantum coherence. The coherence of a quantum state ρ is defined from its distinguishability with respect to the set of states that are diagonal in a given basis 59 . Here such a basis is given by the eigenstates of the operator Ô , and incoherent states are those satisfying ˆρ = O [ , ] 0. In addition to the properties (i) and (v), the QFI does not increase under operations that conser ve Ô , namely These properties make the QFI a reliable measure of asymmetry 60 , a broad notion of quantum coherence 59,61 . Physically, the concept of asymmetry quantifies how much a state ρ satisfying ˆρ QFI and multipartite entanglement. The key property that makes the QFI a multipartite entanglement witness [40][41][42]45 is the convexity in the state [property (i) above]. We recall that a pure state is κ-partite entangled if it can be written as - , N being the total number of parties in the system) that does not factorize. In other words, κ-partite entanglement indicates the number of parties in the largest nonseparable subset. κ-partite entangled states form a convex set and we can indicate with ˆ- a generic element of the ensemble. As a consequence of Eq. (6), every (pure or mixed) κ-partite entangled state satisfies The maximization is done over all possible κ-separable pure states and we have used A theoretical challenge is to calculate the multipartite entanglement bounds (11) for a given operator Ô , which might be local [40][41][42] or nonlocal 45 . The choice of the operator involved in the calculation of the QFI leads to different entanglement bounds κ b O , . While there is no known systematic method to choose the optimal operator Ô (i.e. the one that allows the detection of the largest class of states), an "educated guess" based on some knowledge of the system allows the corresponding QFI to witness multipartite entanglement close to QPTs for different models. For instance, in models showing symmetry-breaking QPTs, the transition is characterized by the divergence of fluctuations of a local order parameter. We thus expect a large QFI at criticality when Ô is given by the order parameter of the transition 37 . In spin models such as the Ising and the bosonic Josephson junction models this is a collective spin operators (given by the sum of Pauli matrices). In this case we have 41,42 N / is the largest integer smaller or equal than N/κ and r = N − sκ. A QFI larger than this bound witnesses (κ + 1)-partite entanglement between spin-1/2 particles. On the contrary, topological QPTs are not detected by a local order parameter. In order to witness multipartite entanglement in topological models it is thus necessary to calculate the QFI with respect to nonlocal operators. For the one-dimensional short-range Kitaev chain discussed below, an optimal choice of operator is suggested by the correspondence, via the Jordan-Wigner transformation, to the Ising model. Indeed, the QFI is able to detect multipartite entanglement in a topological system 48 when choosing, as operator Ô , the Jordan-Wigner transformation of the local order parameter for the Ising chain (see below). Furthermore, this choice leads 48 to the same multipartite entanglement bounds Quantum Fisher information of thermal states. We consider a generic thermal state at canonical equi- , where Ĥ is the many-body Hamiltonian with eigenenergies E n and corresponding eigenstates |ψ n 〉, T is the temperature, is the partition function. The QFI of ρ T , calculated with respect to the operator Ô , iŝˆ∑ 2 at all temperatures. Equation (12) can be rewritten aŝˆˆ∑ Computing the QFI using Eq. (12) or (13) requires the diagonalization of the full Hamiltonian Ĥ . A calculation using a limited manifold of eigenstates (i.e. in a in a Hilbert space given by the most populated states at temperature T) only leads to approximate results, since the matrix element Ô n m , may couple to energy eigenstates outside the manifold. The calculation of the QFI in a Hilbert subspace (as discussed below for the two-mode approximation) leads to accurate results provided that coupling terms between the subspace and the rest of the Hilbert space induced by the operators Ô are negligible.
The QFI (12) can also be rewritten in the useful form 37ˆ is the imaginary part of the dynamical susceptibility χ O , and ℏω n,m = E n − E m . Using the fluctuation-dissipation relation is the dynamic structure factor, ˆˆ . Equation (15) can thus be rewritten asˆˆˆˆ Q T 2 2 which shows that the QFI can be calculated from the knowledge of the time correlation functions ˆ〈 〉 O t O ( ) . These are known, for instance, in the Ising model for certain operators 70 without requiring the full diagonalization of the Hamiltonian. In the specific example of the Ising model, the calculation of ˆ〈 〉 O t O ( ) is time consuming as it requires the computation of the Pfaffian of a N × N matrix but avoids memory limitations required by full diagonalization. Time-dependent two-spin correlators ˆσ σ 〈 〉 t ( ) x i x j ( ) ( ) are exactly known in the free-fermion representation 71 via the Wick-Bloch-de Dominicis theorem 70 . They can be efficiently computed up to N ≈ 100 exploiting a numerical algorithm 72 . A system of size  N 100 is hard to access due to severe computational cost.
Zero-temperature case. At zero temperature, the QFI becomeŝˆˆ∑ Low-temperature limit and two-mode approximation. Here we demonstrate the inequality (1). Let us consider, for simplicity, a nondegenerate spectrum: the equations that we will obtain in this section can be straightforwardly extended to the degenerate case. At low temperature  Δ T , we can neglect the population of high-energy eigenstates (i.e. taking p n = 0 for n ≥ 2). In this case, using the completeness relation ψ ψ ∑ | 〉〈 | = n n n , Eq. (12) becomesˆˆˆˆˆ∑ [ , ] ( ) e e e e tanh (19) We thus obtain the inequality (1) from which we derive Eqs (2) and (3). The inequality is tight in the limit T → 0, when p 1 = 0. It is also tight at all temperature if and only if ˆ| for all m ≥ 2, i.e. when the operator Ô only couples the ground and the first excited state.
Quantum critical scaling. The scaling behavior in Eq. (4) follows from a standard scaling hypothesis for the dynamical susceptibility 73,74 : where χ is the static susceptibility of the operator Ô with respect to a coupled field, ϕ 0 is a suitable scaling function, ξ is the correlation length and L = N 1/d is the linear system size, being d the system dimension. Inserting Eq. (20) into Eq. (14) we obtain t anh , , . We now take into account that, close enough to the critical with h a suitable scaling function. The behavior of the QFI with respect to relevant quantities can be extracted by setting the dominant rescaling factor b up to which the scale invariance of the system is preserved.
, no significant length scale is induced by temperature. Sufficiently far from criticality, ξ  L and scale invariance is preserved up to b ~ ξ. Equation (22) then implies 37 ˆρ : High-temperature limit. For very large temperature, which predicts a universal 1/T 2 scaling. In the limit T → ∞ we have ρ ∝ T : it commutes with Ô and we find ˆρ Thermalization in a subspace of the full Hamiltonian. All the above equations and the analysis in the following sections implicitly assumes thermal equilibrium in the full Hilbert space. However, it might be possible to have a thermalization only in a Hilbert subspace generated, for instance, by a finite subset of the eigenstates of the full Hamiltonian. This scenario may arise from a metastable equilibrium due to different thermalization time scales of different Hilbert subspaces. In this case, ρ ψ ψ = ∑ | 〉〈 | q T n n n n , where q n ≠ 0 if  ∈ ′ n and q n = 0 otherwise, where ′ is a subspace of the full Hilbert space  with a basis given by the states |ψ n 〉. In this case, the QFI writeŝˆˆ i.e. when the operator Ô does not couple states within the subspace ′. Notice that Ô may couple states in ′ with states outside this subspace, but the latter are not populated before the phase-encoding transformation and do not enter into Eq. (24). In this case, the QFI reduces tôˆ .
In other words, we may have, in this case, that multipartite entanglement increases with temperature and remains large even at high temperature, as recently noticed in a spin-1 system 75 .

Applications: Symmetry Breaking QPTs
In the following we witness multipartite entanglement at finite temperature in two paradigmatic models exhibiting a symmetry-breaking QPT. We first discuss the bosonic Josephson junction (BJJ) model, as it allows for analytical calculations of the QFI at zero as well as finite temperature for large particle numbers. We then focus on the Ising model in a transverse field, which is a common testbed of quantum criticality. The BJJ model can be used to describe a Bose gas in two hyperfine levels coupled by a microwave field 44,76 or in a double-well potential in the tunneling regime [77][78][79][80] , whereas the Ising model has been realized experimentally with ultracold atoms in an optical lattice 81 , trapped ions 82-85 and solid-state platforms 22,23,86 . BJJ model. The BJJ consists of N interacting bosonic particles occupying two weakly-coupled modes 39,87,88 |a〉 and |b〉, e.g. two internal levels of an atom or two wells of an external trapping potential. The system is described by the Hamiltonianˆˆ (2) commutation relations, â and b (ˆ † a and ˆ † b ) are bosonic annihilation (creation) operators for the |a〉 and |b〉 modes, respectively. The coefficient  denotes the characteristic energy scale of the system. The control parameter λ ∈ [0, π/2] rules the interplay between particle-particle interaction, described by Ŝ z 2 , and the linear coupling term, given by Ŝ x . In the thermodynamic limit N → ∞, and for T = 0, Eq. (26) exhibits a QPT at λ c = π/4 between a quantum paramagnetic phase (for λ c < λ ≤ π/2) and a ferromagnetic long-range ordered phase 89,90 (for 0 ≤ λ < λ c ). Equation (25) is a special case of a family of models first introduced by Lipkin, Meshkov and Glick in nuclear physics 89,91 . However, while the Lipkin-Meshkov-Glick model consists of N distinguishable spin-1/2 particles and the full 2 N -dimensional Hilbert space is populated at finite temperature, here we restrict to the N + 1-dimensional Hilbert subspace given by all states symmetric under particle exchange, even at finite temperature.
In the following, we calculate the QFI for a thermal state ρ T and optimize with respect to the operators ˆ= O S x y z , , . The optimization procedure consists in taking the maximum eigenvalue, with k, l = x, y, z 92 . In the large-N limit, we find that, for any λ at T = 0 and for λ π < ≤ atan(1/ 2 ) /2 at T > 0, the optimal operator is the order parameter of the model, ˆ= O S z , while for λ ≤ ≤ 0 a tan(1/ 2 ) at T > 0 it is given by the transverse field, ˆ= O S x . Figure 3 shows the phase diagram of the QFI in the λ-T plane. The diagram has the characteristic V-shape illustrated in Fig. 1. Figure 3 The crossover temperature T cross (λ) (solid white line) can be identified by the inflection points ∂ 2 F Q /∂T 2 = 0 and it follows the energy gap Δ (dashed line) apart a constant multiplication factor, Δ/T cross (λ) ≈ 2.4. Figure 3(b) plots the logarithmic derivative of the QFI with respect to temperature,  1] where N a,b is the number of particles in the mode |a〉 and |b〉, respectively.
is an effective Ginzburg-Landau potential 79 , whose profile has a major role in determining the ground-state structure. Due to the term N −2 in the kinetic energy, the ground state and low-energy excited states are sharply localized around the minima z 0 of the potential V λ (z). Thus, for large N we can Taylor expand the Hamiltonian (28) in the particle spectrum. Equation (29) provides a careful description of the system for energies and temperatures s in (2 cot )/4 n , such that the populated eigenstates are only those strongly localized around z 0 = 0 and negligible at the boundaries z ≈ ±1.
At T = 0, only the ground state of the harmonic oscillator is populated and 39,55 . Notice that this variance diverges in the limit λ → λ c where the potential V λ (z) is no longer harmonic. The QFI is extensive: it linearly grows with the system size N. In particular, at λ = π/2 we have , where ξ 2 is the Wineland spin-squeezing parameter 98,99 . The QFI and the spin-squeezing parameter at zero temperature are shown in Fig. 4(a).
At finite temperature we calculate  , the QFI can be expanded in powers of the small quantity e −Δ/T . For T / 2 is independent on k. Only in the limit T → 0 the two-mode approximation agrees with this expansion. A calculation of Eq. (30) for k → ∞ givesρ which is in very good agreement with the numerical results of ρ F [ ] Q T in the temperature range of interest, as shown in Fig. 4(b) Criticality, λ = λ c . At T = 0 the QFI is superextensive: it scales more rapidly than the system size. A scaling ansatz 37 , with critical exponents Δ Q = 1/3 and z = 1/3, reveals that F Q /N ~ N 1/3 as a function of N 37,39,55 , which is confirmed by numerical calculations. Notice also that the spin squeezing is 93,100 ξ −2 ~ F Q /N. We recall that the energy gap Δ 1 vanishes as N −z . At finite temperature Δ  T 1 , we have F Q /N ~ T −1 according to Eq. (4). At T > 0 it is important to distinguish the case of finite N, where the energy gap Δ 1 ∝ exp(−N|1 − cot λ| 4/3 ) damps exponentially, and the thermodynamic limit N → ∞, where Δ 1 = 0. In the latter case, the ground state is doubly-degenerate (μ = 2) and separated from the doubly-degenerate first excited state (ν = 2) by the energy gap . For arbitrary small but finite temperatures, < Δ  T 0 2 , the system is described by the incoherent mixture ρ ψ ψ ψ ψ = | 〉〈 | + | 〉〈 | ( ) /2     [see Fig. 4(c)], as predicted by Eq. (1) for a purely two-mode approximation. In the thermodynamic limit we have that Δ 1 → 0 and we thus find a discontinuous jump of the QFI from its T = 0 value and the plateau described by Eq. (35). This behavior characterizes the TP of Fig. 1.
The QFI in the different regimes is illustrated in Fig. 4(c), where we show a very good agreement between the analytical predictions and the numerical results. Also in this case, for large enough temperature as the leading term in the Taylor expansion of the tanh function.
Multipartite entanglement. In Fig. 4 in the λ-T plane. Multipartite entanglement witnessed by the QFI is found in the colored region that is bounded by the separability condition ρ It is worth pointing out that multipartite entanglement is considered here among distinguishable spin-1/2 particles restricted to occupy permutationally symmetric quantum states. In practical realizations, such as a Bose-Einstein condensate in double-well trap, these spin-1/2 particles are not addressable. In the limit N → ∞, κ-partite entanglement witnessed by the QFI is found at temperatures for λ > λ c , as obtained from Eq. (31), and at for λ < λ c , following Eq. (35). Equations (36) and (37) are shown as dashed lines in Fig. 4(d). In particular, as noticed above, multipartite entanglement in the ground state of the ferromegnetic phase is extremely fragile to temperature. Moreover, in the thermodynamic limit, we find that no entanglement is witnessed by the QFI at Remarkably, at finite temperature, the QFI detects the same amount of entanglement detected by the spin-squeezing parameter: 1/ξ 2 = F Q /N for T > 0 in the limit  N 1: thermal noise is responsible for a loss of coherence entailing a spread of spin fluctuations in any direction. In particular, λ * is the point at which the spin squeezing ceases to detect entanglement (ξ 2 = 1) even at T = 0 because of the vanishing 〈 〉 S x . When maximizing Eqs (36) and (37) over λ, we obtain that entanglement detected by the QFI survives up to  ≈ .
T / 04, see Fig. 2(d). describes N distinguishable spin-1/2 particles interacting via a nearest-neighbor exchange energy  λ sin (open boundaries are assumed) and subject to a transverse magnetic field of strength  λ cos , with λ ∈ [0, π/2]. The interaction term favors ferromagnetic ordering (with all spins aligned along ±z), while the transverse field favors polarization (with all spins aligned along −x). In the thermodynamic limit N → ∞ and for T = 0, Eq. (39) exhibits a QPT at λ c = π/4 between a paramagnetic phase (for λ c < λ ≤ π/2) and a ferromagnetic phase (for 0 ≤ λ < λ c ). The Ising model in a transverse field is a testbed of quantum criticality 4 .

One-dimensional
Phase diagram. Figure 5 shows the phase diagram of the QFI in the λ-T plane, where the QFI is optimized with respect to the collective operator ˆσ = ∑ = O i N x y z 1 2 1 , , . In particular, the black line in Fig. 5(a)  1 (on the right side of the line). The diagram displays the characteristic V-shaped structure radiating from the critical point, as in Fig. 1. In Fig. 5(a) we can recognize the CP (for λ > λ c ) and the TP (for λ < λ c ). Therein, the QFI ρ F [ ] Q T is approximatively constant as a function of temperature and equal to its low-temperature value ρ F [ ] Q 0 -we recall that ρ 0 is given by the ground state |ψ 0 〉 in the CP and by the incoherent superposition of the two lowest energy eigenstates in the TP. We also see that T cross (solid white line) follows the energy gap Δ (dashed line). The finite jump discontinuity of T cross that is visible in the figure is due to the sudden change of optimal operator Ô and prominently manifests only for small N. In Fig. 5 While the spin squeezing agrees with the QFI close to λ = π/2, it sharply decays at λ c . Indeed, a numerical study as a function of N (up to N = 500) for λ = λ c reveals that ξ 2 = 1 and in particular it does not scale with N, see also refs 47,100 . This is in sharp contrast to the results of the BJJ model where the QFI and the spin-squeezing parameters for the ground state have the same scaling at the critical point, see Fig. 4(a). The typical behavior of ρ F [ ] Q T as a function of temperature for λ > λ c is shown in Fig. 6(b). The solid line is , namely Eq. (1) with μ = ν = 1. For λ < λ c , the ground state becomes doubly degenerate (μ = 2) in the thermodynamic limit N → ∞, and so does the first excited state. For finite N the gap Δ 1 is exponentially small: a finite-size analysis extended to = ÷ N 10 10 3 reveals that the energy gap vanishes expo- T as a function of temperature is shown in Fig. 6(c): the decay from the zero-temperature value occurs around a finite in the vicinity of the critical point. In both panels, the white solid line is the crossover temperature T cross (λ). The blue and red dashed lines indicate Δ 1 and Δ 2 , respectively. In both panels, N = 50. (exponentially small in N) temperature  Δ T 1 , while a slower decay takes place for Multipartite entanglement. The multipartite entanglement between spin-1/2 particles detected by the QFI survives at finite temperature in the colored region of Fig. 6(d), bounded by ρ . This region fans out from the zero-temperature noninteracting λ = π/2 corner, where the ground state of Eq. (39) is separable and . We can compare the condition ρ with the spin-squeezing coefficient ξ 2 = 1 (dashed line) 100 . The loss of spin squeezing follows the loss of thermal entanglement only for λ ≈ π/2, while around λ c we have entangled states recognized by the QFI that are not spin squeezed. Furthermore, the multipartite entanglement region in Fig. 6(d) reaches a maximum extension  ≈ .
T / 06 at λ = 0.4π. Notice that this threshold temperature is higher than the one for the BJJ model. Interestingly, the threshold  ≈ .
T / 03 at λ = λ c is consistent with the temperature up to which other thermal signatures of criticality persists 6,23 . Finally, for  λ π . 0 2 multipartite entanglement is no more witnessed by the QFI at any T > 0.

Applications: Topological QPTs
In the following we study the one-dimensional Kitaev model 103,104 for spinless fermions hopping in a tight-binding lattice with p-wave superconducting pairing. With respect to the original model 103 , we consider variable range for the pairing 105,106 . The Hamiltonian iŝˆˆˆ † H J a a n d a a 2 ( Hc ) , that relates fermionic creation and annihilation operators to the Pauli ladder operators σ + and σ − of spin-1/2 particles 71 . Via Jordan-Wigner trasformation, the nearest-neighbor Kitaev chain (α = ∞) maps into the XY model in a transverse field 71,107 permits to recover Eq. (39) for the fully-anisotropic case J = Ω. Within the Jordan-Wigner transformation, the operators in Eqs. (42) and (43) become the local collective spin operators . By mean of this transformation, each lattice site maps into an effective spin-1/2 particle: the z component of the spin is local in the site (the empty lattice corresponding to spin-down, the filled lattice to spin-up), the other x and y components are nonlocal. We can then use the bound discussed above [40][41][42] to witness κ-particle entanglement between the N effective spin-1/2. Specifically, In the following, we set equal pairing and hopping strengths Ω = J and take  λ = J 2 cos and  μ λ = 2 sin in order to describe the whole phase diagram through a bounded control parameter λ ∈ [−π/2, π/2]. The Hamiltonian (40) thus rewrites aŝˆˆˆ † Kitaev model with short-range pairing. We consider the case α = ∞ where pairing occurs only within nearest-neighbor lattice sites. In this case f α (k) = 2 sin k. As shown by Eq. (41), the energy gap between the ground state and the first excited state vanishes as Δ ~ N −1 in the thermodynamic limit at λ c = ±π/4 (for k = π and k = 0, respectively). These quantum critical points separate a different nontrivial phase with W = 1 (for |λ| < π/4) from a trivial phase with W = 0 (for |λ| > π/4). This behavior is common for short-range pairing 48,105 , α > 1.
Phase diagram. As expected, the results of our study are very similar to the case of the quantum Ising model discussed in Sec. III. There is a major difference though: in the Kitaev model the energy gap Δ = E 1 − E 0 remains finite for every λ ≠ λ c , i.e. away from the critical points. Therefore, the system does not host a gapless phase, differently from the ferromagnetic phase of the Ising model. This is a direct consequence of the fact that the Kitaev model is studied here in the closed chain. In the open chain, the Kitaev model hosts a gapless phase for |λ| < π/4, related to the presence of Majorana edge modes. In Fig. 7 for any λ and T. We recognize the presence of plateaus at low temperature and the characteristic V-structure around the critical points. Only QPs are present, due to the nondegenerate nature of the ground state. The phase diagram is invariant under change of sign of the chemical potential λ → −λ, as expected from the particle-hole symmetry of the Hamiltonian 104 . The crossover temperature T cross (λ) (solid white line) follows the energy gap (dashed line) for |λ| > π/4, with Δ/T cross ≈ 2.7. In the region |λ| < π/4, T cross is instead smoothed, due to the quasi-degeneracy of the excited states. In Fig. 7(b) we plot the logarithmic derivative Multipartite entanglement. Figure 8 illustrates the multipartite entanglement witnessed by the QFI. Panel (a) shows the QFI of the ground state, F Q [|ψ 0 〉]/N, as a function of λ. The trivial phase |λ| > π/4 is characterized by an extensive scaling of the QFI for increasing system size N. At λ = +π/2 (λ = −π/2), we find F Q [|ψ 0 〉] = N, according to the fact that the ground state is a separable state of occupied (empty) sites ψ , where {|n〉 i } is the occupation basis and n ∈ {0, 1} is the occupation number at the i-th site. Divergence of multipartiteness F Q /N ~ N is instead observed in the phase with nonzero winding number (|λ| < π/4) 48 . In particular, F Q /N = N at λ = 0. The QPT at λ c is signalled by a sudden change in the scaling F Q /N ~ N 3/4 , that is associated to the specific algebraic asymptotic decay observed for the two-site correlation functions 48 . Figure 8(b) shows the witnessed multipartite entanglement at finite temperature. Within the region |λ| < π/4, superextensive multipartite entanglement in the ground state survives at finite temperature. This robustness is due to the nondegenerate nature of the ground state for any λ ≠ λ c and it is in sharp contrast with the ferromagnetic phase of the BJJ and Ising model, where a superextensive QFI decays exponentially with N at finite temperature. In particular, at λ = 0, where the first excited state is N-fold degenerate, Eq. (1) predicts for low temperaturê  that shows a plateau up to temperatures  ≈ T N / 2/ log . The thermal decay of the QFI is at most logarithmic in N. This behavior is confirmed in Fig. 8(c) where, for a fixed temperature, we plot F Q /N as a function of N. We see that the Heisenberg scaling F Q /N ~ N survives at finite temperature up to   N e T 2 / . For larger system size, temperature is responsible for a softening of the power-law scaling. It is worth noticing that this effect is not related to a vanishing gap, as in the ferromagnetic phase of the Ising model, but it is rather due to the diverging degeneracy of the first excited state. As emphasized in Eq. (1), the robustness of the QFI to temperature depends indeed on this degeneracy.
In Figs 9 and 10 we plot the QFI phase diagram and the witnessed multipartite entanglement, respectively. The operator that maximizes the QFI of the ground state is found to be ˆ+ O x ( ) for λ ≤ λ c and ˆ− O y ( ) for λ ≥ λ c , see Fig. 10(a). The two phases at λ < λ c and λ > λ c are characterized by a diverging multipartiteness F Q /N ~ N 3/4 , while F Q /N ~ N 1/2 at criticality λ = λ c 48 . These scaling behaviors survive at low temperature as shown in Fig. 9(a) where we plot ρ ψ | 〉   characterizes the CP around λ c . In Fig. 9(b), the green region highlights values β ∈ [−0.53, −0.47]. Figure 10(b) highlights the region of the λ-T phase diagram where the QFI witnesses multipartite entanglement (colored region). In particular, in Fig. 10(c) we plot ρ F N [ ]/ Q T as a function of N for λ = 0 and different temperatures. For sufficiently small temperature the QFI is bounded by Eq. (1). In this case, the evaluation of the degeneracy of the first excited state ν is not easily practicable: an analysis of Eq. (41) shows that the number of states in a small interval centered around the energy of the first excited state increases with N. We thus superpose in Fig. 10
In this manuscript we have discussed the universal behavior of the QFI for systems at thermal equilibrium close to a QPT. At low-temperature, the QFI is lower bounded by a simple function that only depends on the structure of the two low-lying energy levels and is factorable in a finite-temperature and a zero-temperature contributions. This feature allows to draw a V-shaped phase diagram for the QFI centered at the critical point, Fig. 1, which is common to both symmetry-breaking and topological QPTs. We showed the existence of a universal low-temperature region -the CP -where thermal decay of the QFI is ruled by few fundamental critical exponents. This region fans out from the critical point and can be identified as a quantum critical regime where quantum coherence has a behavior controlled by the transition and competes with thermal fluctuations. The universal behavior is lost at surprisingly high temperatures.
Finally, the analysis has emphasized the robustness of multipartite entanglement at finite temperature. In particular, a superextensive QFI (with a scaling at the Heisenberg limit F Q ~ N 2 ) survives up to high temperatures, T ∝ 1/log N in topological systems with large finite size. This is an important difference with respect to models showing symmetry-breaking QPTs. In the latter systems multipartite entanglement is generally found at finite temperature in the disordered phase and the superextensive QFI that characterizes the ground state of the ordered phase is exponentially fragile against temperature, being lost for T ∝ e −N . Note added in Proofs: Short before the submission of this manuscript, we became aware of the similar work 108 by I. Frerot and T. Roscilde. There, the quantum variance, a quantity related to the quantum Fisher information, is studied at finite temperature around the critical point of many-body quantum models and used to characterize a quantum critical regime.