Magic in generalized Rokhsar-Kivelson wavefunctions

Magic is a property of a quantum state that characterizes its deviation from a stabilizer state, serving as a useful resource for achieving universal quantum computation e.g., within schemes that use Clifford operations. To date little is known about properties and behaviour of magic in many body quantum systems. In this work, we study magic, as quantified by the stabilizer Renyi entropy, in a class of models known as generalized Rokhsar-Kivelson systems, i.e., Hamiltonians that allow a stochastic matrix form (SMF) decomposition. The ground state wavefunctions of these systems can be written explicitly throughout their phase diagram, and their properties can be related to associated classical statistical mechanics problems, thereby allowing powerful analytical and numerical approaches that are not usually available in conventional quantum many body settings. As a result, we are able to express the SRE of integer Renyi index n > 1 in terms of wave function coefficients that can be understood as a free energy difference of related classical problems. We apply this insight to a range of quantum many body SMF Hamiltonians, which affords us to study numerically the SRE of large high-dimensional systems, unattainable with existing tensor network-based techniques, and in some cases to obtain analytical results. We observe that the behaviour of the SRE is relatively featureless across quantum phase transitions in these systems, although it is indeed singular (in its first or higher order derivative, depending on the first or higher order nature of the transition). On the contrary, we find that the maximum of the SRE generically occurs at a cusp away from the quantum critical point, where the derivative suddenly changes sign. Furthermore, we compare the SRE and the logarithm of overlaps with specific stabilizer states, asymptotically realised in the ground state phase diagrams of these systems. We find that they display strikingly similar behaviors, which in turn establish rigorous bounds on the min-relative entropy of magic.


Introduction
Stabilizer states are an important class of quantum states in quantum information theory [1,2].They have very rich structures [3,4] and can be highly entangled [5,6,7,8].However, it is wellknown by the Gottesman-Knill theorem (and its subsequent extensions) that quantum computation using only stabilizer states and Clifford circuits can be efficiently simulated on a classical computer [1,9,10,11].In quantum computation using the state injection scheme [12,13,14], universal quantum computation is achieved by injection of states outside the stabilizer set, known as magic states, while keeping the set of operations restricted to Clifford operations.Magic thus serves as a fundamental resource that would be required to outperform classical simulations.To quantify the amount of magic resource in a quantum state, the notion of a magic measure has been introduced within the framework of resource theory [15].These measures assess the amount of resource a state can provide in quantum computation by state injection scheme, offering insights into the computational power and quantum capabilities of different states.Most measures that have been introduced require optimization procedures to compute them (see, e.g., Refs.[16,17,18,19,20]), and are thus diffi-cult to evaluate beyond a few qubits [21].More recently, a computable measure of magic has been introduced, called Stabilizer Renyi Entropy (SRE) [22], which is relatively easier to compute as it is expressed only in terms of expectation values of Pauli strings.Nonetheless, examples of analytical results for the SRE remain few and far between, and computational cost often limits numerical studies to relatively small systems.
In recent years, there is an increasing interest in characterizing the role of magic in the ground states of quantum many-body systems [23,24,25,26,27,28,29,30,31,32,33,34,35] -namely, 'how far' such states are from being stabilizer states.Notable progress has also been made in its experimental measurements [36,37,38,31,39,40].The study of quantum information concepts in the realm of many-body theory has a long history, with a prominent example being entanglement, which has proven to be a powerful tool for investigating various many-body phenomena [41,42].The comprehensive exploration of many-body magic thus represents a promising avenue that holds significant relevance for characterizing the quantum complexity of different phases of matter.This includes understanding the computational cost of simulating specific phases and assessing their computational capabilities.In a more practical context, this would also provide insights into the potential of manybody states to be an input in magic state distillation [43], for applications in quantum information processing.
The investigation of many-body magic has been enabled by the recent development of numerical methods based on tensor networks to efficiently compute the SRE [28,29,30,31,33].In particular, these studies have suggested a connection between magic and criticality in onedimensional quantum systems.Notwithstanding these interesting developments, the cost of computing the SRE remains very high, often restricting the study to simple one-dimensional systems.Further studies of highly entangled states, such as higher-dimensional systems, appear to be practically out of reach with current methods (see however Ref. [31]).
In this work, we introduce an approach to compute the SRE with integer Renyi index n > 1 in many-body wavefunction, by expressing it in terms of wavefunction coefficients that make it amenable to computation using Monte Carlo sampling (provided the wavefunction can be gauged to have non-negative coefficients).
We apply this approach to a class of models known as generalized Rokhsar-Kivelson systems [44,45], or Hamiltonians that allow a stochastic matrix form (SMF) decomposition [46].The ground state wavefunctions of these systems can be written explicitly throughout their phase diagram, and their properties can be related to associated classical statistical mechanics problems in thermodynamic equilibrium at temperature T , which plays the role of a parameter in the phase diagram of the original quantum systems.This correspondence allows powerful analytical and numerical approaches to be deployed, that are not usually available in conventional quantum many body settings 1 .
Since the early work of Rokhsar and Kivelson [48], SMF Hamiltonian and wavefunctions have appeared in many different physics contexts.Of late, a resurgence of attention has derived from the fact that they can be naturally realised using tensor networks and PEPS [49], and they can be implemented in measurementprepared quantum states and (monitored) quantum circuits [50,51,52,53,54,55].In this context, the magic of SMF wavefunctions thus directly quantifies the amount of non-Clifford resources required to prepare these systems in the circuits.
We are able to express the SRE of SMF wavefunctions in terms of a free energy difference of related classical problems, which can then be efficiently computed by thermodynamic integration.We apply this insight to a range of quantum many body SMF Hamiltonians, which affords us to study numerically the SRE of large high-dimensional systems, unattainable using existing tensor network-based techniques, and in some cases we obtain explicit analytical results.
We observe that the behaviour of the SRE is relatively featureless across quantum phase transitions in these systems [56], although it is indeed singular in its first or higher order derivative, depending on the first or higher order nature of the transition.On the contrary, we find that the maximum of the SRE generically occurs at a cusp away from the quantum critical point, where the derivative suddenly changes sign.Furthermore, we compare the SRE to the logarithm of overlaps with specific stabilizer states, that are asymptotically realised in the ground state phase diagrams of these systems.We find that they display strikingly similar behaviors, which in turn establish rigorous bounds on the min-relative entropy of magic.
The rest of the paper is structured as follows.In Sec. 2 and Sec. 3 we give a brief review of SMF Hamiltonians and SRE, respectively.We then state our general results in Sec. 4 about the SRE and its upper bounds in SMF systems.In Sec. 5 we then present a study of a range of models, encompassing the Ising ferromagnet in 1D, 2D, 3D, and infinite dimensions; the J 1 − J 2 model on the square lattice; the triangular Ising antiferromagnet; and the Edwards-Anderson model on the cubic lattice.Finally, we conclude in Sec. 6.

