Many-body magic via Pauli-Markov chains - from criticality to gauge theories

We introduce a method to measure many-body magic in quantum systems based on a statistical exploration of Pauli strings via Markov chains. We demonstrate that sampling such Pauli-Markov chains gives ample flexibility in terms of partitions where to sample from: in particular, it enables to efficiently extract the magic contained in the correlations between widely-separated subsystems, which characterizes the nonlocality of magic. Our method can be implemented in a variety of situations. We describe an efficient sampling procedure using Tree Tensor Networks, that exploits their hierarchical structure leading to a modest O (log N ) computational scaling with system size. To showcase the applicability and efficiency of our method, we demonstrate the importance of magic in many-body systems via the following discoveries: (a) for one dimensional systems, we show that long-range magic displays strong signatures of conformal quantum criticality (Ising, Potts, and Gaussian), overcoming the limitations of full state magic; (b) in two-dimensional Z 2 lattice gauge theories, we provide conclusive evidence that magic is able to identify the confinement-deconfinement transition, and displays critical scaling behavior even at relatively modest volumes. Finally, we discuss an experimental implementation of the method, which only relies on measurements of Pauli observables.


I. INTRODUCTION
Over the last two decades, quantum information concepts have revolutionized the way we understand and approach the many-body problem [1].Remarkable insights on quantum matter have been obtained under the lens of entanglement, a measure of separability that has found applications over a wide range of phenomena, from real time dynamics [2], to topological order [3,4] and classification of states [5,6].A pivotal role in establishing these applications has been played by the development of trustful entanglement measures [7], in combination with efficient theoretical methods to explore that in the context of many-body systems [8] -one paradigmatic example being tensor networks [9,10].
On par with entanglement, another quantum information concept that is receiving increasing attention is that of non-stabilizerness, also known as magic [11][12][13].In the context of quantum computing, magic is now understood as a fundamental resource that would be required to outperform classical simulations, and its concrete role in digital simulations has been widely addressed [14][15][16][17][18][19].However, differently from entanglement, there is presently limited understanding of how magic reflects many-body phenomena, and even if it does it at all [20]: a fundamental limiting factor is that, oppositely to entanglement, we lack an array of scalable, efficient methods to actually compute magic -a shortage that severely limits our capability of identifying situations where there can be a direct connection between magic and physical phenomena.
In this work, we present a theoretical framework to measure many-body magic that leverages on a stochastic sampling of the system wave function.Our work builds upon recent developments in the field, in particular, on the recognition of stabilizer Renyi entropies (SREs) as measures of magic (including an experimental demonstration with 4 qubits) [21][22][23][24].While a direct measure of the former is extremely challenging as it requires a number of measurements that grows exponentially with the size of the partition, we introduce a Markov chain on Pauli strings as a tool to distill the most relevant contribution to magic.We show that our protocol returns an unbiased estimator of SREs of all orders, and that it is efficient under several important scenarios: those include both full state magic (that is relevant, e.g., to quantify the overall difference from a stabilizer state), and long-range magic -a quantity that is akin to mutual information and that, crucially, is not plagued by any UV-divergences when applied to field theory.
The estimation of magic via Pauli-Markov chain is a general construction, that is broadly applicable to computations as well as experiments.We explore in detail its capabilities in the context of tree tensor networks (TTN) [25,26].At first, we perform extensive methodological checks, in particular, on the efficiency of Markov sampling and autocorrelations.We then showcase the flexibility of our approach with several applications, to understand advantages and overall comparison with recently introduced direct sampling methods that constitute the state of the art in terms of measuring many-body magic in numerical computations [27][28][29].
Firstly, we consider one-dimensional systems.There, by considering both Ising, Potts and Heisenberg models, we show that full-state magic is not always indicative of quantum critical behavior.In particular, while it works for the conceptually simple cases of Ising (as already ob-served in Ref. [27][28][29][30]) and Potts models, it spectacularly fails detecting any criticality in the case of spin-1 XXZ models.Oppositely, long-range magic (whose computation was not accessible before our algorithm, to the best of our knowledge) displays sharp signatures of critical behavior in all models considered.Our work thus clarifies how, in the context of critical behavior, it is fundamental to construct -and to compute -UV-divergence free estimators to understand the role of magic.
Secondly, we consider two-dimensional interacting systems, where the connection between magic and manybody phenomena is uncharted territory.We focus on the Z 2 lattice gauge theory, for two reasons: its importance as a paradigmatic model for more complicated lattice field theories, as it displays a confinement-deconfinement transition, as well as topological order; and its direct connection to the toric code, an epitome example of quantum memory based on the stabilizer language [31][32][33][34][35][36][37][38][39][40].Thanks to the very modest O(ln N ) size-scaling of our algorithm versus system size N , we are able to consider systems up to 100 spins.Our results show how both confined and deconfined phase have volume-law magic: most remarkably, magic features striking signatures of critical behavior.Close to the transition point, its behavior is akin to that of a Binder cumulant, as magic density displays a crossing as a function of volume, whose functional form is dictated by finite-size scaling theory.Even more remarkably, universal collapses are not only evident at modest volumes, but even at relatively small bond dimensions, signalling that magic might be considerably less affected than other observables by tensor network truncations.At the physical level, our results point out that magic may serve as an order parameter for confinement-deconfinement transitions, even at volumes where other quantities (e.g., order parameters) are of very limited use.
Finally, we give a glimpse of the applicability of our approach to experiments.In that context, we discuss in detail experimental errors as a function of finite sampling, size, and autocorrelations.Our results indicate that the sampling needed to scale to large systems requires very fast repetition rates, which are available in solid state settings, but constitute a challenge for atomic experiments.
The rest of the paper is structured as follows.In Sec.II, we review the basic properties of magic, as well as SREs.In Sec.III, we describe how to sample Pauli strings via Markov chains, discuss the efficiency of various estimators, and detail our implementation with tree tensor networks.In Sec.IV, we present our results on both one-and two-dimensional spin systems.In Sec.V, we detail our experimental protocol, and then conclude in Sec.VI.
Quantum resource theories aim to capture the fundamental aspects inherent in quantum technology.For instance, entanglement is a crucial resource for quantum cryptography and communication.The resource framework for entanglement finds practical application by providing bounds on the efficiency of entanglement distillation protocols.Error-correcting codes play a fundamental role in achieving fault-tolerant quantum computation.These codes enable the storage of quantum information while protecting it from the detrimental effects of noise.
The development of error correcting codes based on the stabilizer formalism -e.g., the toric code -has motivated a resource theory of non-stabilizerness, or magic.Here, we briefly review it, and summarize the main challenges in addressing magic in the context of many-body theory.
We first formally define the notation.We consider a system of N qubits (generalizations to larger Hilbert spaces will be discussed below), with Hilbert space H = ⊗ N j=1 H j .The N -qubit Pauli group P N encompasses all Pauli string operators with an overall phase of ±i or ±1.Mathematically, we define P N as follows: Moving on to stabilizer states, we can establish that a pure N -qubit state falls into this category if it satisfies certain conditions.Specifically, a stabilizer state is associated with an abelian subgroup S ⊂ P N that contains 2 N elements.For every S ∈ S, the stabilizer state |ψ⟩ remains unchanged under the action of S, expressed as S|ψ⟩ = |ψ⟩.Alternatively, we can define stabilizers using Clifford unitaries, which are unitary transformations preserving the Pauli group when conjugated with it, i.e.
The Clifford set C N can be generated using the Hadamard gate, the π/4-phase gate, and the CNOT gate.Notably, stabilizer states are pure quantum states that can be prepared by applying Clifford operations to a canonical trivial state |0⟩ ⊗N .In the framework of resource theory, stabilizer states are considered free states while Clifford unitaries and Pauli measurements constitute free operations.The computation using only free states and free operations can be efficiently classically simulated, whereas universal quantum computation can be achieved through supplying magic (non-free) states.Therefore to enable successful quantum computations, additional techniques are necessary to ensure the fault-tolerant implementation of a universal set of quantum gates.This can be achieved by augmenting the Clifford group with the Toffoli gate or the π/8 phase gate, thus unlocking the potential for universal quantum computation.
In this context, a central task is the quantification of the amount of non-Clifford operations needed to prepare a given quantum state.The properties required to a good measure M of non-stabilizerness are (i) M(|ψ⟩) = 0 ⇐⇒ |ψ⟩ is a stabilizer, (ii) non-increasing under Clifford operations: For many-body systems, previous investigations into magic measures have primarily concentrated on small or weakly correlated systems, leading to a limited understanding of magic in entangled many-body systems.An inherent challenge arises due to the exponential growth of stabilizer states and their increasingly intricate geometric structures as the system size expands.Consequently, the general calculation or numerical analysis of magic measures for large states becomes arduous.In order to enhance our understanding of this phenomenon, several fundamental questions require attention.These include comprehending the extent to which many-body quantum states can exhibit magic, determining the typical amount of magic found in generic states, and developing methodologies for computing the magic associated with manybody states.
From a quantum information viewpoint, the main motivation in understanding and measuring many-body magic stems from its relevance as a resource towards quantum advantage.Recent studies have shed light on the fact that the computational power of a state cannot be solely attributed to its magic density; other characteristics of magic may also play significant roles [20].These properties encompass not only the primary aspect of magic density but also the subleading terms, nonlocal components, topological aspects, and more.Consequently, it is crucial to develop a numerical scheme capable of accessing and analyzing the various features of magic in many-body systems.Such a scheme would facilitate a comprehensive exploration and understanding of the intricate interplay between magic and computational power.This would constitute a major step forward in our endeavor to fully characterize the role played by magic in many-body systems.
Notwithstanding such practical importance, understanding the connection between quantum correlations and physical phenomena is interesting from a broader perspective [20] -especially, given the importance and impact such a connection has had in the case of entanglement.The connection between magic and physical phenomena is presently poorly understood, due to the combined lack of computable measures of magic, and of methods to attack them.
From the point of view of observables, the key result we will exploit is Ref. [21], that demonstrated SREs as a measure of magic (at least in the case of coherent dynamics; in more complicated scenarios, such quantities are not necessarily measures, see Ref. [28]).From the point of view of connection between magic and physical phenomena, three works are serving as a key motivation in this direction [27][28][29].Thanks to the development of novel techniques based on direct sampling of matrixproduct states (MPSs), these works have pointed out strong connections between critical behavior and magic in the context of one-dimensional systems, at precision and volumes never attained before.We will discuss this in more detail over the next section.

