Characterization of an operational quantum resource in a critical many-body system

Quantum many-body systems have been extensively studied from the perspective of quantum technology, and conversely, critical phenomena in such systems have been characterized by operationally relevant resources like entanglement. In this paper, we investigate robustness of magic (RoM), the resource in magic state injection based quantum computation schemes, in the context of the transverse field anisotropic XY model. We show that the the factorizable ground state in the symmetry broken configuration is composed of an enormous number of highly magical H states. We find the existence of a point very near the quantum critical point where magic contained explicitly in the correlation between two distant qubits attains a sharp maxima. Unlike bipartite entanglement, this persists over very long distances, capturing the presence of long range correlation near the phase transition. We derive scaling laws and extract corresponding exponents around criticality. Finally, we study the effect of temperature on two-qubit RoM and show that it reveals a crossover between dominance of quantum and thermal fluctuations.


Introduction
Many-body systems provide a rich playground for different phases of matter [1,2], mediating several tasks including quantum communication [3,4], quantum metrology [5][6][7][8], remote gate implementation [9,10], and refrigeration [11,12]. Ground states of many-body systems are generically highly correlated at criticality [13,14], and usually offer significant amounts of quantum resources, such as entanglement [2,[15][16][17] or quantum Fisher information [5,6,18]. Most of second order quantum phase transitions are accompanied by a spontaneous symmetry breaking between two (or more) degenerate ground states at the critical point [1], with the exception of exotic impurity quantum phase transitions [19][20][21]. In practice, the breaking of symmetry is due to unavoidable environmental fluctuations. The symmetry broken ground states may show interesting properties, including the existence of points on the phase diagram of quantum many-body systems, for which the ground state is factorizable [22][23][24] into uncorrelated pure states. These points may be reached from critical points via local operations and classical communication [25]. So far, the factorizable points are considered worth avoiding due to absence of entanglement. It would be interesting to explore whether these factorizable ground states can be useful in quantum technologies not directly relying on entanglement, which is the most celebrated resource in quantum technologies [26,27]. Many-body systems have been extensively studied for their entanglement content [2], in particular, near quantum phase transitions. In fact, total quantum correlation content for the system is expected to be maximal at the critical point [14]. However, harnessing such entanglement is practically challenging as it demands complex operations over several particles. One practically useful case is the entanglement between individual particles. Remarkably, in the transverse anisotropic XY model which is being investigated in our work, bipartite entanglement in the (thermal) ground state does not persist beyond the second nearest neighbour [14]. One may ask whether many-body systems can provide other quantum resources which survive between individual particles, over longer distances. Whether such persistence can be used to detect the quantum critical region for a symmetry unbroken thermal state in a manner analogous to nearest-neighbour entanglement [28], is also worth investigating.
In this paper, we investigate the phase diagram of the transverse field anisotropic XY spin chain for its single-and two-qubit magic contents. Two regimes will be considered: zero temperature in which one of the symmetry broken ground states naturally emerge as the quantum state of the system, and finite temperature regime in which degenerate symmetry unbroken ground states can equally contribute in the thermal equilibrium state. For the ground state, we observe that the single-qubit magic lies within the stabilizer polytope throughout the disordered phase, and becomes magical immediately after the critical point. We demonstrate the existence of two different scaling behaviours in the vicinity of this point, and extract finite size scaling coefficients. This model allows in-principle exact mining of very large number of perfect unencoded magic states at its factorized ground state (FGS) and is robust against imperfect tuning of the Hamiltonian. The correlation of magic content at criticality for two qubits is shown to persist for long distances. Next, we show that the symmetry-unbroken thermal ground state also possesses magic content between distant qubits with scaling behaviour near the critical point. We then consider the magic content in the thermal states and show that thermal fluctuation eventually wipes out the magic for finite temperatures. As an application of quantification of magic in such spin systems, we furnish evidence that the quantum critical region at finite temperature can be detected by the reduced two-qubit robustness of magic (RoM).
The structure of the paper is as follows. In section 2, we review some preliminary materials. After reviewing the physical model of anisotropic XY spin chain with transverse magnetic field, we provide a brief introduction to the stabilizer computation paradigm and magic as a quantum resource. In section 3, we discuss how magic content in the spin chain behaves near the critical point for the symmetry-broken ground state. In section 4, we show the results for the states in thermal equilibrium. We finally conclude with a brief discussion in section 5.

Preliminaries
We start with a brief overview of the many-body Hamiltonian and its properties at thermal equilibrium that are analytically obtainable. This is followed by a discussion on magic as a quantum resource in the context of stabilizer computation paradigm and its quantification relevant to the many-body model at hand.

Many-body system
We consider a chain of N qubits with the transverse field anisotropic XY Hamiltonian [57,58], where J is the exchange coupling, γ is the anisotropy parameter, h is the magnetic field strength, and σ α i denotes the α-th Pauli operator, with α ∈ {x, y, z}, acting on qubit i. As λ = J/h varies, the ground state exhibits a second order quantum phase transition from an ordered to a disordered phase with the critical point at λ = λ c = 1. The U(1)-symmetry broken ground state of this model is factorizable if λ = λ FGS = 1/ 1 − γ 2 [24]. This model can be solved analytically in the thermodynamic limit (N → ∞) [59,60]. The quantum phase transition can be captured by the longitudinal magnetization σ x as the order parameter of the model, which is non-zero in the ordered ferromagnetic phase for λ > 1, and vanishes in the disordered paramagnetic phase. Notably, the condition γ = 1 corresponds to the transverse Ising model, and the condition γ = 0 corresponds to the isotropic XY chain. The reduced density matrix of a qubit at site i is written as α=0,x,y,z σ α i σ α i , where σ α i is the average of σ α i in the ground state (σ 0 = ). The two-site density matrix of qubits at sites i and j is written as where σ α i σ β j is the average taken in the ground state. In the thermodynamic limit, the two-point correlation functions between two qubits at distance r apart for x-and y-directions are given by the determinant of the Toeplitz matrices in the following way [60,61], where G r = 1 π π 0 dφ cos(φr)(1 + λ cos φ) 1 , and ω φ is given by By assuming the two qubits to be far distant, that is, r → ∞, the longitudinal magnetization σ x can be found in a compact analytic form [60] where g(λ) = 1 in the ordered phase, and vanishes in the disordered phase (see [14] for detailed discussions), and β x = 1 8 . As the density matrix is real, σ y i vanishes for all λ and γ. The transverse magnetization σ z is given in terms of elliptical integrals [60], The two-point correlation function in z-direction is given by These magnetizations and two-point correlators enable one to write down the reduced density matrices of the spin chain which are then used to calculate various objects of interest such as the magic content in the current work. We now go on to provide an introduction to magic as a quantum resource.

Magic state formalism
A universal quantum computer is capable of implementing every unitary transformation on the Hilbert space of N qubits. An important subset of unitary operators is the set {U stab } that stabilize the Pauli operators under their action, such that U stab σ α i U † stab = σ β j . These stabilizer unitaries are the key ingredients for quantum error correction [62]. Although the stabilizer unitaries can generate entanglement [63], they are not universal since their action on a particular input bit-string can only cover a subset of the entire Hilbert space. In fact, a stabilizer circuit is demonstrably efficiently simulable with a classical probabilistic computer with the required depth of the classical circuit scaling polynomially with the desired precision [64], in contrast to an exponential scaling for a generic unitary operation. Remarkably, the injection of some extra quantum resources, known as magic states, allows one to perform universal quantum computation, even with stabiliser circuits [29]. Note that, magic is distinct from other quantum mechanical resources, e.g., entanglement. For instance, product states, or even single qudit states, can be magic states, while certain highly entangled states, e.g., cluster states, contain no magic [63]. In this context, the free, i.e., non-magical states, are called stabilizer states, which are states generated from the action of stabilizer unitaries on |0, 0, . . . , 0 , and any convex mixture. For single qubits, the stabilizer unitaries are the Pauli operators together with the Hadamard and the phase gates [34], whose actions on |0 generate six pure stabilizer states. Therefore, the stabilizer states form an octahedron inscribed within the Bloch sphere [34,65], and states outside the octahedron represent magic states. For multiqubit systems, entangling gates such as CNOTs may also belong to the set of stabilizer unitaries.

Quantification of magic
Quantification of any quantum resource is usually possible in terms of several different monotones. For magic, the most obvious candidates are the distance based measures, e.g., relative entropy of magic, or the trace distance of magic, which have been demonstrated to satisfy monotonicity under stabilizer operations [34]. However, computations of these quantities involve optimizations over the entire Hilbert space, and hence become prohibitively hard for all but very small systems. We mention in this connection that in case of an arbitrary number of odd dimensional qudits, magic monotones based on the negativity of the discrete Wigner function do not require any optimization [34]. However, in our case, where constituent spins are two-dimensional, such monotones do not exist.
Thus, to efficiently quantify magic, we use the RoM [35,37,65], defined for an arbitrary quantum state ρ as, where S is the set of stabilizer states. The reader may note that the formulation of the RoM is in terms of a semi-definite programme, which makes it far more efficient to compute than distance based magic monotones. The RoM R(ρ) quantifies the minimal weight of stabilizer states which, upon mixture with ρ, yields a stabilizer state. Using the overcompleteness of stabilizer states, one may express ρ as a pseudomixture of stabilizer states S k ∈ S such that ρ = k X k S k , where the weights {X k } are in general, arbitrary real numbers satisfying the normalization constraint k X k = 1. For a stabilizer state ρ, it is possible to find at least one such decomposition where all {X k }s are non-negative. However, if ρ is not a stabilizer state, at least one of the weights X k is negative. The RoM can be alternatively expressed as the following linearized minimization problem [36], where A αβ = Tr(σ α S β ) and B α = Tr(σ α ρ), with σ α as the αth Pauli operator. This linear programing problem can be solved numerically via any convex optimization package. It is worth emphasizing that the maximum RoM for a single qubit confined to one of the equatorial planes of the Bloch sphere is √ 2 − 1, attained by the H states [37,65]. For the spin chain we are working with, the single-qubit RoM can be written in the following simple form (see appendix A) In the extreme regime of the disordered phase, i.e., λ → 0, all the spins point towards the magnetic field, and show no magic. Deep in the ordered phase λ → ∞, the ground state takes a GHZ form

Magic in symmetry-broken ground state
We first engage in computing the RoM in the ground state of the spin chain where the global phase flip symmetry (invariance of the Hamiltonian under the action of the operator U = i σ z i ) is spontaneously broken. In a physical realization this usually occurs due to small external perturbations. Therefore, in the ordered phase, where there is a two-fold degeneracy, we choose the ground state with positive order parameter, without loss of generality in the results obtained. The analyses of the behaviour of RoM near criticality for single-qubit and two-qubit subsystems are presented in the following.

Single-qubit magic analysis
The reduced density matrix for a single qubit in the symmetry broken ground state is 1 2 . When λ is close to the critical point λ c in the ordered phase, the order parameter σ x , given by equation (3), scales approximately as where K x is a constant, expressible in terms of the anisotropy parameter γ [14]. The derivative of the transverse magnetization has a logarithmic divergence near criticality, but σ z also can be numerically approximated by the algebraic behaviour, where σ z c is the transverse magnetization at the critical point, and β z , K z are γ-dependent constants. It is noteworthy that β z β x for all possible values of γ ∈ (0, 1] (see appendix B).
For every anisotropy parameter γ, one can locate a value of λ = λ * c (γ), such that the RoM vanishes for every λ λ * c , and finite thereafter. We call this point λ * c (γ) as the magic pseudocritical point (MPP). Although λ * c (γ) is very close to the critical point, it always lies in the ordered phase. Note that the existence of an MPP is a consequence of the definition of RoM in equation (8), and does not imply a new criticality. Nonetheless, we call λ * c the MPP, since R γ (λ) exhibits power law scaling behaviour in the vicinity of λ * c , as shown below. By inserting the algebraic behaviours of σ x and σ z in equation (8), and by using the where δλ c = λ * c − λ c is the distance between the MPP and the critical point. The above result is valid for λ > λ * c . When the parameter λ is in very close vicinity of the MPP, determined by λ − λ * c δλ c , keeping terms up to first order in (λ − λ * c )/δλ c in the expression of RoM in equation (9) results in the following linear scaling behaviour about the MPP where the prefactor is given by Beyond this regime, when λ − λ * c > δλ c , the RoM has contributions from two algebraic behaviours, as given in equation (9), up to a constant which depends on γ. As β x is smaller than β z the dominant behaviour of the RoM stems from the first term in equation (9) We note that this expansion is about the critical point, and as such, valid for all λ > λ * c reasonably near criticality.

Transverse Ising chain
We first focus on the transverse Ising chain (γ = 1) in the thermodynamic limit. Using the analytic form of σ x from equation (3) and numerically computing σ z from equation (5) one can compute R γ=1 (λ) using equation (8). In figure 1(a), we plot the RoM R γ=1 (λ) as a function of the control parameter λ. As figure 1(b) shows, the rising of RoM from zero, namely MPP, is very close to the critical point, and takes place at around λ * c = 1.000 15. After the MPP, the RoM rises very steeply and reaches its maximum around λ = λ max = 1.13, and then decays slowly as λ increases further. Since the MPP is very close to the critical point, the region for linear scaling of magic, as given by equation (10), is very small, and displayed in figure 1(c). In figure 1(d), we plot the derivative ∂ λ R γ=1 (λ) as a function of λ − λ c in the log-log scale, which shows algebraic behaviour of the form ∂ λ R γ=1 (λ) ∼ (λ − λ c ) −μ , with μ = 0.88. The fitting parameter μ is very close to 1 − β x (with β x = 1/8), giving an excellent agreement with the prediction of equation (11).

Anisotropic XY chain
We now consider the more general case of the anisotropic XY chain, where γ ∈ (0, 1]. In figure 2(a), we plot R γ (λ) as a function of both anisotropy parameter γ and the control parameter λ. Again, for a fixed value of γ, as λ varies, there is a sharp rise in the RoM right after the MPP until it reaches its maximum R max γ = R γ (λ = λ max ), and then decays gradually. In figure 2(b), we plot the MPP as a function of γ. As the figure shows, λ * c remains very close to λ c , and monotonously moves towards λ c = 1, as the anisotropy in the Hamiltonian decreases. Interestingly, the deviation δλ c shows power law scaling with γ, which is, δλ c ∼ γ 5.55 .
Approaching the isotropic limit as γ → 0, magic in the system vanishes for any value of the control parameter λ, as the longitudinal magnetization σ x vanishes throughout. The corresponding scaling exponents μ is plotted for various anisotropy parameters γ in figure 2(c) which shows small variations near 1 − β x . The slight variation of μ across the phase diagram is due to sub-dominant corrections from σ z . In figure 2(d), we plot the maximum RoM R max γ as a function of anisotropy parameter γ. Interestingly, the R max γ peaks at γ = γ 0 = 1/3, with the corresponding λ max = λ 0 = 1.06, and reaches its maximum value of √ 2 − 1 (up to the set accuracy of numerical evaluation), which is the maximum attainable magic from a qubit confined to an equatorial plane of the Bloch sphere [37,65]. In figure 2(e), we show that the globally optimal magic in the parameter space (λ, γ) is created when the ground state is factorized [22][23][24]. This is a particularly remarkable result, since the single-qubit FGS is pure, the ground state for this parameter value is demonstrably a pure H-state. A large number of H-state qubits are therefore obtainable at this point in the phase diagram which we denote as (λ 0 , γ 0 ).
We emphasize here that the above quantum many-body system provides excellent raw materials, namely, pure unencoded magic states, for the magic state injection based paradigm of quantum computation. We do not address the question of overheads in a realistic magic state factory [29,36]. While it is possible to create such unencoded magic qubit through, for example, the non-interacting paramagnetic Hamiltonian j − J σ x j + σ z j , the claim of this work is the demonstration that such unencoded pure magic states can be obtained in a more realistic interacting many-body system by focussing on a point in the phase diagram, namely the FGS point, which was previously not strongly considered from a practical standpoint due to lack of other established operational resources like entanglement. We check the stability of our result against small imperfections in the tuning of the Hamiltonian by computing the single-qubit purity and its fidelity with H-state for small deviations from (λ 0 , γ 0 ). Even a relatively large ∼ ±2% error in tuning the external magnetic field or ∼ ±10% error in tuning the anisotropy results in the loss of fidelity of less than ∼ 0.1%.

Finite size scaling
In practice, all systems are finite, for which analytic results are difficult to obtain, since the approach of finding the magnetization σ x from the limiting case of the corresponding two-site correlation functions fails. We use numerical methods based on density matrix renormalization group (DMRG) technique [66][67][68] using the ITensor Library [69], with bond dimension 300, and magnitude of symmetry-breaking field ∼ 10 −8 h, and consider the central site to minimize the boundary effects. We specifically perform the finite size scaling analysis to extract the critical exponents, some of which have been directly calculated in the previous sections. Inspired by equation (11), we suggest a finite size ansatz [70] for the derivative of magic as is an arbitrary function and λ [N] c is the finite size critical point at which the derivative of the order parameter σ x peaks. In figures 3(a) and (b), we plot c ) for various system sizes for the transverse Ising case and for γ = 0.5, respectively. By tuning μ and ν, we collapse the curves corresponding to such systems. The best collapse is achieved with μ = 0.88, ν = 1.00 for the transverse Ising case, and with μ = 0.86, ν = 1.09 for the case when γ = 0.5. In both cases, the exponents are quite close to the value of the scaling exponent μ extracted from the infinite chain behaviour, as demonstrated in figure 2(c).