Brief review of Stochastic Matrix Form (SMF) Hamiltonians
The stochastic matrix form wavefunctions, dependent on the parameter β, are given by [44,45,46] where Z(β) can be seen as a classical partition function at temperature T = 1/β.One can design a quantum Hamiltonian for arbitrary choice of E σ , such that |ψ SMF ⟩ is the ground state of the Hamiltonian.In particular, for a locally interacting E σ , the Hamiltonian also contains only local interactions.The quantum Hamiltonian is said to be SMF decomposable [46].The equal-time correlation function of diagonal operators of SMF wavefunctions are given by the equal-time correlations functions in the associated classical systems in thermal equilibrium.It follows that the ground state phase diagram of the quantum Hamiltonian contains the thermal phase diagram of the classical system in thermal equilibrium.Since the wave function coefficients are known exactly by design, the wave function can be sampled with classical Monte Carlo simulations of the corresponding classical system.

Stabilizer Rényi entropy
Stabilizer Rényi Entropies (SREs) are a measure of nonstabilizerness recently introduced in Ref. [22].For a pure quantum state |ψ⟩ of a system of N qubits (equivalently, spin-1/2 degrees of freedom), SREs are expressed in terms of the expectation values of all Pauli strings in P N : Eq.
(3) can be seen as the Rényi-n entropy of the classical probability distribution Ξ P (|ψ⟩) = ⟨ψ|P |ψ⟩ 2 /2 N , also known as the characteristic function [57].The SREs have the following properties [22] where F STAB is the stabilizer fidelity defined as The following inequality holds [29] M In particular, setting n = 2, and defining for any |ϕ⟩ ∈ STAB.Because the SRE is defined in terms of the characteristic function Ξ P (|ψ⟩), it is natural to estimate it through statistical sampling of Ξ P (|ψ⟩).Indeed, this was the core of previous numerical methods that have been introduced to compute the SRE [29,30,31], which are so far limited to tensor network techniques.To obtain reliable statistics, a large number of samples is required, often resulting in computations being restricted to small bond dimension.This limitation poses a challenge when studying highly entangled systems, such as higher dimensional ones.Unfortunately, the existing methods are not directly suitable for Quantum Monte Carlo approaches due to the inherent difficulty in evaluating the expectation values of high-weight Pauli strings within.

