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.

Introduction

A quantum information approach to the study of quantum phase transitions (QPTs)[1-3] sheds new light on these many-body phenomena[4] and pushes our understanding of the puzzling behavior of strongly-correlated systems[5-7] beyond standard methods in statistical mechanics[8]. Entanglement in the ground state of a many-body Hamiltonian – where and are non-commuting operators and λ is a control parameter – has been extensively investigated close to a quantum critical point λc[1-3,9-13]. Yet, less is known about the survival of entanglement at finite temperature[9,14], including the peculiar quantum critical region that fans out from λc[4,15,16]. This regime is particularly interesting due to the competition of thermal and quantum fluctuations[4,15,16] and plays a key role in interpreting a wide variety of experiments in synthetic matter[17-26].

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-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-43] and it is currently attracting considerable interest[37,44-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-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-58]. The QFI plays a central role in the theory of quantum coherence[59-63]: it quantifies the coherent extent of a generic state over the eigenstates of the operator , vanishing if and only if . Multipartite entanglement is witnessed when the QFI overcomes certain finite bounds: as discussed below[41,42], 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. At low temperature, the QFI satisfies the inequality

Figure 1

Schematic general behavior of the scaling of the QFI in the vicinity of a critical point. 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/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 Tcross (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, β = −Δ/z, according to Eq. (4). For temperatures larger than Tmax (dotted line) – approximatively equal to the maximum energy of the spectrum – the QFI enters the MEP where β = −2. In the crossover grey regions the thermal decay is non-universal.

Here, 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 Egs and first excited state of energy Eex, respectively, and Δ = Eex − Egs is the first energy gap in the many-body spectrum. Equation (1) is valid for and shows that, regardless on the microscopical details of the system, the lower bound to 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 T → 0, where 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) givesand shows that the QFI is bounded from below by a constant for , where Tcross ≈ Δ/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 . 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 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 , where Tcross ≈ Δ/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 . In other words, the QFI of the ground state 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 that remains constant up to Tcross. 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 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), Δ is the exponent[37] that characterizes the finite-size scaling of the QFI with respect to N at T = 0 and λ = λc, i.e. , 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 , 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 .

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–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 , where is the QFI. In terms of the spectral decomposition (with p ≥ 0 and ) we have provided that p + p ≠ 0. The QFI has key mathematical properties[65-68], that allow the derivation of relevant bounds, see Fig. 2:

Figure 2

Bounds of the QFI. For unitary phase-encoding transformations, if and only if the state is incoherent. Among quantum coherent states, , we can find bounds to the QFI depending on the entanglement properties of the state: for separable states [orange region, where b1 is also indicated as shot-noise (SN) limit], for κ-partite entangled states with 1 ≤ κ ≤ N − 1 [green region], and for all possible states [where b is also indicated as Heisenberg limit (HL)]. The bounds b depend, in general, on the operator .

Convexity. The QFI is nonnegative and convex in the state:for any state and q ≥ 0.

Additivity. The QFI is additive under tensor product:

Monotonicity. The QFI always decreases under arbitrary parameter-independent completely positive trace-preserving map Λ: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:

The QFI satisfieswith equality for pure states.

The QFI vanishes if and only if and can be diagonalized simultaneously:

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 . In addition to the properties (i) and (v), the QFI does not increase under operations that conserve , namely for maps Λ satisfying . 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 changes when applying the unitary transformation . The changes in the state can be used to estimate the phase ϕ with a nonvanishing sensitivity[39,69] in a sensor implementing the transformation .

QFI and multipartite entanglement