Two-qubit magic analysis
We now go on to analyse the RoM for two qubits at arbitrary distance r in the symmetry-broken ground state. Unfortunately, computing the two-qubit reduced density matrix requires the knowledge of σ x 0 σ z r , which is analytically intractable. Therefore, we again resort to numerical techniques based on infinite DMRG algorithm [71]. The two-site reduced density matrices are obtained using the ITensor Library [69], with the same maximum bond dimension as before and a small magnetic field in the x-direction to break the symmetry. The two-qubit RoM is computed by optimizing equation (7), using the convex programme solver CVX [72,73]. As shown in figure 4(a), the two-qubit RoM is present for arbitrary distances between  The maximum value of global magic Q R at MPP stays decreases very slowly with increases in inter-site distance, with larger rate of decrease observed for more anisotropic spin chains. The inset shows similar trend displayed by the relative entropy between the reduced two-qubit state and product state of two single-qubits for the transverse Ising case. the qubits, contrary to the short-ranged behaviour of two-qubit entanglement. It increases sharply as one goes from the disordered to the ordered phase. This steepness of the growth at criticality increases as the distance between the two qubits increases, and should reproduce similar qualitative behaviour as displayed by the single qubit in the limit of infinite distance. Expectedly, the maximum two-qubit RoM value of (3 √ 2 − 2)/3 [37] is attained at FGS point (λ 0 , γ 0 ), as shown in figure 4(a). At this point, a natural question to ask is whether the two-qubit RoM reveals anything about the correlations in the symmetry broken ground state. For single-system quantum properties like coherence or entropy, a possible measure of correlation for a bipartite quantum system can be put forward as the difference between the magnitude of that property for the bipartite state, and that for the uncorrelated state which is a product of reduced individual states of the bipartite system. To wit, if C(ρ 12 ) is the value of such a property for a bipartite quantum state ρ 12 , then the following quantity C global may tell us about the correlation in the bipartite state C global = C(ρ 12 ) − C(ρ 1 ⊗ ρ 2 ).
If C is the quantum coherence, then this quantity is known in literature as the global coherence [74]. If C is the von-Neumann entropy of a state, then this quantity is simply the mutual information (MI). In our analysis, we assume C to be the log-RoM [35,36], that is, C = log(1 + R), and concentrate on the global part of the magic defined on a two-qubit state ρ 12 as This quantity is certainly nonzero for a bipartite correlated state, and its invariance under local Clifford unitaries can be easily shown [34]. Plotting this quantity in figure 4(b) reveals several interesting features. Firstly, we know that any quantum system near criticality is scale invariant, and is accompanied by presence of long range correlations. For entanglement, its derivative diverges at the criticality, but the absolute magnitude does not attain an extremum. Moreover, since bipartite entanglement does not persist beyond the next nearest neighbour in this model, it cannot be used when searching for presence of long range quantum correlations [14]. Discord and similar quantities also have a diverging derivative at the critical point and not an extremum [75]. However, the global magic Q R attains its maxima very close to the critical point, and in fact the point at which the maxima is obtained is exactly the MPP, irrespective of the inter-site distance r. This is in stark contrast with the case for other quantum correlation measures like entanglement, or quantum discord, which exhibit maxima for different values of the order parameter depending upon the inter-site distance [14,76,77]. This endows the MPP with another physical property. It is also interesting to note that at the MPP, the global magic shows a non-analyticity in the sense that the derivatives from the left and right do not match. From figure 4(c), we also observe another interesting fact, that is, the value of global magic Q R at MPP in particular, and near criticality in general, is a slowly varying function of the distance between the qubits. This supports our claim that the global magic Q R captures long range correlations in critical quantum systems very well. In the scenario of considering pairwise qubits, a particular manifestation of the long range correlation present at criticality is observed when considering the geometric distance between the correlated two-qubit reduced state and the uncorrelated product state composed of local single-qubit density matrices. This distance, quantified in terms of relative entropy (or, fidelity), does not decay sharply in the spin chain at criticality, which can be seen in the inset of figure 4(c) for the Ising chain. However, more familiar measures of correlation fail to capture this feature, as may be observed in figure 4(b) for the case of MI, entanglement, or discord. On the contrary, the global magic Q R successfully captures the long range correlation present among qubits in the critical transverse field XY model. It is beyond the scope of the present investigation to delineate the contributions from quantum or classical nature of the correlations.