B. Stabilizer Renyi entropy
Stabilizer Rényi Entropies (SREs) are a measure of nonstabilizerness recently introduced in Ref. [21].For a pure quantum state ρ, SREs are expressed in terms of the expectation values of all Pauli strings in P N : with d is the local dimension of the Hilbert space of N qudits and P N is the generalized Pauli group of N qudits [41].The SREs have the following properties: [21] (i) faithfulness: , and (iii) additivity: The SREs are thus a good magic measure in the point of view of resource theory, where the free states are defined as the stabilizer states while the free operations are the Clifford unitaries.This definition is a straightforward generalization to general local dimension d from the one given in Ref. [21].For d > 2, the Pauli operators are no longer Hermitian, and thus the expectation values can be complex.In Eq. (3), we take the absolute values of the expectation values | Tr (ρP ) |. Eq. ( 3) can be seen as the Rényi-n entropy of the classical probability distribution: It has the following properties: [21] (i) faithfulness: , and (iii) additivity: The SREs are thus a good magic measure in the point of view of resource theory, where the free states are defined as the stabilizer states while the free operations are the Clifford unitaries.
Moreover, the definition of SREs can be extended to mixed states by properly normalizing Ξ P .For example, for n = 2, the mixed state SRE is given by [21] which can be seen as the Rényi-2 entropy of apart from some offset.Here, the free states are defined as the mixed states that can be obtained from pure stabilizer states by partial tracing [21].Furthermore, the long-range magic can be quantified by where A and B are two separated subsystems (see Fig. 2 (a)-(b)).A similar quantity has been considered previously in the context of mana [42,43] and robustness of magic [44,45].L(ρ AB ) measures how magic is contained in the correlation between the subsystems, and thus it quantifies the degree to which magic cannot be removed by finite-depth quantum circuits [42].Indeed, due to the additivity of SRE, L(ρ AB ) vanishes for a product state ρ A ⊗ ρ B .On the other hand, a non-vanishing value of L(ρ AB ) effectively quantifies the extent of deviation from the additivity in the case of entangled subsystems.The long-range magic is directly reminiscent of mutual information, that has played a major role in characterizing the distribution of both classical information and quantum correlations in many-body systems [1,[46][47][48][49][50][51][52][53][54][55][56][57].On the lattice, the main motivation for looking at functionals such as in Eq. ( 7) is that they are much more meaningful than simple bipartition properties from a field theory standpoint.Indeed, these quantities are expected to be free of UV divergences, and thus solely dominated by infrared, universal properties of the lattice theory.This parallels the f-functions used in field theory [58].
As discussed above, SREs have attracted recent interest due to their computability.The first technique was introduced in [27], which expressed the SREs of integer index n > 1 as the norm of a "2n-replica" MPS.Although this technique yields an exact value of M n within a given MPS, its computational cost scales as a large power of the bond dimension χ, specifically O(N χ 6n ).Thus, although the method is efficient in principle, in practice it can only access bond dimension up to χ = 12, which limits its applicability to investigate many-body physics.
A different approach based on sampling of Pauli strings according to the probability distribution Ξ P was proposed very recently in [28,29].In those works, the Pauli strings are sampled directly via the perfect sampling scheme with Matrix Product State (MPS) introduced in [59].For the case of open boundary conditions (OBCs), the cost scales as O(N χ 3 ), thus enabling access to larger bond dimensions, which opens the door for investigating magic in entangled states.However, as we discuss in more detail in the next section, this method provides only an efficient estimation of M 1 .It has been demonstrated that for 0 < n < 2, M n violates monotonicity under measurements followed by conditioned Clifford transformations [28].Thus, it is important to develop an efficient scheme to efficiently compute M 2 , which also has the nice property of being experimentally measurable [23,60].Furthermore, M 2 is directly linked to the average over the Clifford orbit of entanglement spectrum flatness in an arbitrary bipartition [61] and participation  8), and for single qutrit state defined in Eq. ( 9) (b).
We also note that the aforementioned two methods have inherent limitations when it comes to evaluating magic within a subsystem of a state -for instance, none can access long-range magic.As a result, the existing techniques are unable to provide insights into how magic is distributed within a given state.