Magic in SMF ground states
In this section we show how the special form of the ground state wave functions of SMF systems allows for analytical and numerical routes into the calculation of their magic, that are not afforded to general quantum many body systems.In doing so, we develop the machinery that will later be used to study a broad range of model systems, to gain insight in the behaviour of this intriguing quantity.
For the SMF systems introduced in Sec. 2, and in particular their ground state wavefunctions, these known expressions for the magic can be further manipulated upon substituting c σ = e −βEσ/2 / √ Z in Eq. ( 8) to obtain: where 1) , σ (2) , σ (3) , σ (4)   exp −β One can interpret Z M as a classical partition function constructed from four copies of the original classical degrees of freedom.The second term in the square bracket describes the energy of a configuration obtained from the spin product of three out of four copies.The expression in Eq. ( 9) can thus be seen as (proportional to) the difference between the free energy of the classical system described by Z M and four non-interacting copies of the original classical system described by Z.This manipulation is not only interesting from a conceptual point of view, but also from a pragmatic one: it provides a new angle to compute the magic of a quantum (SMF) state using classical statistical mechanics tools in the same number of dimensions.As we demonstrate later, it affords us the ability to access significantly larger system sizes and higher dimensional lattices than previously possible [29,30,31].
In practice, rather than computing ratios of partition functions or differences in free energies, it is convenient to notice that the derivative of M 2 with respect to temperature reduces to: Thus, M 2 can be more efficiently obtained by computing the energies ⟨E⟩ and ⟨E M ⟩ M and then proceding to integrate the r.h.s. of Eq. (11).

Upper bound of M 2
As discussed in Sec. 3, M 2 is bounded from above by D(|ψ⟩, |ϕ⟩) = − log |⟨ϕ|ψ⟩| 2 for any |ϕ⟩ ∈ STAB.In the following sections, we compute D(|ψ⟩, |ϕ⟩) for several ad hoc choices of states |ϕ⟩, specific to the system being considered.Once again, in the case of SMF wavefunctions, these overlaps can be reduced to classical statistical mechanical objects, amenable to corresponding analytical or numerical estimates.
Here, we illustrate the procedure to compute D(|ψ⟩, |ϕ⟩) in a couple of cases that will often be used in the following.Consider for example |ϕ⟩ = | + + + ...⟩, where |+⟩ = |↑⟩+|↓⟩ 2 (i.e., a spin-1/2 state polarized in the x direction).In the context of SMF wavefunctions, this is the ground state at β = 0. We denote = 1 we then obtain In the ferromagnetic Ising model, this is the ground state at T = 0 (namely, the exact ground state, ignoring spontaneous symmetry breaking effects).We denote where E zz is the energy of the configuration σ i = 1(−1), for all i, we obtain In all cases, the overlaps D(|ψ⟩, |ϕ⟩) are related to the partition function of the classical model.Therefore, similarly to M 2 , we compute them using direct thermodynamics integration.

SMF models
Armed with the tools developed above, we proceed to investigate a broad range of models, in the attempt to deepen our understanding of stabilizerness and magic in many body systemsalbeit of the fine tuned SMF kind -and its relation to quantum phase transitions.For this purpose, we consider in the first instance quantum Ising ferromagnets in 1D, 2D, 3D, and infinite dimensions (mean field); we also consider the J 1 -J 2 model tuned to exhibit a first order phase transition, for comparison.We then move on to more exotic examples, such as the quantum triangular Ising antiferromagnet (which is fully frustrated in the SMF realisation) and the Edwards-Anderson model (exemplifying a spin glass transition in the droplet picture).
Since the SREs are generally an extensive quantity, we focus on the SRE densities M n /N , where N = L d is the total number of sites in the system with linear size L and dimensionality d.For the numerical simulations at finite volume, we impose periodic boundary conditions.