Magic at thermal equilibrium
After studying the behaviour of RoM in the symmetry-broken ground state of the spin chain, we now look at the features appearing for states in thermal equilibrium. These states now retain the phase-flip symmetry, since the presence of thermal fluctuations allow degenerate ground states to eventually become equally populated. RoM for the single qubit would be always zero as σ x = 0. Therefore, in this section we will only focus on two-qubit RoM.

Thermal ground state
We start with the zero-temperature limit where the density matrix is the pure ground state in absence of degeneracy, and an equal mixture of the two degenerate ground states otherwise. Taking advantage of the global phase-flip symmetry and the additional symmetries of the problem as well, namely, the translational invariance and the reality of the Hamiltonian, one can write the reduced density matrix for the sites i and i + r as, Using the expressions for the single-and two-site correlations from section 2 we compute the density operator, and subsequently determine the RoM R(γ, λ, r) for two qubits at a distance r by optimizing equation (7) as before.
In figure 5(a) we show result for the transverse Ising case i.e., γ = 1. A closer look at this figure reveals a diverging nature of ∂ λ R near criticality which is studied next and is shown in figure 5(b). Indeed, we find a logarithmic divergence, The scaling exponent μ is positive for r = 1 and negative for r > 1, as is shown with the help of the numerical fits in the inset of figure 5(b) for r = 1 and r = 10. We now calculate the two-qubit RoM R in the γ-λ plane for various distances. The results are qualitatively similar to that for the transverse Ising case. The peak of R is usually reached in the ordered phase, except for the cases of γ > 0.56 where the peak for only r = 1 is found in the disordered phase. The divergence at criticality is analysed again, and absolute value of the exponent |μ| is found to increasing with r for a given γ. Moreover, we found a fit for such behaviour: |μ| ≈ c 1 tanh(c 2 r)r 0.8 , where c 1 , c 2 are constants that depend on the anisotropy parameter γ. This is depicted in figure 5(c). Another interesting kind of scaling appears if we consider the magic rising point (λ MRP ) of the system which we define as the lowest value of λ for which R is non-zero. From figure 5(a), it is evident that MRP increases with increasing r. In figure 5(d) we show that the distance between λ MRP from the critical point goes as lnr, for large values of r. We note in passing that the two-qubit discord also shows persistence over long distance for the transverse XY model. However, for quantum discord, the first derivative does not show a pronounced peak at criticality, and a logarithmic divergence at criticality is observed only for the second derivative with respect to the order parameter λ [76]. In this respect, two qubit magic for arbitrary distances behaves rather similarly to nearest neighbour entanglement than discord, in the sense that the first derivative of the two qubit magic does show a logarithmic divergence at criticality [13].