The key property that makes the QFI a multipartite entanglement witness[40-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 [38], where |ψ〉 is a state of N ≤ κ parties (with , 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 , where

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-42] or nonlocal[45]. The choice of the operator involved in the calculation of the QFI leads to different entanglement bounds . 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] , where 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 equilibrium, , where is the many-body Hamiltonian with eigenenergies E and corresponding eigenstates |ψ〉, T is the temperature, and is the partition function. The QFI of , calculated with respect to the operator , iswhere . Notice that at all temperatures. Equation (12) can be rewritten as

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 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]where is the imaginary part of the dynamical susceptibility χ, and ℏω = E − E. Using the fluctuation-dissipation relation we can writewhere  =  is the dynamic structure factor, , and we have used the property . Equation (15) can thus be rewritten aswhich shows that the QFI can be calculated from the knowledge of the time correlation functions . 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 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 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 is hard to access due to severe computational cost.

Zero-temperature case

At zero temperature, the QFI becomesin case the ground state has a degeneracy μ (we have indicated as the degenerate eigenstates, with d = 1, … μ), such that , and reduces to in absence of degeneracy (μ = 1).

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 , we can neglect the population of high-energy eigenstates (i.e. taking p = 0 for n ≥ 2). In this case, using the completeness relation , Eq. (12) becomes

Notice that the second and third terms in Eq. (18) are always positive (at all temperatures), which impliesWe thus obtain the inequality (1) from which we derive Eqs (2) and (3). The inequality is tight in the limit T → 0, when p1 = 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 = N1/ is the linear system size, being d the system dimension. Inserting Eq. (20) into Eq. (14) we obtainwhere . We now take into account that, close enough to the critical point, Δ ~ δ, ξ ~ δ− and χ ~ δ−, where δ = |λ − λc| and z, ν and γ are critical exponents. Under coarse-graining transformation on the system, lengths scale as l → l′ = b−1l (b > 1), while N → N′ = b−N. A dimensional analysis reveals that Δ → Δ′ = bΔ and δ → δ′ = b1/δ, whereas the scaling function g only depends on adimensional variables and does not scale under length rescaling. Thus, under coarse graining the QFI transforms according towith 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.

At small temperatures , no significant length scale is induced by temperature. Sufficiently far from criticality, and scale invariance is preserved up to b ~ ξ. Equation (22) then implies[37] for . Conversely, at the critical point , the constituents of the system are correlated on a scale : the system experiences finite-size effects and it remains scale invariant up to b ~ L. Equation (22) gives[37] the scaling of the QFI with N for : , where we have used the Fisher relation γ/ν = 2 − η and defined Δ = 2 − η − z.

On the contrary, thermal fluctuations dictate a dominant length scale if . Scale invariance is expected to be broken at the scale b ~ T−1/. Thus, Eq. (22) provides the scaling of the QFI with temperature valid for : , namely Eq. (4).

High-temperature limit

For very large temperature, , we can expand . Equation (12) becomeswhich predicts a universal 1/T2 scaling. In the limit T → ∞ we have : 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, , where q ≠ 0 if and q = 0 otherwise, where is a subspace of the full Hilbert space with a basis given by the states |ψ〉. In this case, the QFI writes

Interestingly, the second term in Eq. (24) can vanishes. This occurs when , 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 to

This equation predicts results that are completely different from the ones discussed above. For instance, if the excited states are characterized by values of that are larger than those for the low-lying states, then may increase with temperature. Furthermore, in the large-T limit, taking , we obtain that the QFI saturates a finite constant value, . 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-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 Hamiltonianwhere , and satisfy SU(2) commutation relations, and ( and ) 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 , and the linear coupling term, given by . 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-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 and optimize with respect to the operators . The optimization procedure consists in taking the maximum eigenvalue, , of the matrixwith k, l = x, y, z[92]. In the large-N limit, we find that, for any λ at T = 0 and for at T > 0, the optimal operator is the order parameter of the model, , while for at T > 0 it is given by the transverse field, .

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(a) plots . The crossover temperature Tcross(λ) (solid white line) can be identified by the inflection points ∂2F/∂T2 = 0 and it follows the energy gap Δ (dashed line) apart a constant multiplication factor, Δ/Tcross(λ) ≈ 2.4. Figure 3(b) plots the logarithmic derivative of the QFI with respect to temperature, , giving the scaling exponent for the thermal decay of . We clearly distinguish regions characterized by constant values of β and corresponding to the different plateaus of Fig. 1: β = 0 in the TP and QP, β = −1 in the CT and, finally, β = −2 in the MEP.