1D SMF Ising ferromagnet
The quantum SMF Hamiltonian for the 1D Ising ferromagnet is related to the classical 1D Ising model with energy where the spins σ i live on a chain, and it reflects its thermodynamic behaviour.In particular, there is no phase transition and the system orders only in the limit β → ∞ (we stress that β plays the role of inverse temperature for the classical system whereas it is merely a tunable parameter in the SMF quantum Hamiltonian, which is considered to be at zero temperature in its ground state, for the purpose of this work).Despite its simplicity, this model serves as a useful warmup example and the magic can be computed analytically using Eq.(9).Indeed, we recall that the partition function of the 1D (nearest-neighbour) Ising model can be computed with transfer matrix techniques: where V η,η ′ = e βηη ′ is a 2 × 2 matrix.The eigenvalues of V are λ 1 = 2 cosh β and λ 2 = 2 sinh β.Thus, We can similarly compute Z M by interpreting the four layers of the 1D chain as a 1D chain with 16 states per site: , where V M is a 16 × 16 matrix.
More generally, for an integer index n > 1, the transfer matrix V M,n is a 2 2n × 2 2n matrix with elements To compute Z M,n = Tr(V L M,n ) we work in the thermodynamic limit L → ∞, where we only need to find the largest eigenvalue of V M,n .One can verify that the column vector with all elements equal to 1 is an eigenvector of V M,n .By the Perron-Frobenius theorem, the corresponding eigenvalue is the unique largest eigenvalue: Finally, using Eq. ( 22) and (20), we find Furthermore, D x and D zz can be computed directly by plugging in the partition function Eq. (20) into Eq.( 15) and (17), respectively.In Fig. 1, we show the SREs, D x , and D zz of the SMF 1D Ising model, observing the expected asymptotic agreement in the limits T → 0 and T → ∞.
To avoid confusion, we remark here that the SMF Hamiltonian and corresponding GS wavefunction related to the classical 1D Ising ferromagnet are strikingly different from the conventional 1D Ising ferromagnet in a transverse field.Most notably, in the SMF case there is no phase transition and ordering occurs only asymptotically in the limit T → 0.

2D SMF Ising ferromagnet
The quantum SMF Hamiltonian of the 2D Ising ferromagnet is related to the 2D classical Ising model with energy where the spins σ i are taken without loss of generality to live on the 2D square lattice.There is a phase transition between the ferromagnetic and the paramagnetic phase at To study the thermodynamics properties of Z M , we perform Monte Carlo simulations augmented with Wolff cluster algorithm [58] and parallel tempering [59,60].We study the energy ⟨E M ⟩ M as a function of temperature for a range of system sizes as shown in Fig. 2 (top panel).We observe a behaviour compatible with a first-order transition from a high-temperature paramagnetic phase to a low-temperature ordered phase, at some value T * ̸ = T c (whereas we know ⟨E⟩ from the classical 2D Ising model to be smooth, with a singularity in its first derivative at T c ).The presence of a first order transition in the classical system described by Z M will be a common feature in all examples considered in our work; however, while T * > T c for the 2D Ising case, we will find for example that T * < T c in higherdimensional systems.
Integrating the energy as in Eq. (11) from high temperature, we obtain the SRE M 2 .In contrast to the results of previous studies [28,31], we see that M 2 is continuous and does not exhibit a maximum nor minimum at the transition point.In fact, we know that at T c , M 2 inherits a singularity in its second derivative from the singularity in the first derivative of the energy ⟨E⟩ of the associated classical system described by the partition function Z (which undergoes a second order phase transition).The maximum of M 2 occurs instead away from the quantum critical point (into the paramagnetic phase), at a cusp that can be related to the first order transition point of the classical system described by Z M .Furthermore, we observe that the bounds 4D x and 4D zz lie very close above M 2 .By Eq. (7), this also establishes strict upper and lower bounds for D min .This is again a common feature that we consistently observe in all examples considered in this work.
Although the SRE M 2 appears relatively featureless across T c , we know from Eq. (11) that it must inherit any singularity present in ⟨E⟩ and in ⟨E M ⟩ M .The well-known critical behaviour of the 2D Ising model gives a singularity in the second derivative of M 2 , which is related to the specific heat of the classical 2D Ising model: d 2 M 2 /dT 2 exhibits a peak at T c which diverges logarithmically, as shown in Fig. 3.A (negative) peak is observed at T * , due to the the specific heat of Z M ; here the transition is first order and the peak diverges as N , much faster than the known logarithmic divergence at T c .
We note that, by Wegner duality [61], the SMF groundstates corresponding to the 2D Ising model are dual to a wavefunction deformation of the toric code studied in Ref. [62].Therefore, the SREs of the two wavefunctions are identical (up to a constant shift), since the SREs are preserved by Wegner duality [31].