Thermal behaviour at criticality
So far we have focussed our attention on the quantum critical point at T = 0. In the finite temperature scenario, thermal noises would dominate at higher temperatures and a classical phase transition would be  [78]. The effect of quantum fluctuations is to be captured by the deviation of the Hamiltonian parameter λ from unity, while the strength of the thermal fluctuations is indicated by the temperature T. The quantum critical region in the (T, λ) parameter-space is supposed to fan out conically from the single quantum critical point at (T, λ) = (0, 1). From the perspective of pairwise quantum entanglement, this effect was indeed observed in [28,79]. In the final part of the present work, we demonstrate that two-qubit RoM in the spin chain can also detect the existence of the quantum critical region, retaining largely the same qualitative features. As before, the fact that pairwise magic persists for long distances, unlike entanglement, which is absent in qubit pairs more than two sites away, may present usefulness in a less constrained characterization of the quantum critical region.
Let us first concentrate on the effect of thermal noise on magic content at the critical point λ = λ c = 1. In addition to the significance in characterizing the nature of thermal fluctuations for finite temperatures, this is obviously important from the operational standpoint as well, provided quantum computing gates utilizing these qubits as magic state ancillae are being considered. We demonstrate in figure 6(a) that the magic content of the two qubits decreases monotonically with increasing temperature, and vanishes at a certain finite temperature. The situation is roughly similar to that of nearest neighbour entanglement except for one detail. At very low temperature, the entangled ground state may get mixed with an even more entangled excited state, thus leading to an increase in entanglement vis-a-vis the T = 0 scenario, a phenomenon that is often called entanglement mixing [13]. However, for the case of bipartite magic, we could not find any evidence of such magic mixing, and the quantity of magic monotonically decreases with temperature, finally vanishing at a sudden death point at some finite temperature. It is noticeable that the temperature at sudden death point monotonically decreases as the distance between the two qubits is increased. Similar results are obtained if more general transverse XY chains are considered.
The above observation is quantified in figure 6(b), where we demonstrate that, especially for larger distances, the sudden death temperature T c at criticality follows a scaling law with the distance between the two qubitss r, expressed as where κ is a function of the anisotropy parameter γ of the transverse XY chain. We also demonstrate in the inset of this figure that κ decreases slowly as the anisotropy of the transverse XY chain is increased, and its magnitude is maximum for the Ising case. Low temperature scaling. We have previously demonstrated that at zero temperature, the derivative of magic with respect to the order parameter λ diverges at λ = λ c = 1, signifying the quantum phase transition. It is thus a natural question to ask what happens to the quantum fluctuations, signified by the behaviour of the derivative of magic at criticality, if the temperature T being considered is finite, yet lower than the sudden death temperature of magic T c . We observe the following general behaviour of the derivative of magic in these cases where the constants c and ξ depend on the anisotropy parameter γ, as well as the distance r between the two qubits. This is similar to the behaviour of entanglement [14]. Figure 6(c) shows this for the transverse Ising case, although this is observable for general anisotropy parameter values.