Figure 3

Phase diagram of the BJJ model. (a) QFI normalized to its low-temperature value, (color scale) as a function of λ and T. The region where the low-temperature behavior survives is highlighted by the orange color. The black dotted line at separates the regions where the optimal parameter is (on the left) and (on the right). (b) Scaling coefficient (color scale) as a function of λ and T. The dotted line is the upper bound of the spectrum, . The inset shows β as a function of T at λ/π = 0.4: the different plateaus are clearly visible. In both panels N = 2000, the solid white curve is the crossover temperature Tcross(λ) following the energy gaps Δ1 (dashed blue line) and Δ2 (dashed red line).

These results can be fully understood analytically in the large-N limit via an Holstein-Primakoff approach[93,94]. An expansion in powers of 1/N allows to rewrite Eq. (26) as[95-97]

Here, z = (N − N)/N ∈ [−1, 1] where N 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 z0 of the potential V(z). Thus, for large N we can Taylor expand the Hamiltonian (28) around z0 and retain only the quadratic terms in z − z0.

Paramagnetic phase, λ > λc

In this case, Eq. (28) reduces to the Hamiltonian of an effective harmonic oscillator centered at z0 = 0,with gap in the particle spectrum. Equation (29) provides a careful description of the system for energies and temperatures , such that the populated eigenstates are only those strongly localized around z0 = 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 F[|ψ0〉] = N, consistently with the fact that the ground state is separable and given by the coherent spin-polarized state . Furthermore, we have squeezing of the spin fluctuation below the projection-noise limit and sufficiently high such that the state is also spin squeezed [39,55], 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).

Figure 4

QFI for the BJJ model. (a) Fisher density F[|ψ0〉]/N (blue line) and inverse spin-squeezing parameter 1/ξ2 (orange line) as a function of λ for the ground state of Eq. (26). The two lines superpose for . The vertical dashed line signals the critical point λc = π/4. Panels (b and c) show the Fisher density (dots) as a function of T for (b) λ = 0.3π > λc and (c) λ = 0.24π < λc. Solid lines are analytical curves, Eqs (30) and (34), for different values of the cutoff k. The vertical dashed lines indicates T = Δ1,2. In panels (a–c) the shaded area indicates multipartite entanglement F/N > 1. (d) Fisher density (color scale) in the λ-T phase diagram. Multipartite entanglement is witnessed at nonzero temperature in the colored region, where , while corresponds to the white region. The dashed line is the analytical boundary of in the thermodynamic limit, given by Eq. (36) for λ > λc and Eq. (37) for λ < λc. In all the panels, numerical data are obtained for N = 2000.

At finite temperature we calculate taking into account the first k harmonic oscillator modes and using with :

For k = 2 we recover the two-mode approximation leading to Eq. (1), which agrees with the case k > 2 when . In fact, for low temperatures , the QFI can be expanded in powers of the small quantity e−Δ/. For k = 2 we obtain . For k > 2, the result 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 → ∞ giveswhich is in very good agreement with the numerical results of in the temperature range of interest, as shown in Fig. 4(b). The Fisher matrix Eq. (27) reads . Thus, the optimal operator is for all T. For we find at any value of λ. However, Eq. (31) is valid for all T only in the thermodynamic limit: for a finite system size, at we recover the MEP regime, where .

Criticality, λ = λc

At T = 0 the QFI is superextensive: it scales more rapidly than the system size. A scaling ansatz[37], with critical exponents Δ = 1/3 and z = 1/3, reveals that F/N ~ N1/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/N. We recall that the energy gap Δ1 vanishes as N−. At finite temperature , we have F/N ~ T−1 according to Eq. (4).