Examples
To familiarize with the behavior of SREs in many-body systems, here we provide some examples of SREs in simple wave functions.First of all, we stress that the SREs are basis-dependent, i.e., it is not invariant under local basis change.In particular, the SREs of a single-qubit state may be non-trivial.For example, consider the following one-parameter family of single-qubit states Note that |ψ(π/4)⟩ corresponds to the canonical T-state.The SREs can be computed easily by evaluating the expectation values of P ∈ {I, X, Y, Z}, and then plugging it in Eq. ( 3).The result is shown in Fig. 1 (a).As can be seen, the SREs are non-zero apart from some special points θ = mπ/2 with integer m.Now, the SREs of a product state of N copies of |ψ(θ)⟩ can also be computed straightforwardly, utilizing the additivity property of SRE, For an example of qudit states, we consider the follow-ing family of single-qutrit states Here, |ϕ(2π/9)⟩ corresponds to the canonical qutrit Tstate.We now need to compute the expectation values of 3 2 single-qutrit Pauli operators.To define the Pauli operators, we first define the shift and clock operators for d-level system as where ω d = e 2πi/d , and the addition is defined modulo d.
For qutrits, we have d = 3.The qudit Pauli operators are defined as Computing the expectation values of the Pauli operators in Eq. ( 11), we can compute the SREs of |ϕ(θ)⟩ using Eq. ( 3).The result is shown in Fig. 1 (b).In this case, the SREs are non-trivial apart from some special points θ = m2π/3 with integer m.

III. MARKOV CHAIN MONTE CARLO SAMPLING OF PAULI STRINGS
In this work, we investigate the SREs using Monte Carlo sampling of Pauli strings according to some probability distribution Π P , which only depends explicitly on the expectation values of Pauli strings.For example, for the calculation of M n , we get Π P = Ξ P (Eq.( 4)), while for M2 we have Π P = ΞP (Eq.( 6)).Here we focus on Metropolis algorithm, although other sampling methods, such as heat bath, may also be employed.Since Π P only depends on the expectation value of P , this method is applicable to any numerical methods in which expectation values of (non-local) operators can be accessed, such as exact diagonalization and tensor network methods.Furthermore, this method can also be utilized to experimentally measure SREs (see Sec. V).

Algorithm 1 Monte Carlo sampling of Pauli strings
Input: a quantum state ρ and number of sampling NS Propose a candidate Pauli string P ′ .

6:
Accept the move with probability: min 1, Measure the estimators.8: end for Output: a Markov chain of P with probability ΠP .