3D SMF Ising ferromagnet
The discussion of the 3D Ising ferromagnet goes along the same line as in 2D, with the spins σ i living without loss of generality on the 3D cubic lattice.The model is known to exhibit a second-order phase transition between the ferromagnetic and the paramagnetic phase.Through large-scale Monte Carlo simulations, the critical point was found to be at T c ≈ 4.5115 [63].
The results are shown in Fig. 4. The energy ⟨E M ⟩ M once again exhibits a first-order transition that induces a cusp in the behaviour of M 2 at T * , where it reaches its maximum value.At the quantum phase transition, M 2 is once again continuous, with a singularity in its second derivative.Differently from the 2D case, the cusp (and maximum) of M 2 occurs in the ferromagnetic phase instead of the paramagnetic one.Once again, the upper bounds 4D x and 4D zz lie very close to M 2 .

Infinite-range Ising model
For completeness, we consider the case of an infinite-range Ising ferromagnet, with classical energy which becomes again analytically tractable.Hereafter, we shall neglect the trivial constant energy shift in the last line.
The partition function of the infinite-range model can be evaluated by first recasting it to a Gaussian integral and then performing a saddlepoint approximation, which is exact in the thermodynamic limit L → ∞ (see e.g., Ref. [64]).Explicitly, the free energy is given by where the saddle point magnetisation is found by solving the self-consistency equation The system exhibits a second-order phase transition at T c = 1 between the ferromagnetic and the paramagnetic phase.The evaluation of βF M = − log Z M follows along the same lines.First, we write the partition function as Then, we make use of the identity √ N amx (30) to obtain where In the limit L → ∞, we can evaluate the above integral using the saddle-point approximation, such that the free energy is given by Taking partial derivative with respect to all m a and q b , we obtain the self-consistent equations ) and ) respectively, for a, b = 1, 2, 3, 4. The outer summations above are over η (1) , η (2) , η (3) , η (4) = ±1.The quantity m a corresponds to the magnetization of the a-th layer, while q b corresponds to The procedure outlined above can be straightforwardly generalized to higher (integer) index n > 2. If we assume, by symmetry, that the solution satisfies m 1 = m 2 = ... = m 2n = q 1 = ... = q 2n = m 2 , then the self-consistent equations sim- 2 The symmetry between ma and q b may not be immediately obvious, but it can be seen as follows: for each site i in the layer a, define spin s plify to m = cosh(βm) 2n−1 sinh(βm) + sinh(βm) 2n−1 cosh(βm) 1 + cosh(βm) 2n + sinh(βm) 2n (38) One can verify that Eq. (38) always admits m = 0 as a solution, which minimizes F M at high temperature; it also admits the solution m = 1, which minimizes F M at T = 0, and the transition between the two is first-order.
Furthermore, D x and D zz can be computed directly from Eq. ( 15) and (17), respectively, using the free energy in Eq. (27).We plot them along the SREs in Fig. 5.
Similarly to the 3D Ising ferromagnet, and contrary to the 2D case, the cusp (and maximum) of M 2 occurs well within the ferromagnetically ordered phase.In fact, the behavior of M 2 along with D x and D zz are very similar to the 3D case.Interestingly, the stabilizer bound to the magic is exactly met by D x for any T larger than the cusp value.This is because in this regime log Z M,n = 2n log 2, while log Z(β/2) = log 2 in Eq. (15), which implies M n = 2n n−1 D x .By Eq. ( 6), it follows that D min = n−1 2n M n for any T larger than the cusp value of an index n.For T ≥ 1, all M n and D x vanish, i.e., the states are asymptotically close to the stabilizer state | + + + ...⟩.