Ferromagnetic phase, λ < λc

For sufficiently large N, we can calculate the QFI using the semiclassical model of Eq. (28). The effective potential V(z) has a symmetric double-well profile[96] with two minima located at . Equation (28) takes the formwhen locally approximating each well as a harmonic oscillator.

At T = 0 the QFI, calculated for the ground state of Eq. (32), is superextensive[39,96], , see Fig. 4(a). In particular, in the limit λ → 0, the ground state is the NOON state , given by a coherent symmetric superposition of N particles in mode |a〉 and |b〉, that has the highest possible value of the QFI[40], F[|ψ0〉]/N = N.

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, , the system is described by the incoherent mixture . Its QFI iswhere arises from the optimization of the operator: for λ ≤ λ, while for λ > λ. We can calculate the QFI using a k-mode approximation (i.e. taking the first k states in each harmonic well). By means of , where and , we can evaluate the QFI in Eq. (12), obtainingthat becomeswhen taking the limit k → ∞. It should be noticed that Eq. (33) is a factor N smaller than the zero-temperature value F[|ψ0〉]. For large but finite N, the rapid decay from F[|ψ0〉] to is described by the decaying function [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 we recover as the leading term in the Taylor expansion of the tanh function.

Multipartite entanglement

In Fig. 4(d) we plot 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 temperaturesfor λ > λc, as obtained from Eq. (31), and atfor λ < λ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 T > 0 for λ ≤ λ*, whereRemarkably, at finite temperature, the QFI detects the same amount of entanglement detected by the spin-squeezing parameter: 1/ξ2 = F/N for T > 0 in the limit : 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 . When maximizing Eqs (36) and (37) over λ, we obtain that entanglement detected by the QFI survives up to , see Fig. 2(d).

One-dimensional Ising model in a transverse field

The one-dimensional quantum Ising chain in a transverse field[101,102],describes N distinguishable spin-1/2 particles interacting via a nearest-neighbor exchange energy (open boundaries are assumed) and subject to a transverse magnetic field of strength , 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].

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 . In particular, the black line in Fig. 5(a) marks a region where the optimal operator is (on the left side of the line) and the one where the optimal operator is the order parameter of the transition, (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 is approximatively constant as a function of temperature and equal to its low-temperature value – we recall that 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 Tcross (solid white line) follows the energy gap Δ (dashed line). The finite jump discontinuity of Tcross 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(b) we plot the logarithmic derivative in the vicinity of the critical point, which provides the scaling of the QFI with temperature. According to the scaling ansatz (see Sec. II), using Δ = 3/4 and z = 1[37], we find in the QP. Numerical results are plagued by finite size effects and we do not observe a clear plateau for β. We argue (supported by a finite-size study, yet limited to ) that the CP, where β = −0.75, approximatively coincides with the green region in the figure, which highlights values of β ∈ [−0.78, −0.72].

Figure 5

Phase diagram of the Ising model in transverse field. (a) QFI normalized to its low-temperature value, (color scale), in the λ-T phase diagram. (b) Scaling coefficient in the vicinity of the critical point. In both panels, the white solid line is the crossover temperature Tcross(λ). The blue and red dashed lines indicate Δ1 and Δ2, respectively. In both panels, N = 50.

The behavior of the QFI is further inspected in Fig. 6. In panel (a) we show the QFI (blue line) for the ground state of Eq. (39) and compare it to the spin squeezing (orange line). Similarly to the BJJ model, for λ > λc the QFI is larger than N – signaling multipartite entanglement – but extensive, i.e. the Fisher density F[|ψ0〉]/N does not scale with N[47]. For λ ≤ λc, the QFI is superextensive. It scales as F[|ψ0〉]/N ~ N for λ < λc[47] and as F[|ψ0〉]/N ~ N3/4 at λ = λc[37]. 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 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 reveals that the energy gap vanishes exponentially for λ < λc, with nΔ = a|1 − cot λ| (a ≈ 1.5, b ≈ −0.85). The behavior of as a function of temperature is shown in Fig. 6(c): the decay from the zero-temperature value occurs around a finite (exponentially small in N) temperature , while a slower decay takes place for . The constant plateau found for defines the TP.

Figure 6

QFI for the Ising model in transverse field. (a) Fisher density F[|ψ0〉]/N (blue line) and inverse spin squeezing (orange line) for the ground state of Eq. (39) as a function of λ. The vertical dashed line signals the critical point λc. Panels (b and c) show the typical decay of the Fisher density as a function of T in the paramagnetic (b) and ferromagnetic (c) phase. The solid lines are or . In panels (a–c) the shaded area indicates multipartite entanglement. (d) Fisher density (color scale) in the λ-T phase diagram. The dashed line is the spin-squeezing boundary ξ2 = 1. In all panels N = 50.

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 at λ = 0.4π. Notice that this threshold temperature is higher than the one for the BJJ model. Interestingly, the threshold at λ = λc is consistent with the temperature up to which other thermal signatures of criticality persists[6,23]. Finally, for 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 iswhere is a fermionic creation operator at i-th site (satisfying the anticommutation relation ) and counts the number of fermions in the i-th site. The amplitude of the hopping between different lattice sites is J and the chemical potential of the chain is μ. The superconducting pairing has strength Ω and range specified by , where is a site-to-site distance and α > 0: α → ∞ corresponds to nearest-neighbor pairing, while α = 0 accounts for infinite-range pairing. For a closed ring, () if (). In Eq. (40) we consider antiperiodic boundary conditions ().

The Hamiltonian (40) can be diagonalized exactly by a Bogoliubov transformation[71] for any α. The quasiparticle spectrum reads[105]where are the quasimomenta of the excitations (n = 0, 1, …, N − 1) and . The first energy gap in the many-body spectrum corresponds to the energy necessary to create one elementary excitation. The ground state of the Kitaev chain Eq. (40) reads , with , , being the Fourier transform of and |0〉 denoting the vacuum of quasiparticles. The ground state hosts different topological phases that can be characterized by the winding number . At zero temperature, the mean number of fermions in the system is given by .

Here, we study the QFI calculated with respect to the nonlocal operators (ρ = x, y)and the local operators

This choice is suggested by the Jordan-Wigner transformation and , 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] . A local collective rotation 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 , being , and . 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-42] to witness κ-particle entanglement between the N effective spin-1/2. Specifically, signals (κ + 1)-partite entanglement. Recently, operators have been used to demonstrate the superextensivity of the QFI at zero temperature in the different phases of the Kitaev model (40), for both short-range and long-range pairing[48].