A. Algorithm theory
The scheme is summarized in the Algorithm 1.If we sample according to Ξ P , M n can be estimated using the unbiased estimators for n > 1 and for n = 1, where ⟨...⟩ Ξ P is the average over Ξ P obtained with sampling.For n < 1, a better estimation can be done by reversing Eq. ( 12), i.e., where Let us analyze the efficiency of these estimators.
➩ SRE with n = 1.-For n = 1, the variance of M 1 is shown to be at most quadratic in N in Ref. [29].Thus, the estimator for M 1 is efficient.Actually, we can even make a stronger statement, if we make the assumption that the SREs are linear in N , i.e., M α = N f (α) + O(1), where f (α) is a function that does not depend on N .Using the relation [63]: we see that Var(M 1 ) is linear in N .It follows that the variance (standard deviation) of the SRE density, m 1 = M 1 /N , scales as 1/N (1/ √ N ).➩ SRE with n ̸ = 1.-For n > 1, the variance of Eq. ( 12) is given by Now, by second-order approximation Var (log x) ≈ Var (x) /x 2 , we have For n < 1, Then, In both cases, if the SREs grow at most logarithmically in N , the variance grows at most polynomially.Thus, by Chebyshev's inequality, the number of samples needed for a fixed error ϵ is polynomial in N , i.e, the estimator is efficient.On the other hand, if the SREs are linear in N , as is typically the case in many-body systems [27,30,42], the variance grows exponentially with N when n ̸ = 1.Thus, the estimator for M n , n ̸ = 1 is efficient only if the SREs are at most O(log N ).One can also see this intuitively by noting that the quantity being estimated is exponentially small in N when M n is linear, and thus we need exponentially small precision.We note in passing that states with logarithmically growing SREs can arise in many-body systems in the frustrated regime [24].Note, however, that the SREs are typically linear in N .Therefore, using the estimators in Eq. ( 12), the estimation of M n , n ̸ = 1 will almost always be exponentially costly.Nevertheless, the cost typically grows much more slowly than d 2N which is the cost for exact computation.Thus, in practice, using this estimator is still beneficial to extend the system sizes we can study, as we shall illustrate in Sec.IV.Importantly, using Monte Carlo sampling, we are not restricted to sample the Pauli strings according to Ξ P .An alternative approach is to sample Pauli strings according to the probability distribution Π P,n ∝ Tr(ρP ) 2n .We then need to estimate the normalization constant of Π P,n to estimate M n .This is a non-trivial task, equivalent to estimating the partition function, for which a wealth of sophisticated methods have been put forward [64][65][66][67][68][69][70][71][72][73][74][75].
➩ Long-range magic.-Inaddition, we are interested to estimate the long-range magic as quantified by L(ρ AB ) in Eq. ( 7).While we can in principle compute the individual M2 for ρ C , C ∈ {A, B, AB}, this is not optimal, as we have seen that the estimation for M2 is not efficient when M2 grows linearly with N .Moreover, we expect that the leading term of M2 will be canceled out in L(ρ AB ).In this case, it is more desirable to estimate L(ρ AB ) directly, without having to resort to inefficient estimation of M2 .To do this, we first rewrite Eq. ( 7) as follows: where Two widely-separated partitions for the calculation of longrange magic in Eq. ( 7).(c) Subleading term as in Eq. ( 24), as well as a cartoon depicting the increment trick discussed in the main text.
Π P AB ∝ Tr(ρ AB P AB ) 4 , we can estimate W (ρ AB ) by where P AB is decomposed as P AB = P A ⊗ P B .Similarly, we have Therefore, as a byproduct, our scheme can be applied to compute the Rényi mutual information for disjoint subsystems.
➩ Subleading term.-Theprevious scheme can be straightforwardly modified to extract the subleading term in the expansion M n (N ) = D N N + c N [27].Here we consider 1D systems for simplicity.Specifically, the subleading term is approximated by the quantity c N = 2M n (N/2) − M n (N ) (see Fig. 2 (c)), which expands as for n ̸ = 1, where ρ N,N/2 is the density matrix for a 1D system of size N and N/2, respectively.For simplicity, we have assumed translational invariance in Eq. ( 24), but the procedure can be straightforwardly generalized to any system.Here, denoting P = P 1 P 2 ...P N , where P i is a Pauli operator acting on site i in the N −site system, we choose P (1) = P 1 P 2 ...P N/2 and P (2) = P N/2+1 P N/2+2 ...P N .Note that, differently from Eq. ( 22), here we consider two pure states of different sizes N and N/2.For the subleading term in 1D systems, the term inside the log in Eq. ( 24) does not decay exponentially, and thus the estimation can be done more efficiently than the estimation of the leading term in Eq. (12).
➩ Increment trick for SRE.-The extraction of the subleading term in Eq. ( 24) presents an alternative strategy to estimate M n , which circumvents the problem of exponential variance for the estimator in Eq. (12).The key idea is that, if the estimation in Eq. ( 24) is efficient, then we can estimate c N , c N/2 , ..., until the size is small enough that M n can be evaluated exactly.The number of c M 's that needs to be computed scales as O(log N ) (assuming translational invariance).Then, we can determine M n (N ) by considering a proper linear combination of c M 's.This strategy is reminiscent of the increment trick employed in estimation of Rényi entanglement entropies in Quantum Monte Carlo simulations [76][77][78], which considers the difference of Rényi entropies of smaller and smaller regions, to compute the Rényi entropy of a large entangling region with high precision.However, in this case, the form of c N is specifically designed to cancel out the volume-law term of M n , differently from entanglement entropy which exhibits area law.
The above strategy is effective in 1D systems because the subleading term c N is expected to either remain independent of system size or exhibit at most logarithmic growth.However, in higher-dimensional systems, c N may exhibit area-law scaling, leading to growth with size.In this case, more complicated linear combination of M n 's shall be considered to eliminate the area-law term (while, at the same time, also keeping the volume law one vanishing).For example, in 2D systems, the form of linear combination used in extracting the topological entanglement entropy with Kitaev-Preskill [3] or Levin-Wen scheme [4] will cancel both the volume-law and area-law term.It is convenient to partition the system into four subsystems as proposed in [79], which is also suitable with 2D TTN geometry.With this scheme, the estimation of M n (L×L) is reduced to M n (L/2 × L), M n (L/2 × L/2), ..., such that only O(log N ) computations are required, as in 1D case1 .