Identifying quantum critical region through pairwise magic
Quantum criticality occurs at T = 0, which is unreachable in principle. However, even at finite temperatures, it is possible to surmise the existence of quantum criticality through a crossover phenomenon between the long-range ordered, i.e., renormalized-classical, and the disordered regimes. At low temperatures, in addition to quantum fluctuations of the order parameter σ x , the thermal noise due to temperature becomes important. For T = 0, there are only two regions in the phase diagram, namely the long range ordered region associated with λ > 1 and the disordered regime associated with λ < 1. However, for finite temperatures, these regimes cross over into each other through the so-called 'quantum critical region', characterized by the crossover temperature T cross = |λ − λ c | zν , where z = ν = 1 for the Ising universality class [80]. If the temperature is much lower than T cross , the thermal excitations, characterized by the thermal de-Broglie wavelength, are not sufficient to induce transitions from ground to excited states on their own. Hence the correlation function separates out into two contributions coming from quantum and thermal fluctuations, respectively. On the other hand, when T T cross , the thermal wavelengths are of the same order as the spacings of excitations, and such a separation can no longer be maintained. The latter regime is the quantum critical region, where the physics is dominated by the interplay between these two fluctuations.
By the theory of quantum critical phenomena [80], the scaling ansatz for the derivative of magic with respect to the order parameter λ for two qubits at sites r-distance apart is given by and the derivative of magic with respect to temperature, signifying the thermal fluctuations are to obey the following scaling ansatz where f 1 and f 2 are the scaling functions and the exponent n is non-universal. Function f 2 is characterized by a maxima that collapses for all those points in the (λ, T) phase diagram, whose distance from quantum critical point is given by d = √ T 2 + (λ − λ c ) 2 . The maxima, at the crossover temperature T * = aT cross , is verified in figure 6(d) to be almost independent of the distance r between the two qubits in transverse Ising chain. For example, at r = 10 (shown in figure 6(d)), a ≈ 0.55, and for r = 1 (not shown), a ≈ 0.56. Therefore the magic experiences the maximal variation in thermal fluctuation at T = T * , signifying the crossover phenomenon.
We have previously mentioned that the quantum critical region is characterised by the interplay between the thermal and the quantum fluctuations, even when T T cross . The strength of these fluctuations, characterized by the Grüneisen parameter [81] or the ratio of the absolute values of ∂ T R(r), and ∂ λ R(r), is depicted in figures 7(a) and (b). As is expected for the quantum critical region, the ratio is very small, signifying the dominance of quantum fluctuations for temperatures far in excess of the crossover temperature T cross , which is the signature of the quantum critical region.
Let us now concentrate on the response of the quantum fluctuations upon variations in temperature, i.e., ∂ T ∂ λ R(r). In the quantum critical region, since T T cross , the thermal de Broglie wavelength is of the same order as the spectral gap, and the system responds to any infinitesimal change in the temperature. We observe, similar to entanglement, that this second derivative attains an extremum along the line T M = bT cross . For a representative (r = 10) case, we observe in figure 7(b), that for T > T * , any reduction in the temperature T leads to the increase of the effects of the criticality, while for T < T * , the reduction in temperature leads to a net decrease in the effects of the criticality. This is along similar lines to behaviour of entanglement. Another interesting aspect is the observation that as one moves away from the critical point, the quantum critical region identified by magic in nearest neighbour qubits (r = 1), is skewed in a different direction than the critical region identified by magic in more distant neighbours owing to the opposite signs of the derivative of magic near criticality for nearest neighbour qubits. These, in conjunction, serve to capture the quantum critical region in detail. We add that although the results shown in this subsection pertain to the transverse Ising chain, they similarly hold for differing magnitudes of anisotropy for the more general transverse XY chain. Similar to entanglement, magic is a purely quantum feature of a system, hence our result supports the claim of reference [28] that purely quantum effects like entanglement serves to characterize the crossover between quantum critical and renormalised classical regions better than merely classical effects.