In the following, we set equal pairing and hopping strengths Ω = J and take and in order to describe the whole phase diagram through a bounded control parameter λ ∈ [−π/2, π/2]. The Hamiltonian (40) thus rewrites as

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.

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 Δ = E1 − E0 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(a) we plot the λ-T phase diagram for . The optimal operator maximizing the QFI is found to be 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 Tcross(λ) (solid white line) follows the energy gap (dashed line) for |λ| > π/4, with Δ/Tcross ≈ 2.7. In the region |λ| < π/4, Tcross is instead smoothed, due to the quasi-degeneracy of the excited states. In Fig. 7(b) we plot the logarithmic derivative around the critical point λ = π/4. The QPT is characterized by the same critical exponents as the Ising model, Δ = 3/4 and z = 1, and we thus expect a thermal decay , according to Eq. (4). The region where β ∈ [−0.8, −0.7] is highlighted in the figure.

Figure 7

Phase diagram of the Kitaev chain with short-range pairing. (a) QFI normalized to its low-temperature value, (color scale), in the λ-T phase diagram. (b) Scaling coefficient as a function of λ and T. In both panels the white line is Tcross and the blue dashed line is the energy gap Δ. Here N = 50 and α = 100.

Figure 8 illustrates the multipartite entanglement witnessed by the QFI. Panel (a) shows the QFI of the ground state, F[|ψ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[|ψ0〉] = N, according to the fact that the ground state is a separable state of occupied (empty) sites (), where {|n〉} is the occupation basis and n ∈ {0, 1} is the occupation number at the i-th site. Divergence of multipartiteness F/N ~ N is instead observed in the phase with nonzero winding number (|λ| < π/4)[48]. In particular, F/N = N at λ = 0. The QPT at λc is signalled by a sudden change in the scaling F/N ~ N3/4, that is associated to the specific algebraic asymptotic decay observed for the two-site correlation functions[48].

Figure 8

QFI for the Kitaev chain with short-range pairing. (a) Fisher density F[|ψ0〉]/N as a function of λ for the ground state of Eq. (44) with N = 50 and α = 100. The vertical dashed lines signal the critical points λc. The shaded area marks entanglement, F[|ψ0〉] > N. (b) Fisher density (color scale) in the λ-T plane for N = 50. The colored area corresponds to . (c) Scaling of as a function of N for different temperatures. The thick black line is the Heisenberg limit F = N2, the dashed lines are the bound in Eq. (45).

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 temperaturewhere we have used F[|ψ0〉] = N2, , μ = 1 and ν = N. For large N the right-hand side of Eq. (45) can be well approximated by that shows a plateau up to temperatures . 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/N as a function of N. We see that the Heisenberg scaling F/N ~ N survives at finite temperature up to . 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.

Kitaev model with long-range pairing

We study the Kiteav model with α = 0 where pairing involves fermions in arbitrarily-distant sites. In this case f(k) = cot(k/2), which diverges at k = 0. The energy gap vanishes as Δ ~ N−1 at λc = π/4 (for k = π). The winding number is[48,105] W = +1/2 for λ < λc, and W = −1/2 for λ > λc. The symmetry under λ → −λ is lost, due to the loss of particle-hole symmetry.

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 for λ ≤ λc and for λ ≥ λc, see Fig. 10(a). The two phases at λ < λc and λ > λc are characterized by a diverging multipartiteness F/N ~ N3/4, while F/N ~ N1/2 at criticality λ = λc[48]. These scaling behaviors survive at low temperature as shown in Fig. 9(a) where we plot in the λ-T phase diagram. For low temperatures , we recognize two QPs at both sides of the critical point. Since the transition is characterized by Δ = 1/2 and z = 1, the thermal decay for 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 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(c) the numerical data (dots) for to the curve (dashed lines) as suggested by Eq. (1), where c = 1.4 is a fitting parameter and (in the thermodynamic limit). For , the QFI grows as .

Figure 9

Phase diagram of the Kitaev chain with long-range pairing. (a) QFI normalized to its low-temperature value, (color scale), in the λ-T phase diagram. (b) Scaling coefficient . In both panels the white line is Tcross and the blue dashed line is the energy gap Δ. Here N = 50 and α = 0.

Figure 10

QFI for the Kitaev chain with long-range pairing. (a) Fisher density F[|ψ0〉]/N as a function of λ for the ground state of Eq. (44) with N = 50 and α = 0. The QFI is calculated using both the operators (dark blue line) and (light blue line). The vertical dashed line signals the critical point λc, while the shaded area marks multipartite entanglement. (b) Fisher density (color scale) in the λ-T plane for N = 50. (c) Scaling of for increasing N. The dashed lines are Eq. (45) for different values of T.

Conclusions

The QFI, as a multipartite entanglement witness, allows to study strongly-correlated systems from a quantum information perspective and is thus attracting increasing interest[37,44-50,56]. Differently from bipartite/pairwise entanglement measures the QFI[37] (or close lower bounds[44,53,54]) can be extracted experimentally in arbitrary large systems of atomic ensembles and solid-state platforms.

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 ~ N2) 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−. 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.