Stabilizer entropies and nonstabilizerness monotones

We study different aspects of the stabilizer entropies (SEs) and compare them against known nonstabilizerness monotones such as the min-relative entropy and the robustness of magic. First, by means of explicit examples, we show that, for R\'enyi index $0\leq n<2$, the SEs are not monotones with respect to stabilizer protocols which include computational-basis measurements, not even when restricting to pure states (while the question remains open for $n\geq 2$). Next, we show that, for any R\'enyi index, the SEs do not satisfy a strong monotonicity condition with respect to computational-basis measurements. We further study SEs in different classes of many-body states. We compare the SEs with other measures, either proving or providing numerical evidence for inequalities between them. Finally, we discuss exact or efficient tensor-network numerical methods to compute SEs of matrix-product states (MPSs) for large numbers of qubits. In addition to previously developed exact methods to compute the R\'enyi SEs, we also put forward a scheme based on perfect MPS sampling, allowing us to compute efficiently the von Neumann SE for large bond dimensions.


Introduction
Stabilizer states and Clifford operations are a pillar of quantum information theory [1,[19][20][21]36].For instance, the stabilizer formalism is at the basis of many prototypical quantum error-correcting codes [16,30], where Clifford operations can be implemented fault tolerantly [42,47].A prominent feature of Clifford circuits and stabilizer states is that they can be efficiently simulated on a classical computer [1,[19][20][21].While this makes them easy to study, it implies that they have to be supplemented with suitable nonstabilizer ancillary states in order to achieve universal quantum computation [4,11].
It is an important task to quantify the degree to which a quantum state cannot be prepared by Clifford gates, i.e., how much it deviates from a stabilizer state.This property, known as nonstabilizer-Tobias Haug: tobias.haug@u.nus.eduness or magic, has received significant attention in the quantum information literature.Similar to quantum entanglement, nonstabilizerness allows for a resource theory [13], where Clifford operators and stabilizer states are free resources.In this context, different nonstabilizerness monotones have been proposed as measures to quantify it [3, 6-8, 28, 29, 35, 51, 52].
More recently, nonstabilizerness was considered in the context of many-body quantum physics [35,38,44,46,53].Here, a major difficulty is that measures of nonstabilizerness are typically very hard to compute (although an efficient measurement protocol for quantum computers has been recently reported [24]).This is especially true in the prototypical case of qubit systems and, more generally, when the local Hilbert space dimension is even [2,9,10].
In order to make progress, useful functions to quantify nonstabilizerness, the so-called stabilizer entropies (SEs), were introduced in Ref. [32], cf. also [12,33,34,37,38].They are expressed in terms of the expectation values of all Pauli strings and, differently from several known monotones, they do not require minimization procedures, making their evaluation simpler.In addition, while the computation of the SE of a typical state is exponentially costly in the system size [38], it can be computed efficiently for the important class of matrix product states (MPSs) [25].This allowed for a detailed study of SEs in one-dimensional systems [25], substantiating previously suggested connections between nonstabilizerness and criticality [53].It is also worth mentioning that SEs can be probed experimentally by randomized measurement protocols [39] or Bell measurements [24,26].
The SE of a pure state |ψ⟩ satisfies the following properties [32]: (i) it is zero if and only if |ψ⟩ is a stabilizer state; (ii) it is invariant under Clifford unitaries; (iii) it is additive under tensor product.SEs can be also defined for mixed states.However, if we define the mixed stabilizer states as the convex hull of the pure stabilizer states, then property (i) does not hold anymore.While the latter can be recovered if a restricted definition of mixed stabilizer states is used [32], in this work we will restrict to the case of pure states.
From the resource-theory point of view, an important open question is whether SEs are monotones with respect to the set of stabilizer protocols [27,51].
Together with Clifford unitaries, these include Pauli measurements and the possibility to discard qubits or operations conditioned on measurement outcomes.More generally, it is also natural to ask how the SE compares to known magic monotones, for instance in terms of order relations between them.
In this work we address these questions, making further progress on the characterization of SEs.We first present a series of general results about monotonicity properties of SEs and their relation with known monotones, namely the min-relative entropies [5,35] and the robustness of magic [28,43,45].For Rényi index 0 ≤ n < 2, we show that the SE is not a monotone for general stabilizer protocols, including protocols involving only pure states.In addition, we show that the SEs do not satisfy a strong monotonicity condition with respect to computational-basis measurements for any Rényi index.This result contradicts a previous statement in the literature stating that the Rényi-1/2 SE is a strong monotone [22], cf. also [23].
Next, we study SEs in different classes of manybody states and compare them against known magic monotones.We compare the SEs with other measures, either proving or providing numerical evidence for inequalities between them.Finally, we discuss tensor-network methods for the efficient computations in MPSs, extending the results of Ref. [25].In particular, we introduce a new computational scheme based on perfect MPS sampling [18].This allows us to access the von Neumann stabilizer entropy, which was outside the reach of the exact method developed in Ref. [25].
Finally, it is important to mention that, independent of their monotonicity properties, SREs have already found numerous applications in different areas of quantum information theory.For instance, they allow us to bound the cost of certain quantumcertification [34] and purity-estimation [33] protocols, they are useful to characterize quantum chaos in many-body states [32,40], and are naturally related to important features of the so-called entanglement spectrum [49,50].Therefore, we expect that the numerical techniques introduced in our work will be practically useful to compute the SRE in different contexts.
The rest of this manuscript is organized as follows.In Sec. 2 we introduce the stabilizer protocols, the SEs and the other monotones of interest in this work.In Sec. 3 we present our counterexamples to (strong) monotonicity for the SEs.In Sec. 4 we study order relations between the SEs and the other measures, both using analytical arguments and providing numerical evidence for small numbers of qubits.Finally, in Sec. 5 we discuss efficient evaluation methods for SEs of MPSs, putting forward a novel approach to compute the von Neumann SE.Our conclusions are consigned to Sec. 6.

Stabilizer Protocols and monotonicity 2.1 Stabilizer protocols
We consider a system of N qubits, with Hilbert space H = ⊗ N j=1 H j , and H j ≃ C 2 .We denote by {σ α } 3 α=0 the Pauli matrices (σ 0 = 1) and by {|0⟩, |1⟩} the local computational basis.Denoting by PN = {i α0 N k=1 σ α k } α0,α1,...,α N the Pauli group, i.e. the group of all N -qubit Pauli strings with phases ±1, ±i, we define the Clifford group to be the set of unitary operators U such that U W U † ∈ PN for all W ∈ PN .Then, pure stabilizer states are the states generated by applying elements of the Clifford group to the reference state |0⟩ ⊗N .Finally, we also intro- ..,α N as the set of Pauli strings with trivial phase +1 only.
According to standard resource theory [13], in order to define nonstabilizerness (or magic) monotones, one needs to first specify the set of free operations.This choice is not unique [27,35], but a minimal set, typically included in the set of free operations, is that of stabilizer protocols.These are quantum channels consisting of the following elementary operations: Given a stabilizer protocol E, a nonstabilizerness monotone M then satisfies In addition to the above, one can define a further property which is sometimes referred to as strong monotonicity [13].
Consider performing computational-basis measurements on a subset Λ containing m qubits.We denote by λ = {λ 1 , . . ., λ m } the set of outcomes for a single measurement (λ j = 0, 1), by ρ λ the post-measurement state, and by p λ = Tr [Π λ ρΠ λ ] the corresponding probability, where (2) The strong-monotonicity condition, with respect to computational-basis measurements, read This relation states that on average the nonstabilizerness cannot increase due to computational-basis measurements.Eq. (3) is not required for a function to be a magic monotone.For instance, the min-relative entropy (introduced in the next section) does not satisfy it.Still, it seems particularly desirable in the many-body setting, as we will discuss later on.As mentioned, this work is restricted to pure states.Therefore, given a quantum state |ψ⟩, in order to test monotonicity of a measure M , we need to consider stabilizer protocols such that for some state |ϕ⟩.This corresponds to a deterministic state-transformation protocol.Crucially, E is not necessarily a unitary channel, and could feature measurements and subsequent Clifford operations conditioned on the outcomes.These feedback operations are necessary in order to make the transformation deterministic, given the randomness of the measurement outcomes.We will discuss an explicit example in Sec.3.1.

SEs and monotones
We now introduce the SEs and two known nonstabilizerness monotones, the min-relative entropy and the robustness of magic.Given a pure state |ψ⟩ ∈ H, the Rényi SE of order n reads [32] M where we use the natural logarithm.We note that it can be rewritten as [32] M where Ξ P (|ψ⟩) = ⟨ψ|P |ψ⟩ 2 /2 N .Since Ξ P (|ψ⟩) ≥ 0 and P ∈P N Ξ P (|ψ⟩) = 1, it can be interpreted as the Rényi-n entropy of the classical probability distribution Ξ P (|ψ⟩), up to an off-set −N log 2. The von Neumann SE is obtained by taking the limit n → 1.
Next, we introduce the so-called log-free robustness of magic [28,35], which is defined as [28] LR where S = {σ i } is the set of pure N -qubit stabilizer states.Note that ρ can be either a pure or a mixed state.The robustness of magic is then defined as Finally, we consider the min-relative entropy of magic [5,35].For a given pure state |ψ⟩, it reads where we introduced the stabilizer fidelity where the maximum is taken over the set of stabilizer states |ϕ⟩.Roughly speaking, D min measures the distance between |ψ⟩ and its nearest stabilizer state.
The robustness and the min-relative entropy of magic are known to be genuine nonstabilizerness monotones [28,35].In addition, the robustness of magic satisfies the strong monotonicity condition [43], while the min-relative entropy does not.Both monotones are natural objects from the resource-theory point of view, and it was argued in Ref. [35] that any "good" magic measure M should satisfy So far, the question of whether the SEs are monotones with respect to stabilizer protocols remained open.In the next section, we will provide an explicit counterexample showing that, for Rényi index 0 ≤ n < 2, they are not.In addition, we will show that, for any Rényi index, the SEs are not strong monotones.Despite these results, the SEs could still be of interest from the resource-theory point of view.For instance, we will discuss in Sec. 4 order relations between the SEs and known monotones.Given the possibility to efficiently compute the SEs, cf.Sec. 5, this could allow one to obtain useful bounds on the nonstabilizerness resource of a given state.

Counterexamples to monotonicity 3.1 Violation of monotonicity for Rényi index 0 ≤ n < 2
We start by presenting our counterexample to the monotonicity condition (1) of the SEs for Rényi index 0 ≤ n < 2. In Appendix A, we will discuss in detail how this counterexample was found.
Consider a set of N = 4 qubits, and take the ordered computational basis We will consider the qubits to be ordered from left to right, so that qubit 1 is the leftmost, etc.We define the state |φ * ⟩ by its coordinates in the basis B, which read so that we have, for instance, etc. Next, consider the protocol defined by the following set of instructions: 1. Perform a projective measurement (in the computational basis) of the first (leftmost) qubit; 2. If the outcome is 1, do nothing; 3. If the outcome is 0, apply the following unitary operator where and Here Note that this is a stabilizer protocol, because all the unitary operators applied after the measurement are Cliffords.In particular, the operator (σ x −σ y )/ √ 2 is a Clifford operator, as it can be seen from the decomposition into Hadamard (H) and S gates, where S = diag(1, i).
Let us apply the protocol E to the input state |φ * ⟩.A priori, E(|φ * ⟩ ⟨φ * |) is a mixed state, because of the random outcome of the measurement.However, the feedback unitary operator (16) has been chosen to make the output state pure (this is true only when the input state is |φ * ⟩), i.e. the output state is independent of the measurement outcome and the state-transformation protocol is deterministic.In particular, by an explicit calculation, one can show and Note in particular that where Finally, we define the difference between the SEs of the input and output state: In the second line, we have used the factorization (21), the fact that the SE is additive and that M n (|1⟩) = 0.The function ∆M (n) can now be straightforwardly evaluated numerically.We plot it in Fig. 1 for n ∈ [0, 15], from which we see that ∆M (n) < 0 for 0 ≤ n < 2, implying a violation of the monotonicity condition (1).
A few comments are in order.First, note that, while not being a stabilizer state, |φ * ⟩ is an eigenstate of the Clifford operation U in (16), i.e.
Ultimately, this makes it possible to find a Clifford feedback operation which makes the statetransformation deterministic, because where we used U σ z 1 = −σ z 1 U .For n ≥ 2, we were not able to find examples with ∆M (n) < 0 for N = 4 qubits.As we explain in Appendix A, the state |φ * ⟩ was found by maximizing numerically the violation of the strong monotonicity condition (3).For N = 4, we find violations for 0 ≤ n < 2. Increasing N , we find violations even for n ≥ 2, cf.also Sec. 3.2.For instance, for N = 5 we find that the strong monotonicity condition can be violated up to n ≃ 3.2.The latter, however, does not imply violation of (1).In particular, none of the nontrivial states that we have found for N = 5 violating the strong monotonicity condition are eigenstates of Clifford operators, so that we could not devise a deterministic protocol as the one presented in this section.Therefore, the question of whether the SEs are monotones for Rényi index n ≥ 2 (at least when restricted to pure states) remains open.

Violation of strong monotonicity
We now show that, for any value of n, the SEs do not satisfy the strong monotonicity condition (3).First, we show this for 0 ≤ n < 2. To this end, it is enough to consider the counterexample (13).Indeed, let us define where Therefore, because of (23), we have Using the results of Sec.3.1, we have M n (|ψ * ⟩) > M n (|φ * ⟩) for 0 ≤ n < 2. This proves that strong monotonicity is violated for 0 ≤ n < 2.
Next, let us show that the condition (3) is also violated for n ≥ 2. To this end, let us consider a system of N qubits, and define the state where |χ⟩ is the magic state [4] |χ⟩ = e −i(π/4) cos β|0⟩ + sin β|1⟩ , (31) with cos(2β) = 1/ √ 3, and Note that We claim that, for finite ε > 0 and n > 1, we have where c n (ε) is a constant which depends on n and ε, but not on N .We prove this in Appendix A, where we also support this statement numerically.Consider now performing a projective measurement of the first qubit.As usual, we denote by |ψ ε 0 ⟩, |ψ ε 1 ⟩ the post-measurement states associated with the outcomes 0 and 1, and by p 0 , p 1 the corresponding probabilities.When the measurement outcome is 1, the post-measurement state is Because of additivity and the fact that M n (|1⟩) = 0, we have so that Finally, we note that p 1 remains finite in the thermodynamic limit N → ∞, Combining Eqs.(34), (37), and (38), we obtain a violation of the strong monotonicity condition for sufficiently large N and n > 1 (and thus also for n ≥ 2).
In passing, we mention that this example can also be used to show that the min-relative entropy D min is not a strong monotone.Moreover, it illustrates how, in the many-body setting, strong monotonicity appears to be particularly important.Indeed, if strong monotonicity is absent, we see here that measuring a single qubit is enough to increase nonstablizerness from a O(1) constant to a O(N ) value with finite probability.
Finally, let us mention that the results presented in this section contradict previous claims in the literature.This is because the Rényi-1/2 SE coincides with the measure introduced in Ref. [22], which was claimed to satisfy the strong monotonicity condition (3), cf. also [23].

Relations to other monotones
As mentioned, even if they fail to be monotones, SEs could still be of interest from the resource-theory point of view, based on their relation with known monotones.In this section we discuss elementary inequalities between the SEs, the min-relative entropy and the robustness of magic (as usual, restricting to pure states).Note that the latter two are known to satisfy the general inequality [35] First, we recall the following inequality [28,32] For n > 1 we now derive an upper bound also in terms of the min-relative entropy, namely To see this, we note that for any state |ψ⟩ with stabilizer fidelity F STAB (|ψ⟩), we can find a Clifford unitary where This follows immediately from In Eq. (42), we have denoted by {|n⟩} the set of states in the computational basis.Next, we have Here P is the set of Pauli strings while P z is the set of Pauli strings containing 1 and σ z only.In the second line we used the convexity inequality Now, using and Eq. (44) yields (48) Finally, using that the SE is invariant under unitary Clifford operations, we have M n (|ψ⟩) = M n (|ϕ⟩) and, combining with Eq. (48), we arrive at Eq. (41).
Given the inequalities (39) and (41), it would be very useful to also provide an N -independent lower bound of the SE in terms of either the log-robustness or the min-relative entropy.In the following, we provide a simple example showing that this is not possible for the log-robustness of magic if n > 1/2.Consider the single-qubit state For small s, we can compute where a n , b n are n-dependent constants.Now, introducing the N -qubit state and using additivity of the SE, we have Therefore, using (40), we have that, for n > 1/2 (55) This shows that it is not possible to derive an Nindependent lower bound for n > 1/2.
On the other hand, since D min (|Ω(s)⟩) ∝ s 2 , and given the sub-additivity of the min-relative entropy, this example does not rule out the possibility that the SE with n > 1 can be lower-bounded by D min .Since M n vanishes for n → ∞, this bound should necessarily involve a n-dependent factor.In order to avoid dealing with this, we have performed extensive numerical minimization of the fraction M n (|ψ⟩)/D min (|ψ⟩) for small values of n up to N = 4 qubits.We have found numerical evidence that the minimized value of the ratio appears to be a constant, which only depends on n.In particular, we found M n (|ψ⟩) ≳ 1.7D min (|ψ⟩) for all states |ψ⟩ of up to N = 4 qubits and 1 ≤ n ≤ 2. For higher n, the minimized value for N ≤ 4 appears to scale as O(1/n), which matches the overall scaling of M n ∼ O(1/n) for n ≫ 1.Overall, our results hint at the possibility that such a lower bound exists.We leave this as an open question for future research.

Numerical methods for the SE
One of the reasons why SEs are appealing is the fact that they can often be computed much more efficiently than previously known nonstabilizerness monotones.In fact, even if for generic states the computational cost grows exponentially in N [32,38], the evaluation of Pauli expectation values is, in practice, simpler than carrying out the minimization procedure involved in the definition of the robustness and minrelative entropy of magic.In addition, for certain classes of many-body states, the computational cost is polynomial in the number of qubits N .This is the case for MPSs [14,15,41], the simplest example of tensor-network states [48].This was pointed out in Ref. [25], where an efficient computational method was put forward for the computation of Rényi SEs with integer index n > 1.
The aim of this section is to review the result of Ref. [25], and complement it putting forward an approach to the computation of the von Neumann SE.Combined with the bounds discussed in Sec. 4, these methods could be practically useful from the resourcetheory point of view, when the number of qubits makes the computation of the min-relative entropy or robustness of magic unfeasible.
First, we recall the definition of an MPS (with open boundary conditions) where A s k are χ k × χ k+1 matrices, with χ 1 = χ N +1 = 1.We call χ = max j χ j the bond-dimension.The method developed in Ref. [25] starts from the simple identity where 3 α=0 (σ α j ⊗ σ α j ) ⊗n , while (•) denotes complex conjugation.This formula allows to one replace the sum over exponentially many terms with the computation of the norm of an MPS of bonddimension χ 2n , which can be carried out efficiently in N , at a cost O(N χ 6n ) [25].This mapping is also convenient from the analytic point of view.For instance, it allows one to show that the density of SE can be extracted locally for MPSs without long-range correlations.
This method does not allow for a numerical evaluation of the SE with arbitrary Rényi index n.However, for n = 1 an efficient numerical scheme is possible exploiting once again the properties of MPSs.To this end, we start from the expression (58) and interpret Ξ P (|ψ⟩) as a probability distribution.Since each term in the sum is positive, we can evaluate M 1 (|ψ⟩) by sampling the probability distribution Ξ P (|ψ⟩).Crucially, for an MPS |Ψ N ⟩, we do not need a Monte Carlo approach, but we can sample from Ξ P (|Ψ N ⟩) exactly, slightly generalizing the algorithm of perfect MPS sampling of Ref. [18].
For an MPS with bond dimension χ, we can use the standard prescriptions [17,48] to evaluate the optimal contraction cost of the tensor network associated with (60), which yields O(χ 3 N ).Therefore, the conditional probabilities can be computed efficiently in N .Note that the dependence with χ is favorable compared to that needed in the exact computation of the Rényi SEs for integer n [25].
Given (59), we can use directly the algorithm detailed in Ref. [18] to generate strings of Pauli operators associated with α, according to a probability which is exactly Ξ Pα .Then, we can evaluate M 1 (|Ψ N ⟩) by averaging the value log Ξ Pα (|Ψ N ⟩) over all the generated samples.
The accuracy of the method increases with number of samples, which we denote by S (in order not to distinguish it from the number of qubits N ).Denoting by ε S the difference between the average value, obtained by sampling the distribution Ξ P , and the actual value of M 1 , standard arguments yield ε S = O(1/ √ S).In turn, this implies that M 1 can be computed efficiently, i.e. at a computational cost scaling polynomially in N , for fixed accuracy.
We briefly illustrate this method for the groundstate of the interacting XXZ Heisenberg Hamiltonian where ∆ is the anisotropy parameter, while we choose open boundary conditions.Note that the to- In Fig. 2, we report an example of our data for the SE of the ground state of the Heisenberg model at half-filling, i.e. the ground state within the subspace of states with C |ψ⟩ = 0.In Fig. 2a, we plot the SE density m n = M n /N against ∆ for the von Neumann SE density m 1 and the Rényi-2 SE density.Qualitatively, both m 1 and m 2 behave similar as a function of ∆, where m is first increasing with ∆, then decreases.Further, both m 1 and m 2 increase in value with N , where we expect them to converge to a constant for large N as the SE is extensive.Quantitatively, the decrease in m towards ∆ → 1 is more pronounced for m 2 .In Fig. 2b, we show m 1 for larger values of N , illustrating the efficiency of the method.We see that m 1 decreases near ∆ = 1 more sharply with higher N .
In Appendix B we provide further details on the numerical accuracy on the values m 1 and m 2 , as a function of the bond dimension χ and number of samples S.
Finally, we note that one could in principle extend the perfect sampling method to arbitrary Rényi SEs.However, we expect that, in general, the number of samples required to obtained M n (|ψ⟩) up to fixed accuracy could scale exponentially in the number of qubits N .

Outlook
We have studied different aspects of the SEs, comparing them against the min-relative entropy and the robustness of magic.For Rényi index 0 ≤ n < 2, we have shown that the SEs are not monotones with respect to generic stabilizer protocols, not even when restricting to pure states.In addition, we have shown that, for any Rényi index, the SEs do not satisfy the strong monotonicity condition with respect to computational-basis measurements, contradicting a previous claim in the literature [22], cf. also [23].
Next, we have discussed, both analytically and numerically, some inequalities between the SEs, the minrelative entropy and the robustness of magic.Finally, we have presented available tensor-network methods for the efficient computation of SEs in the context of MPSs, reviewing the results of Ref. [25] and putting forward a new method for the von Neumann SE.Our work raises several questions.First, it would be interesting to either generalize our counterexample to show rigorously that the monotonicity condition (1) is violated also for Rényi index n ≥ 2, or, conversely, to prove it.Second, as we discuss in Sec. 4, we believe it would be important to understand whether one can derive a N -independent lower bound for the SE in terms of the min-relative entropy, at least for some Rényi index.Indeed, even if the SE is not a monotone, this result, together with those presented in Sec. 4, could be used to provide upper and lower bounds to a genuine nonstabilizerness monotone, the min-relative entropy.
Finally, our results show that the SEs are not monotones with respect to general stabilizer protocols (at least for Rényi index 0 ≤ n < 2), and in any case not strong monotones.Therefore, our work highlights the need to find genuine nonstabilizerness monotones which are efficient to compute, at least in restricted classes of pure states such as MPSs.In this context, a very natural question is whether such a monotone can be defined in terms of the Pauli spectrum [3], similarly to the SEs.We hope that our work will motivate further research in this direction.
Note added : While finalizing this manuscript, the work [31] appeared on the arXiv, where a similar method to compute SEs in MPSs was put forward.
Next, we note that the state |ψ ε ⟩ is invariant under permutation of qubits.Therefore, we can simplify (65) In addition, the expectation value of Pauli matrices can be evaluated explicitly, yielding The number of terms in the sum of Eq. (65) scales polynomially in N , and can be evaluated in a few seconds for up to N ∼ 50 qubits.We have done this for different values of n > 1, and confirmed that M n (|ψ ε ⟩) is bounded in N .In fact, it is not difficult, although a bit cumbersome, to prove this explicitly, as we now sketch.
We can write G n (|ψ ε ⟩) ≃ (I)/N 2n ε + (II)/N 2n ε , where and First, we note that where δ N = |2ε cos(β) N cos(πN/4)| is an exponentially vanishing term.Next, we show that (II) is exponentially small.Indeed, using the triangular inequality for the absolute value, and the convexity inequality For n > 1 we see that (II) is exponentially small in N .Therefore, using log N ε ≃ log(1 + ε 2 ) (up to exponentially small terms in N ) we have where δN is an exponentially small correction.We have verified (72), based on exact evaluation of (65).
We report an example of our numerical data in Fig. 3, showing excellent agreement.

B Additional numerical results
We provide further details about the numerical results discussed in Sec. 5. First, we present additional numerical data about the perfect-MPS sampling method for the von Neumann SE.In Fig. 4 we show how the approximation error decreases with the number of samples S. We plot the average difference | m1 − m 1 |, where m1 is an estimate for m 1 , computed by averaging over a given number of samples S, while m 1 is our best prediction, computed by averaging over a very large value of samples, S max (in this case, S max = 10 5 ).The plot clearly shows the scaling | m1 − m 1 | ∼ S −1/2 , as expected.Next, we study the dependence of the estimated values of m 1 and m 2 with the bond dimensions for the ground-state of the Heisenberg model at half-filling.Our data are reported in Fig. 5.In Fig. 5a, we show the Rényi-2 SE density M 2 /N = m 2 as function of bond dimension χ, which was computed using the method of Ref. [25].We find that the value of m 2 decreases and converges to a constant for increasing χ.In Fig. 5b, we show a similar plot for the SE, computed with the method presented in Sec. 5. Once again, we see convergence as the bond dimension is increased.Finally, in Fig. 5c, we show the fidelity F = | ⟨Ψ N (χ)|GS⟩ | 2 between the MPS approximation with bond dimension χ, |Ψ N (χ)⟩, and the true ground state |GS⟩.The true ground-state |GS⟩ is estimated by taking a very large bond dimension (in this case, χ ≃ 200).

Figure 4 :
Figure 4: Average error in estimation of the von Neumann SE density | m1 − m1| against number of samples S. m1 is the estimated SE density, while m1 is a highly accurate estimation from S = 10 5 samples.The error is averaged over 100 random instances.We show the ground state of the Heisenberg model for half-filling, ∆ = 0.95 and different qubit number N .