J1-J2 model and first order behaviour
Up to now we considered quantum SMF Hamiltonians that exhibit continuous phase transitions.
Here we investigate what happens at a first order quantum phase transition by looking at the SMF J 1 -J 2 Ising model on the square lattice, related to a classical model with energy where J 1 , J 2 > 0. The first term corresponds to a ferromagnetic Ising nearest-neighbour interaction, while the second term is an antiferromagnetic interaction across the diagonals of the square plaquettes.For an appropriate choice of the system parameters, e.g., when the ratio g = J 1 /J 2 = 0.55, the system exhibits a first-order transition between a high-temperature paramagnetic phase and a low-temperature stripe phase at T c ≈ 0.772 (setting J 1 = 1 as the reference energy scale) [65].In the stripe phase, the ground states are fourfold degenerate, and can be understood as two decoupled antiferromagnetic ground states.
We also observe a first order phase transition in the associated coupled layered system Z M , but at a higher temperature T * , well into the paramagnetic phase (as in the 2D Ising ferromagnet, and contrary to 3D and infinite-range).Therefore, in this system we expect two discontinuities in the first derivative of M 2 with respect to T : One at the quantum phase transition (T c ), where the slope is positive on both sides and approximately doubles across it; the other at T * , where the slope changes sign abruptly, giving rise once again to a maximum in M 2 , where a cusp occurs.
In Fig. 6 we also compare M 2 with the bounds provided by the paramagnetic (D x , asymptotically accurate for T → ∞) and stripe (D stripe , asymptotically accurate for T → 0) phases.The latter appears to be remarkably close for all T < T * .

Antiferromagnetic triangular Ising model
We now proceed to a more exotic model where the quantum SMF Hamiltonian is related to the classical antiferromagnetic triangular Ising model [66,67], with energy where the spins σ i live on the sites of a triangular lattice, i = 1, . . ., N .The model features an extensive ground state degeneracy with algebraically decaying correlations at T = 0, while it is disordered at all temperatures T ̸ = 0 [66].We first show that the classical system described by Z M also features an extensive ground state degeneracy.To this end, we provide a lower bound on the zero-point entropy by explicitly constructing an extensive set of states with lowest energy.To do so, let us divide the triangular lattice on each layer into three sublattices.Let us then fix the spins on two of the sublattices to 1 and −1, respectively, with the same choice for all four layers.One can then straightforwardly verify that the spins on the remaining sublattice on each layer can be chosen arbitrarily without affecting the energy E M of the system, and that the latter is indeed minimised.The number of such states is 2 4N/3 , which implies As one can straightforwardly think of other configurations that minimize the energy, this bound is not tight.We find that the corresponding Z M features a phase transition that appears to be first-order, as displayed in Fig. 7, albeit of a less strong nature than in the cases discussed previously.Once again, M 2 exhibits a cusp at the transition point T * of the classical system described by Z M .
Unlike in the other models considered so far, in this case the ground state at T = 0 is not a stabilizer state.We note that the configurations with the three sublattice structure given above is also known as the clock state, which arises as the quantum ground state of the quantum triangular Ising antiferromagnet at small magnetic field [68,69,70,71].While it is not the exact ground state for T = 0, it constitutes a significant part of the ground state.Thus, we compare M 2 with D clock (|ψ⟩) = D(|ψ⟩, |ϕ⟩) where |ϕ⟩ is the clock state defined above, which is a stabilizer state.D clock is obtained in a similar way as D zz (see Sec. 4.2).At T = 0, D clock is given by where S(0) ≃ 0.3383 N is the zero-point entropy of the antiferromagnetic triangular Ising model [66].On the other hand, M 2 is given by This is obtained by setting T = 0 in Eq. (9).In this case, the observation that 4D clock is strictly larger than M 2 can be attributed to the fact that the zero-point entropy S M (0) is strictly larger than 4N 3 ln 2. In turn, this is also a manifestation of the fact that the ground state of the classical system described by Z(T ) is not a stabilizer state for T = 0.