Conclusion
In this paper, we have studied transverse field anisotropic XY model near criticality for single-and two-qubit magic content. In the physically relevant symmetry broken ground state, we quantify scaling behaviours of the single-qubit RoM and identify that the factorization point in the phase space provides perfect unencoded magic states. This endows the FGS with an operational significance in terms of presence of a resource for fault tolerant quantum computation for the transverse XY model. It will be interesting to investigate in future whether this phenomenon of FGS being associated with maximal magic content is persistent for other quantum many-body systems. Furthermore, we show that, similar to discord, the two-qubit RoM persists over long distances, which is a direct consequence of the presence of long-range correlation in this critical system. We show that the sharp maximum attained by the global component of magic is located at the exact point separating the magical and non-magical regimes irrespective of the inter-site distance, which is extremely near the critical point (at λ ≈ 1.000 15). These properties are in marked contrast with the smooth peak observed for entanglement and discord. In the symmetry unbroken case, which is relevant for the thermal states, scaling behaviours have been identified for two-qubit RoM in the zero temperature limit. The effect of the thermal fluctuation is shown to be destructive resulting in a sudden death of RoM. The competition between quantum and thermal fluctuations extends the quantum critical point into a critical region which can also be identified using RoM of two distant qubits.
With our understanding of the behaviour of RoM in the anisotropic spin chain, we can now compare magic with other established resources, namely, entanglement, discord, and coherence, in this system, which  [82,83] and magic show long range survival, in contrast to bipartite entanglement which does not survive beyond the next nearest neighbour. However, logarithmic divergence at criticality is shown by the second derivative of discord, compared to the first derivative in case of magic. The similarity of behaviour of magic with entanglement is observed in terms of thermal sudden death, which is not as pronounced for discord and coherence. It may be noted in this context that although discord does not undergo sudden death in this spin system, a sudden transition induced by environmental interaction during Markovian evolution has been reported [84][85][86]. In the symmetry broken case, any correlation including discord and entanglement is absent at the FGS, and the single site quantum coherence is finite but non-maximal. In contrast, the single site magic attains its global maximum at the FGS. This study can be extended for other many-body systems and there are various avenues to explore including dynamical phase transitions and open system evolution. We add in this context that in previous years, there has been a growing interest in dynamical evolution of quantum properties in spin systems, as well as quantification of magic content of quantum channels [87]. While we focus on the static properties of magic in the spin chain in this work, a thorough study of open system evolution of magic content can shed further light on how decoherence affects fault tolerant quantum architectures in realistic scenarios. Figure B1. Scaling of deviation of transverse magnetization σ z from its magnitude at λ c , with deviation from criticality λ − λ c for an infinite transverse XY chain with different γ (points), and the exponent β z extracted from the corresponding linear fits (solid straight lines).
Now, let us note that only the coefficients of |1 1|, and |− −| may be negative here. If μ < 1, then the new decomposition with two negative coefficients leads to a lower RoM value of 2|1 − μ| + 2 , which is less than the assumed optimal decomposition if is arbitrarily small. If μ 1, then the state lies within the stabilizer polytope since all the coefficients are now positive. Hence the proof is complete.

Appendix B. Scaling of transverse magnetisation
In this subsection, we present numerical evidence that the transverse magnetization has an algebraic behaviour near criticality, i.e., σ z σ z c + K z (λ − λ c ) β z . (B.1) As is depicted in figure B1, the algebraic behaviour of transverse magnetization close to the critical point is clear. The scaling exponent β z broadly lies in the range ∈ (0.8, 0.9). The detailed values of scaling exponents β z for different values of anisotropy, along with the fitting lines are shown in the corresponding figure here.