B. Efficient sampling with tensor networks: the example of tree tensor networks
The probability Π P of a given Pauli string P only depends on the expectation value of P , and thus it is efficiently computable in TTN (or any loopless tensor network [26]).Following the convention introduced in Ref. [25], each tensor in the TTN structure is denoted by the pair of zero-indexed integers [l, n], where l corresponds to the layer index (starting from the top root tensor) and n denotes the tensor at a particular layer l counted from the left (see Fig. 3(a)).Obviously, in this notation, the top root tensor is represented by The algorithm to sample Pauli strings for the ground state of a quantum many-body system is described below.
➩ After performing the adaptive variational groundstate search [25] for a many-body Hamiltonian, we arrive at the TTN representation of the manybody ground state wavefunction |ψ⟩.We start by bringing the TTN into the central canonical form, where the [0, 0]-tensor is the orthogonality center (see Fig. 3(a)).
➩ Given the initial Pauli string P = P 1 P 2 . . .P N , where P i is a Pauli operator at site i, we construct the coarse-grained effective "link" operators O [l,n] at each link iteratively from the physical sites to the top-most links, where at the bottom-most (i.e., the physical) layer these link operators are identified with the Pauli operators (see Fig. ➩ At each sampling step, we either propose a singlesite update P ′ = P 1 . . .P ′ i . . .P N , or a two-site update P ′ = P 1 . . .P ′ i . . .P ′ j . . .P N , following Algorithm 1.The updated sites i and j are chosen randomly. ➩ We observe (Fig. 3(d)) that the effective link operators for P ′ only differ with those of P on the links that lie on the path from the site i (or j) to the root [0, 0]-tensor.The number of such links scales only logarithmically in system size.This implies that computing ⟨ψ|P ′ |ψ⟩ can be done very efficiently with a computational cost of O(log(N )χ 4 ), as opposed to O(N χ 4 ) for a generic many-body operator for the TTN.
The heart of our efficient sampling procedure lies within the above observation for TTN.We exploit this scaling property to perform efficient Monte Carlo sampling of Pauli strings by standard Metropolis algorithm, where the candidate Pauli string for the next configuration only differs at a few sites with the previous Pauli string configuration.Crucially, the sites can be chosen arbitrarily, and this does not change the log N scaling of the TTN sampling, provided that the number of modified sites does not scale with system size.This allows for flexible sampling strategy, which can be designed by taking into account our knowledge about the state that we want to sample -very much like Monte Carlo methods are designed to probe partition functions.
The final step for calculating the expectation value of a proposed candidate Pauli string at each Metropolis iteration is the following.
➩ The link operators, that reside in the path from the updated site i (or j) to the root [0, 0]-tensor, are updated by the coarse-graining step.The expectation value ⟨ψ|P ′ |ψ⟩ is now calculated by tensor contractions of the root tensor and top-most (updated) link operators (see Fig. 3(d)).
. At this stage, it is important to discuss the efficiency of the more widely used MPS tensor network structure in relation to our sampling strategy.The computational cost for direct sampling of Pauli strings using MPS with OBCs scales as O(N χ 3 ) [28,29,59], that also holds for Monte Carlo sampling using MPS2 , as opposed to the O(log(N )χ 4 ) that we get utilizing TTN.Consequently, our method with TTN for obtaining SREs becomes increasingly efficient as the number of qudits N grows large, particularly when N/ log N ≳ χ.Specifically, since the MPS or the TTN bond dimension χ saturates to a constant value with N in 1D quantum systems with gapped spectrum due to the area-law of entanglement entropy, our approach involving TTN vastly outperforms MPS based methods in terms of efficiency for large N .Most importantly, the enhanced connectedness inherent in the TTN structure allows for efficient exploration of higherdimensional (2D and even 3D) many-body systems (see e.g., [80][81][82][83][84]).This paves the way to investigate SREs in higher-dimensional systems, as we present in Sec.IV.
Finally, we mention that our can also be used to compute the SREs of any partition of the system.To do this, we only need to restrict the Pauli strings to have support on the sites in the partition.Using the estimator for n = 2, the same Monte Carlo procedure will yield M2 in Eq. ( 5).Moreover, the algorithm is easily generalized to Tree Tensor Operator (TTO) [85], which represents many-body density operator for mixed states.
We note that the use of Monte Carlo techniques in tensor network has been considered before [86][87][88] to compute the expectation value of a local operator.Instead, here the expectation values are computed exactly, while the sampling is done at the level of operators being computed.

IV. APPLICATION TO QUANTUM MANY-BODY SYSTEMS
We apply the TTN based sampling method in Sec.III B using the estimators in Eq. ( 12) and Eq. ( 13) to investigate the SREs in various many-body systems, especially near quantum critical points, both in 1D and 2D geometries.Unlike MPS, the structure of TTN allows for efficient exploration of systems under periodic boundary conditions (PBC) with similar computational cost as the open boundary conditions [25].Therefore, we consider the periodic many-body systems, i.e., ring and torus geometry in 1D and 2D, respectively, to avoid boundary effects.For the analysis of statistical errors and the autocorrelation times in the Markov chain samples, we refer to the Appendix A, whereas for the analysis of convergence with bond dimension of the TTN, we refer to the Appendix B.
To obtain the TTN representation of the ground state of many-body systems we perform variational minimization with TTN sweeping algorithm [25,26], and then employ the sampling scheme in Sec.III B to estimate the SREs of the ground state.In particular, since the SREs are generally linear in the number of qudits N , we focus on the SRE densities m n = M n /N .
All of the models we consider possess Z n symmetry, with n = 2 or 3, and thus, a two-site update scheme is required to sample only the Pauli strings that preserve the symmetry.The Pauli strings that preserve the Z n symmetry, generated by i Z i , are generated by Z i and X † i X j (up to a phase constant).Here, X and Z are the shift and clock operators defined in Eq. (10) To ensure that only the Pauli strings that obey the Z n symmetry are considered, we generate the candidate Pauli string P ′ by randomly multiplying the current Pauli string P with either Z i or X † i X j .It is easy to see that the update scheme is ergodic.For d = 3, we set the probability to multiply with Z i or Z † i to be equal, so as to satisfy detailed balance.For d = 2, when there is time-reversalsymmetry, the Pauli strings are additionally constrained to those with even numbers of Y = iZX.As such, the Pauli strings with odd numbers of Y can be directly re- The SRE density m2 for the 1D quantum Ising chain near the critical point hc = 1 computed using the increment method using different subleading terms.(Inset) The sampling errors for m2 at hc = 1 for various system-sizes L (in log-scale).Clearly, the errors show even slower than than logarithmic growth for the efficient sampling scheme.Here we consider TTN bond dimension χ = 30 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.

jected.
A. Non-stabilizerness in 1D many-body systems The behavior of SREs in quantum Ising chain in 1D, i.e, with σ x,z being the spin-1/2 Pauli matrices, has been studied in Refs.[27,30], where it has been shown that the SRE densities peak at the critical point h c = 1, and follow universal critical finite-size scaling hypothesis.In Fig. 4, we show the results for Rényi-2 SRE M 2 , estimated efficiently using the subleading term c L = 2M 2 (L/2) − M 2 (L) as described in Sec.III.Surprisingly, the sampling errors of the SRE density m 2 scales slower than log L, with L being the system-size, even at the critical point h c = 1.Therefore, unlike the MPS-based 2n-replica method employed in Ref. [27] that suffers from a computational cost of O(χ 12 ), our Monte Carlo method for estimating m 2 provides accurate results without being severely limited by χ.Moreover, the computation of m 2 using the perfect sampling of MPS [28,29] will necessarily incur statistical errors that are exponential in system-size as the direct estimation of the subleading term c L = 2M 2 (L/2) − M 2 (L) is not feasible by perfect sampling.
In the following, we extend the studies of SREs in 1D quantum many-body systems to qutrit systems by considering the three-state Clock model and the spin-1 XXZ model in 1D.

Three-state Clock model
The quantum Clock model is a generalization of the quantum Ising model with d states per site.Here we focus on the case d = 3, where the Hamiltonian is given by where X, Z are the shift and clock operators in Eq. ( 26) with d = 3.The model is equivalent to the three-state Potts model [89].There is a transition from the ferromagnetic phase to the paramagnetic phase at h c = 1, as in the quantum Ising model.The critical point is described by Z 3 parafermion CFT, with central charge c = 4/5.The exact correlation length exponent is ν Potts = 5/6 [89].It is to be noted that, since the system obeys Z 3 symmetry, a two-site update scheme (see Sec. III B) is required to sample the Pauli strings that preserve the symmetry.Indeed, the Pauli strings that preserve the Z 3 symmetry, generated by i Z i , are generated by Z i and X † i X j (up to a phase constant).
In the three-state Clock model, the magic density displays similar behavior as in the quantum Ising model [27,30], as shown in Fig. 5(a).Namely, m 1 displays maximum at the critical point h c = 1.We further investigate the finite-size scaling of m 1 , that has been done for the quantum Ising chain [27], using the finite-size scaling hypothesis: where m 1,m is the maximum SRE density at h c = 1.In Fig. 5(b), we show the data collapse corresponding to the finite-size scaling relation of Eq. ( 27), where we obtain the critical exponent ν ≈ 0.844, close to the expected theoretical value ν Potts = 5/6.We extract the critical exponent ν ≈ 0.844 and γ ≈ 0.66.The correlation-length exponent ν is close to the known νP otts = 5/6.We used bond dimension up to χ = 36 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.

Spin-1 XXZ chain
Next, we consider a spin-1 XXZ chain with single-ion anisotropy, whose Hamiltonian reads where S α 's, α = x, y, z, are the spin-1 operators, ∆ is the easy-axis anisotropy, and D is the single-ion anisotropy.The model has a global U (1) symmetry corresponding to the conservation of total magnetization i S z i , and here we consider the scenario of zero total magnetization.
The phase diagram of the model has been studied in previous works [90][91][92][93].For ∆ > 0, the model hosts three phases (with increasing D): the antiferromagnetic Néel order, the symmetry-protected topological (SPT) Haldane phase, and the large-D trivial phase.The Néel to Haldane transition is an Ising transition, while the Haldane to large-D transition is a Gaussian transition.
Here, we focus on the isotropic case, i.e., ∆ = 1.In this case, the transition is known to be at D ∼ −0.3 and D ∼ 0.97 for Néel-Haldane and Haldane-large D transitions, respectively [91][92][93].Fig. 6(a) shows the SRE density m 1 .We observe that m 1 is large and rather constant in the topological Haldane phase, while it becomes smaller in the neighboring phases.Note that the maximum value of m 1 for a product state is 2 3 log(4) ≈ 0.92, achieved by the tensor product of single-qutrit states, each of which has ⟨P ⟩ 2 = 1/4 for all P ̸ = I.Thus, it is seen that the magic in the SPT Haldane phase almost saturates the maximum value.

Long-range SRE
In the spin-1 XXZ chain, while the onset of the topological Haldane phase is rather apparent from the magic density, there is no clear peak at the transitions, rendering the determination of the critical point difficult.
The long-range magic, for the the spin-1 XXZ chain, as plotted in Fig. 6(b) shows clear extremums at the two transitions.Although L(ρ AB ) is still non-zero for small L away from criticality, it quickly decays to zero as the system size is increased.The peak at the Gaussian transition is very close to D ∼ 0.97, as obtained with DMRG up to L = 20000 spins [92].Notably, our results are obtained with only moderate sizes, and without any prior knowledge of the order parameter.At the Ising transition, the extremum occurs at a negative value as a minimum.Unlike entanglement, the SRE is not known to satisfy subadditivity, meaning that it is not always the case that L(ρ AB ) ≥ 0. Nevertheless, the non-trivial value at criticality is a useful indicator for detecting criticality.
The decay of long-range SRE away from criticality can be understood through a simple physical argument.Within a gapped phase characterized by a finite correlation length, when considering two subsystems A and B separated by a distance exceeding the correlation length, A and B are approximately uncorrelated.More formally, ρ AB ≈ ρ A ⊗ ρ B , which implies L(ρ AB ) ≈ 0. In contrast, at criticality, the correlation length becomes infinite, such that A and B are always correlated regardless of their distance.This results in a non-trivial value of L(ρ AB ).
We also come back to the quantum Ising chain (Eq.( 25)), and investigate the long-range magic across the Ising transition.We observe that L(ρ AB ) peaks at the critical point, as shown in Fig. 7. Furthermore, we plot L(ρ AB ) at h c = 1 in the inset of Fig. 7, where we see that the long-range magic grows logarithmically in L. In contrast, L(ρ AB ) quickly decays away from criticality (not shown).We note that, at the critical point, we observe long autocorrelation times between samples, which is the reason for the growing errors for larger sizes.This is reminiscent of the problem of critical slowing-down in the Monte Carlo simulations at criticality [94].It is thus interesting to develop a cluster update, akin to Wolff cluster update [95], that may overcome this issue, which we leave for future studies.
Based on the favourable scaling of our scheme with system size, we investigate the non-stabilizerness in 2D systems, which so far have not been properly explored in the literature.In particular, we consider a Z 2 lattice gauge theory, with Hamiltonian: where the spin-1/2 Pauli operators, τ α , α = x, z, live on the links of the square lattice.The first term is the plaquette term that flips the four spins on an elementary square plaquette of the lattice.We are interested in the charge-free sector, that satisfies the Gauss' law on each vertices of the lattice.It is well known that the Hamiltonian in Eq. ( 29) is dual to the 2D transverse-field Ising model on the square lattice by Wegner duality [96].Here, the spin-1/2 Pauli operators, σ α , α = x, z, live on the lattice sites of the dual square lattice.It can be shown that the duality transformation preserves SREs (see Appendix C).This enables us to compute the SREs of the the Z 2 gauge theory ( 29) by considering the ground state of the transverse-field Ising model, which is computationally more convenient for TTNs.At the same time, our results also shed light on the transition point of the Ising model: there.the transition from ferromagnetic phase to the paramagnetic phase is known to be at h c ≃ 3.04, as obtained with Quantum Monte Carlo [97].In the lattice gauge theory framework, such transition corresponds to confined to deconfined transition, where the behavior of Wilson loops turns from area to perimeter law.The results for magic density for n = 1, 2 are presented in Fig. 8.It is seen that both quantities detect the transition.However, the observed behavior is very different from the 1D quantum Ising chain, which exhibits a peak at the transition.Instead, here we observe that the curves exhibit crossings at the transition.
In Fig. 9(a), we depict m 1 close to the critical point, using a fixed bond dimension χ = 30.Remarkably, we observe that m 1 detects the transition point very well: all the curves cross near the critical point at h c = 3.04 (1).We should highlight at this point that the TTN ansatz with such a low bond dimension of χ = 30 can not approximate the ground state wave function accurately near the critical point, particularly in 2D critical systems.Consequently, the standard phase transition detectors, such as the Binder cumulant, calculated from the TTN state with χ = 30, not exhibit the expected critical crossing behavior -see Appendix D for a direct comparison in the present case.Therefore, the remarkable observation of the perfectly crossing behavior in m 1 near the critical point underscores the significant value of magic in detecting and characterizing quantum phase transitions.This is particularly relevant in situations where other quantities are prone to significant errors, e.g., due to limited bond dimensions in tensor network states.While we believe that a further characterization of what the scaling resources (e.g., size and bond dimension) to detect a transition point are is outside the scope of our paper, this would be very much worth pursuing based on the Ising model results we presented.
Furthermore, we show excellent data collapse for m 1 in Fig. 9(b), using the finite-size scaling relation of Eq. ( 27), from which we extract the correlation length exponent ν = 0.64 ± 0.05, that is close to the known ν 3D = 0.63 for 3D (classical) Ising universality [97].

V. EXPERIMENTAL PROTOCOL
The numerical method described above can be easily adapted for experimental measurements of SREs.In particular, we can sample Pauli strings according to Ξ P using Monte Carlo sampling.We note that, although the probability distribution Ξ P can be sampled directly through measurements in the Bell basis [60,98], the method requires preparation of two copies of a state and joint operations on them.In practice, this may not be feasible in some experimental platforms, or difficult to scale up to larger sizes and higher-dimensional systems.Moreover, the method only works for real wavefunctions [99].Instead, our proposal relies solely on measurements in the computational basis on a single instance of a state, and it is applicable to generic quantum states.
In experiments, the Pauli strings are measured from N M copies of ρ where the measurement outcomes are A i ∈ {+1, −1}.The expectation value is then given by the average taken over the random measurement outcomes.The sampling of Pauli strings can be performed with Metropolis algorithm, similar to our numerical calculations.However, it is important to note that in experimental setups, the candidate Pauli string is not restricted to few-site updates, as is the case of TTN.This flexibility allows for multi-site updates and can potentially reduce the autocorrelation time associated with the sampling process enormously.
For a finite number of measurements N M , we have that is an estimate for ⟨P ⟩.The total number of resources is thus N M × N S , where N S is the number of sampled Pauli strings.In view of Eq. ( 16), when the SREs are at most O(log N ), the required N S is polynomial in N .Note that N M may still be exponential, but it is expected to be no larger than O(d N ), with d being the local dimension.As a result, the number of resource required in our protocol is significantly lower than the protocol in Ref. [23] when the SREs are at most O(log N ).Moreover, our protocol offers a possibility to measure M 1 , in which case N S is always polynomial 3 .
The variance of the estimator in Eq. ( 32) is given by Var(P ) = 1 − ⟨P ⟩ 2 .Thus, the standard error reads For large N M , the random variable P approximately has a Gaussian distribution with average ⟨P ⟩ and standard deviation ∆P .Note that this will introduce bias to the estimators in Eq. ( 12) and Eq. ( 13).This bias can be made smaller by increasing N M , where the estimators become unbiased in the limit N M → ∞.
Here, we simulate this situation numerically by perturbing the computed ⟨P ⟩ with ϵ, where ϵ is a random number chosen from a Gaussian distribution centered at zero and with standard deviation ∆P .We would like to investigate the effects of taking finite N M and N S .Here, we consider the ground state of 1D transverse-field Ising chain at h = 1 for concreteness.An example of the results of such a protocol is shown in Fig. 10 for L = 8 with N M = 500 and N S = 10000.11.We see that, for fixed N S , the error first increases for small N M , before it eventually decreases.We expect this is due to the bias with finite number of Pauli measurements, as mentioned above.Indeed, as shown in Fig. 11(b), we see that increasing N S while fixing N M does not result in vanishing δm n .

VI. CONCLUSIONS AND OUTLOOK
We have proposed a Markov chain Monte Carlo approach to compute magic in many-body systems.We have discussed how the full state magic M n can be estimated for different values of n, and demonstrated the corresponding efficiency in several scenarios.Moreover, long-range magic can be estimated efficiently in general.The implementation of our algorithm is flexible and compatible with various wave-function based methods.Specifically, we have provided detailed insights into the efficiency and flexibility of our method when applied to tree tensor networks.
Through our algorithm's flexibility, we have gained valuable insights into the role of magic in many-body systems.In one-dimensional systems, we observed that full state magic is not universally associated with critical behavior.While it displays criticality signatures in certain cases like Ising and Potts models, it does not in others.However, long-range magic overcomes this limitation and consistently exhibits indications of critical behavior across all scenarios we investigated.We speculate that the functional form of long-range magic, similar to mutual information, is free of potential UV-divergences in a field theory framework.
The very mild volume scaling cost of our sampling has also enabled us the exploration of two-dimensional Z 2 lattice gauge theories.There, we have found that magic displays finite-volume crossings in correspondence of the confined-deconfined phase transition, and it also follows universal scaling behavior up to the volumes (100 spins) we were able to treat.Remarkably, magic was well converged even at modest bond dimensions.
Our numerical results suggest a deep connection between (long-range) magic and many-body properties, highlighting the direct links between stabilizer Renyi entropies and physical phenomena such as quantum critical behavior and confinement-deconfinement transitions.To complement our theoretical findings, we have proposed an experimental protocol for measuring stabilizer Renyi entropies solely using measurements in the computational basis.
In terms of future investigations, our technique can be extended to explore nonstabilizerness in finitetemperature scenarios by generalizing it to tree-tensor operators that efficiently represent low-temperature many-body states.In particular, it would be interesting to study the behavior of stabilizer Renyi entropies at finite-temperature phase-transition and compare it with other information-theoretic quantities, such as entanglement [85,[100][101][102], quantum discord [103], and quantum coherence [104].Along the same lines, another possible scenario would be applying our tools to faulty quantum circuits, recently discussed in the context of magic in Ref. [105].It would also be instructive to perform a systematic investigation of magic within topological phases, extending our analysis of the Haldane phase.Another interesting perspective is to understand the role of magic in many-body quantum dynamics of closed quantum systems, whose investigation in the context of Ising models has been the subject of recent works [106].In particular, our method allows for the investigation of genuine long-distance magic, that might be instrumental in establishing the presence or absence of propagation bounds for magic.
At the methodological level, our work opens a series of questions.The Markov chain Monte carlo approach could be extended to investigate other magic measures that depend only on expectation values, such as mana.Moreover, so far, we have only employed very basic sampling strategies.It would be worth exploring how different ones, such as heatbath or non-local updates, can be used to design better magic estimators since, in terms of experimental applicability, having shorter autocorrelations could considerably improve realistic implementations.In terms of efficiency of the increment trick in 2D models, it would be interesting to study whether a one-dimensional projection of 2D systems such as the one introduced in [107] would be beneficial.Finally, it would be interesting to understand the finer structure of sampling Pauli strings in many-body systems, that could reveal both useful insights into novel algorithms, and potentially deeper connections between many-body properties and magic.The duality transformation between Eq. ( 29) and Eq. ( 31) is defined with the following transformation, More precisely, the transformation maps the charge-free sector of Eq. ( 29) to the even sector of Eq. (31).
It is easy to see that the mapping in Eq. (C1) maps Pauli strings in the Ising model to Pauli strings in Z 2 gauge theory, because the Pauli operators on both sides of the equation generate the Pauli group in the corresponding models.Since the SREs depend only on the expectation values of Pauli strings, it follows that the SREs are preserved by the duality transformation.Therefore, the SREs in the Ising model are identical to the SREs in Z 2 gauge theory.It should be, however, noted that equivalence relation in case of the subsystem mixed-state SRE (e.g., M2 defined in Eq. ( 5)), and the long-range magic thereof, is non-trivial because of the non-local nature of the transformation (C1).Consequently, the distribution of magic within the subsystems may differ in these two theories.
It is worth nothing that the same conclusion evidently holds for other dualities that map Pauli strings to Pauli strings, such as the Kramers-Wannier duality which maps h → h −1 in Eq. ( 25) and Eq. ( 26).As previously discussed, the long-range magic is not preserved under the duality.This is reflected in the distinct behavior of L(ρ AB ) for h > 1 and h < 1 in Fig. 7.In Sec.IV B, we have demonstrated the ability of the magic density to accurately detect and characterize the quantum critical point in the 2D Z 2 gauge theory, and thereby in 2D quantum Ising model.Notably, the curves of m 1 for different linear system-sizes exhibit a clear critical crossing behavior near the critical point h c = 3.04, even with a modest TTN bond dimension of χ = 30.However, the same level of accuracy is not achieved when utilizing the Binder cumulant, defined as for the 2D Ising model (31).Due to the inability of the TTN state with a small bond dimension of χ = 30 to faithfully represent the ground state in the vicinity of the critical point, the calculation of the Binder cumulant U yields erroneous results.Consequently, the curves of U for different linear system-sizes L do not exhibit a clear crossing behavior near the critical point (Fig. 14).For instance, while the curves for L = 4 and 5 intersect at h = 2.98, the intersection for L = 7 and 8 occurs around h = 3.14.As such, if one attempts to perform finite size scaling on the Binder cumulant data, the resulting critical point and the correlation-length critical exponent ν will be erroneous.

FIG. 2 .
FIG. 2. Schematics of partitions.(a) Full partition.(b) Two widely-separated partitions for the calculation of longrange magic in Eq. (7).(c) Subleading term as in Eq. (24), as well as a cartoon depicting the increment trick discussed in the main text.

FIG. 3 .
FIG. 3. Efficient Monte Carlo sampling using Tree Tensor Network.(a) Tree Tensor Network (TTN) representation of a many-body wavefunction |ψ⟩, where tensors are depicted as circles arranged in a binary-tree structure.Each tensor is identified by a pair of zero-indexed integers [l, n], representing its layer index l and tensor index n at that layer.The red circle at the top-most layer represents the root tensor having index [0, 0], where the isometry center of the TTN is taken.(b) To evaluate the expectation value of a tensor-product of single-site operators ⟨O1O2 . . .ON ⟩, we first place each operator Oi at the physical site it acts on in the TTN representation.Then, we compute the effective link operators which live at the virtual links by the coarse-graining procedure as shown in the figure.The coarse-graining is performed iteratively from the physical sites to the top-most virtual links, which are directly connected to the root tensor.At each step, the link operators O [l+1,2n] and O [l+1,2n+1] are combined into O [l,n] by the [l, n]-tensor.The resulting link operator O [l,n] acts on the [l − 1, ⌊n/2⌋]-tensor one layer above in the TTN structure.(c) The expectation value ⟨O1O2 . . .ON ⟩ is calculated from the contraction of the root [0, 0]-tensor and the top-lost link operators as shown in the figure.(d) Considering a modified operator which differs only at a single site from the previous one, O1O2...O ′ i ...ON , we only need to recompute the link operators in the path from the modified physical site i to the topmost link.

FIG. 4 .
FIG.4.Efficient estimation of Rényi-2 SRE density in 1D quantum Ising chain.(a) The subleading term for the Rényi-2 SRE, cL = 2M2(L/2) − M2(L), directly estimated using the efficient scheme specified in Sec.III, for various system-sizes in 1D quantum Ising chain.(b) The SRE density m2 for the 1D quantum Ising chain near the critical point hc = 1 computed using the increment method using different subleading terms.(Inset) The sampling errors for m2 at hc = 1 for various system-sizes L (in log-scale).Clearly, the errors show even slower than than logarithmic growth for the efficient sampling scheme.Here we consider TTN bond dimension χ = 30 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.

96 FIG. 5 .
FIG. 5.Magic density in 1D quantum three-state Clock model.(a) The SRE density m1 in the ground-state of the three-state Clock model as a function of h.(b) Finitesize scaling for m1.Here.m1,m is the maximum m1 at hc = 1.We extract the critical exponent ν ≈ 0.844 and γ ≈ 0.66.The correlation-length exponent ν is close to the known νP otts = 5/6.We used bond dimension up to χ = 36 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.

64 FIG. 6 .
FIG. 6. Magic density and long-range magic in spin-1 XXZ chain.(a) the magic density m1 and (b) long-range magic L(ρAB) of the ground-state of the spin-1 XXZ model with ∆ = 1 as a function of D. We consider bond dimension up to χ = 60 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.The dashed vertical lines represent the best estiamtes available for the transition points.

2 FIG. 7 .
FIG.7.Long-range magic in 1D quantum Ising chain.The long-range magic L(ρAB) as in Eq. (7) in the ground state of 1D quantum Ising chain as a function of the transverse field h.It peaks at the critical point hc = 1.(Inset) L(ρAB) at hc = 1 for various system sizes L (in log-scale).We consider TTN bond dimension up to χ = 30 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.

10 FIG. 8 . 8 FIG. 9 .
FIG. 8. Magic densities in 2D Z2 gauge theory.The SRE desnities (a) m1 and (b) m2 of the ground-state of Z2 gauge theory on L × L square lattice as a function of h.We use TTN bond dimension up to χ = 60 and the number of sample is NS = 10 6 .Error bars represent 95% confidence interval.

8 FIG. 14 .
FIG. 14.The Binder cumulant across the critical point in the 2D quantum Ising model.Here we approximate the ground state of 2D quantum Ising model with TTN having bond dimension χ = 30, in parity with Fig.9.
Appendix D: Binder cumulant in 2D quantum Ising model with TTN