Edwards-Anderson model
Finally, we consider an example of a disordered system, where the quantum SMF Hamiltonian is related to the Edwards-Anderson (EA) model, with energy  where the spins σ i live on the 3D cubic lattice.
Here, the couplings J ij are independently drawn from a Gaussian distribution with zero mean and unit variance.The system is known to undergo a continuous transition from the high-temperature phase to the low-temperature spin glass phase at T c ≈ 0.95 [72].This model has a unique ground state for any realization of J ij (up to global spin flip).For system sizes up to L = 10, the exact ground states and their energy can be readily obtained using the McGroundstate server [73].
We show the energy ⟨E M ⟩ M , the specific heat C v,M and the magic M 2 in Fig. 8. Again, the maximum occurs well into the paramagnetic phase, at the transition point of the coupled layered system Z M .Unlike in the previous cases, where the first order behaviour of the Z M transition was self-evident because of the noticeable discontinuity in ⟨E M ⟩ M , the situation is less clear-cut here.While the maximum of the specific heat appears to grow slower than N , suggesting a second order transition, we are unable to identify a clear scaling of the specific heat within the accessible system sizes.We further compute equilibrium energy histograms at different temperatures around T * , shown in Fig. 9.The behaviour closely resembles a trade off between two different peaks, whose positions are approximately temperature-independent (although we are unable to see the minimum in between them scale to zero as a function of system size, within the systems accessed in this work).Overall, we suggest that the transition in this model is weakly first-order.
In Fig. 8 we also compare M 2 with bounds from D x and D GS , where D GS is obtained from the overlaps with the exact ground state for each realization, computed by thermodynamic integration as discussed in Sec.4.2, using the McGroundstate server [73] to obtain the exact ground state energy.In this case, the bounds are somewhat higher than encountered in previous cases.Nevertheless, the crossing between D x and D GS still occurs close to the maximum of M 2 .
Finally, we investigate the nature of the lowtemperature phase of Z M .A natural candidate is a spin glass phase, accompanied by replica symmetry breaking (RSB), akin to the lowtemperature phase of the 3D EA model.To detect RSB, we compute the spin overlap where α and β represents two copies of the system with the same disorder.We show the spin overlap in the coupled system Z M in Fig. 10.It can be seen that the spin overlap is vanishing in the high-temperature paramagnetic phase, while it becomes non-zero in the low-temperature phase, signifying RSB.

Conclusions
We introduced a way to compute the SRE [22] with integer Renyi index n > 1 in terms of wavefunction coefficients in many body systems, that make it amenable to efficient computation using Monte Carlo sampling.We applied this approach to generalized Rokhsar-Kivelson systems whose   Hamiltonians allow a stochastic matrix form decomposition [46].Thanks to the known correspondence between ground states of these systems and associated classical statistical mechanics problems, we have been able to express the SRE in terms of classical free energy differences, which can be efficiently computed by thermodynamic integration.Crucially, temperature plays the role of a tunable parameter in the quantum Hamiltonians, allowing us to drive these systems across quantum phase transitions and study the behaviour of their SRE.With this approach we were able to study the SRE of large highdimensional systems, unattainable using existing tensor network-based techniques, and in some cases obtaining explicit analytical results.
We applied this insight to a range of quantum many body SMF Hamiltonians, encompassing the Ising ferromagnet in 1D, 2D, 3D, and infinite dimensions; the J 1 − J 2 model on the square lattice (exhibiting a first order transition); the triangular Ising antiferromagnet (fully frustrated, devoid of ordering); and the Edwards-Anderson model on the cubic lattice (which undergoes a glass transition).Generally, we observed that the behaviour of the SRE is relatively featureless across quantum phase transitions in these systems, although it is indeed singular in its first or higher order derivative, depending on the first or higher order nature of the transition.We found that the maximum of the SRE generically occurs at a cusp away from the quantum critical point, where the derivative suddenly changes sign.Curiously, the cusp appears to occur in the disordered phase in two dimensions, and in the ordered phase in higher dimensions, suggesting that it may be altogether unrelated to the ordering behaviour of the quantum system.
We further compared the SRE to the logarithm of overlaps with specific stabilizer states, that are asymptotically realised in the ground state phase diagrams of these systems.We find that they display strikingly similar behaviors, which in turn establish rigorous bounds on the min-relative entropy of magic.
In our work we were able to make some progress in understanding the behaviour of the magic and its maximum in many body quantum (SMF) systems, throughout their phase diagrams, in terms of partition functions and thermodynamic properties of associated classical problems, and by comparing it with overlaps of asymptotic stabilizer states.One wonders whether further progress could be made using field theoretic approaches for the associated classical problems, in particular ε-expansions just above 2D or just below 3D to shed light on the location of the SRE maximum with respect to the quantum phase transitions.We shall leave these and other interesting open questions for future work.
As we discussed at the end of Sec.5.2, our results for the 2D SMF Ising ferromagnet straightforwardly extend to the toric code [74] and SMF variations thereof [62], which are some of the simplest examples of Z 2 lattice gauge theories.It will be interesting to consider other quantum SMF Hamiltonians constructed from classical systems that exhibit an emergent gauge symmetry, such as dimer [75,76], (spin) ice [77,78,79] and vertex models in general.Simulating these systems in their low temperature phases typically requires the use of loop updates, which pose a nontrivial challenge for the partition function Z M introduced in Sec.

Figure 1 :
Figure 1: Behaviour of various measures of magic, including two stabilizer bounds, for the SMF 1D Ising ferromagnet.For each n, we plot n−1 2n M n /N , which are upper bounded by D x and D zz by Eq. (6).

Figure 2 :
Figure 2: Magic and stabilizer bounds for the SMF 2D Ising ferromagnet.The top panel shows the behaviour of ⟨E M ⟩ M introduced in Sec.4.1 as a stepping stone to compute M 2 /N (bottom panel).The vertical dashed lines indicate the location of the quantum phase transition, whereas the vertical dotted lines indicate the location of the transition in the coupled layered system Z M (resulting in a cusp in the magic M 2 ).

Figure 3 :
Figure 3: Second derivative of M 2 for the 2D SMF Ising ferromagnet.The vertical dashed line indicates the location of the quantum phase transition, whereas the vertical dotted line indicates the location of the transition T * in the coupled layered system Z M .The inset shows the full extent of d 2 (M 2 /N )/dT 2 near T * , which is truncated in the main plot for visualisation purposes.

Figure 4 :
Figure 4: Magic and stabilizer bounds for the SMF 3D Ising ferromagnet.The top panel shows the behaviour of ⟨E M ⟩ M , and the bottom panel shows the SRE density M 2 /N (same system sizes).The vertical dashed lines indicate the location of the quantum phase transition, whereas the vertical dotted lines indicate the location of the transition in the coupled layered system Z M (resulting in a cusp in the magic M 2 ).

=Figure 5 :
Figure 5: Magic and stabilizer bounds for the SMF infinite-range Ising ferromagnet.For each n, we plot n−1 2n M n /N , which are upper bounded by D x and D zz by Eq. (6).The system undergoes a quantum phase transition at T c = 1.

Figure 6 :
Figure 6: Magic and stabilizer bounds for the SMF 2D J 1 -J 2 Ising model.The top panels show the behaviour of ⟨E⟩ and ⟨E M ⟩ M , and the bottom panel shows the SRE density M 2 /N (same system sizes).The vertical dashed lines indicate the location of the quantum phase transition, whereas the vertical dotted lines indicate the location of the transition in the coupled layered system Z M (resulting in a cusp in the magic M 2 ).

Figure 7 :
Figure 7: Magic and stabilizer bounds for the SMF triangular antiferromagnetic Ising model.The top panel shows the behaviour of ⟨E M ⟩ M , and the bottom panel shows the SRE density M 2 /N (same system sizes).The vertical dotted lines indicate the location of the transition in the coupled layered system Z M (resulting in a cusp in the magic M 2 ).

Figure 8 :
Figure 8: Magic and stabilizer bounds for the 3D SMF EA model.The top panels show the behaviour of ⟨E M ⟩ M (left) and of the specific heat C v,M (right), and the bottom panel shows the SRE density M 2 /N (all for the same set of system sizes).The vertical dashed lines indicate the location of the quantum phase transition, whereas the vertical dotted lines indicate the location of the transition at T * in the coupled layered system Z M (resulting in a maximum in the magic M 2 ).The results are averaged over 100, 50, and 30 realizations for linear sizes L = 6, 8, 10, respectively.

Figure 9 :
Figure 9: Energy histograms of the coupled layered system Z M corresponding to the 3D EA model, in thermodynamic equilibrium at different temperatures for a system of size L = 12.The vertical dashed lines are guides to the eye tracking the (same) location location of the two peaks across the panels.

8 Figure 10 :
Figure 10: Spin overlap ⟨q⟩ M for the coupled layered system Z M corresponding to the 3D EA model, for system sizes L ∈ {4, 6, 8}.

i
4.1, and is beyond the scope of the present work.′c σ c * σ ′ ⟨σ ′ |P a,a ′ |σ⟩ a j a ′ j ⟨σ ′ j |X a j Z a ′ j |σ j ⟩