# Stabilizer entropies and nonstabilizerness monotones

Tobias Haug¹ and Lorenzo Piroli²

¹QOLS, Blackett Laboratory, Imperial College London SW7 2AZ, UK

²Philippe Meyer Institute, Physics Department, École Normale Supérieure (ENS), Université PSL, 24 rue Lhomond, F-75231 Paris, France
August 22, 2023

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ényi index $0 \le 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 \ge 2$). Next, we show that, for any Rényi 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ényi 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.

## 1 Introduction

Stabilizer states and Clifford operations are a pillar of quantum information theory [1, 19–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–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.edu

ness 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 $|\psi\rangle$ satisfies the following properties [32]: (i) it is zero if and only if $|\psi\rangle$ 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 \le 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 many-body 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 quantum-certification [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.

# 2 Stabilizer Protocols and monotonicity

## 2.1 Stabilizer protocols

We consider a system of $N$ qubits, with Hilbert space $\mathcal{H} = \bigotimes_{j=1}^{N} \mathcal{H}_j$, and $\mathcal{H}_j \simeq \mathbb{C}^2$. We denote by $\{\sigma^\alpha\}_{\alpha=0}^3$ the Pauli matrices ($\sigma^0 = 1$) and by $\{|0\rangle, |1\rangle\}$ the local computational basis. Denoting by $\tilde{P}_N = \{i^{\alpha_0} \otimes_{k=1}^{N} \sigma^{\alpha_k}\}_{\alpha_0, \dots, \alpha_N}$ the Pauli group, i.e. the group of all $N$-qubit Pauli strings with phases $\pm 1$, $\pm i$, we define the Clifford group to be the set of unitary operators $U$ such that $UWU^\dagger \in \tilde{P}_N$ for all $W \in \tilde{P}_N$. Then, pure stabilizer states are the states generated by applying elements of the Clifford group to the reference state $|0\rangle^{\otimes N}$. Finally, we also introduce $\mathcal{P}_N = \{\otimes_{k=1}^{N} \sigma^{\alpha_k}\}_{\alpha_1, \dots, \alpha_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:

1. Clifford unitary operations, $\rho \mapsto U\rho U^\dagger$, where $U$ is in the Clifford group;

2. Composition with stabilizer states, $\rho \mapsto \rho \otimes S$, where $S$ is a stabilizer state;

3. Measurements in the computational basis;

4. Discarding of qudits;

5. The above operations conditioned on the outcomes of measurements.

Given a stabilizer protocol $\mathcal{E}$, a nonstabilizerness monotone $M$ then satisfies

$$M[\mathcal{E}(\rho)] \le M(\rho). \qquad (1)$$

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 $\Lambda$ containing $m$ qubits. We denote by $\lambda = \{\lambda_1, \dots, \lambda_m\}$ the set of outcomes for a single measurement ($\lambda_j = 0, 1$), by $\rho_\lambda$ the post-measurement state, and by $p_\lambda = \text{Tr}[\Pi_\lambda \rho \Pi_\lambda]$ the corresponding probability, where

$$\Pi_\lambda = |\lambda\rangle\langle\lambda| \otimes \mathbb{1}_{N\setminus m}. \qquad (2)$$

The strong-monotonicity condition, with respect to computational-basis measurements, read

$$M(\rho) \ge \sum_{\lambda} p_{\lambda} M[\rho_{\lambda}]. \qquad (3)$$
This relation states that on average the nonstabilizer-ness 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 $|\psi\rangle$, in order to test monotonicity of a measure $M$, we need to consider stabilizer protocols such that

$$ \mathcal{E}(|\psi\rangle\langle\psi|) = |\phi\rangle\langle\phi|, \qquad (4) $$

for some state $|\phi\rangle$. This corresponds to a deterministic state-transformation protocol. Crucially, $\mathcal{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.

## 2.2 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 $|\psi\rangle \in \mathcal{H}$, the Rényi SE of order $n$ reads [32]

$$ M_n(|\psi\rangle) = (1-n)^{-1} \log \sum_{P \in \mathcal{P}_N} \frac{\langle\psi|P|\psi\rangle^{2n}}{2^N}, \qquad (5) $$

where we use the natural logarithm. We note that it can be rewritten as [32]

$$ M_n(|\psi\rangle) = (1-n)^{-1} \log \sum_{P \in \mathcal{P}_N} \Xi_P^n(|\psi\rangle) - N \log 2, \qquad (6) $$

where $\Xi_P(|\psi\rangle) = \langle\psi|P|\psi\rangle^2/2^N$. Since $\Xi_P(|\psi\rangle) \ge 0$ and $\sum_{P \in \mathcal{P}_N} \Xi_P(|\psi\rangle) = 1$, it can be interpreted as the Rényi-$n$ entropy of the classical probability distribution $\Xi_P(|\psi\rangle)$, up to an off-set $-N \log 2$. The von Neumann SE is obtained by taking the limit $n \to 1$.

Next, we introduce the so-called log-free robustness of magic [28, 35], which is defined as [28]

$$ \mathrm{LR}(\rho) = \log \left[ \min_x \left\{ \sum_i |x_i| : \rho = \sum_i x_i \sigma_i \right\} \right], \qquad (7) $$

where $S = \{\sigma_i\}$ is the set of pure $N$-qubit stabilizer states. Note that $\rho$ can be either a pure or a mixed state. The robustness of magic is then defined as

$$ R(\rho) = e^{\mathrm{LR}(\rho)}. \qquad (8) $$

Finally, we consider the min-relative entropy of magic [5, 35]. For a given pure state $|\psi\rangle$, it reads

$$ D_{\min}(|\psi\rangle) = -\log[F_{\mathrm{STAB}}(|\psi\rangle)], \qquad (9) $$

where we introduced the stabilizer fidelity

$$ F_{\mathrm{STAB}}(|\psi\rangle) = \max_{|\phi\rangle} \{|\langle\psi|\phi\rangle|^2\}, \qquad (10) $$

where the maximum is taken over the set of stabilizer states $|\phi\rangle$. Roughly speaking, $D_{\min}$ measures the distance between $|\psi\rangle$ 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

$$ D_{\min}(\rho) \le M(\rho) \le \mathrm{LR}(\rho). \qquad (11) $$

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 \le 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.

# 3 Counterexamples to monotonicity

## 3.1 Violation of monotonicity for Rényi index $0 \le n < 2$

We start by presenting our counterexample to the monotonicity condition (1) of the SEs for Rényi index $0 \le 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

$$ \mathcal{B} = \{|0000\rangle, |0001\rangle, |0010\rangle, |0011\rangle, |0100\rangle, \\ \dots, |1110\rangle, |1111\rangle\}. \qquad (12) $$

We will consider the qubits to be ordered from left to right, so that qubit 1 is the leftmost, etc. We define the state $|\varphi^*\rangle$ by its coordinates in the basis $\mathcal{B}$, which read

$$ \frac{1}{2\sqrt{6}} \{0, 0, 0, 2i, 1-i, -1-i, -1+i, -1+i, \\ \sqrt{2}, \sqrt{2}i, -\sqrt{2}, -\sqrt{2}, 0, 0, 0, \sqrt{2}(1+i)\}, \qquad (13) $$

so that we have, for instance,

$$ \langle 0000 | \varphi^* \rangle = \langle 1110 | \varphi^* \rangle = 0, \qquad (14) $$

$$ \langle 0011 | \varphi^* \rangle = i/\sqrt{6}, \qquad (15) $$
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

$$U = V_2 V_1, \qquad (16)$$

where

$$V_1 = \sigma_1^x \otimes \frac{(\sigma_2^x - \sigma_2^y)}{\sqrt{2}} \otimes \mathbb{1}_3 \otimes \mathbb{1}_4, \qquad (17)$$

$$V_2 = \mathbb{1}_1 \otimes \mathbb{1}_2 \otimes w_{34}, \qquad (18)$$

and

$$w_{34} = (\sigma_3^x \otimes \mathbb{1}_4) \text{CZ}_{34}(\sigma_3^x \otimes \mathbb{1}_4). \qquad (19)$$

Here CZ$_{34}$ is the (Clifford) controlled-Z gate. Namely, in the basis $\{|00\rangle_{34}, |01\rangle_{34}, |10\rangle_{34}, |11\rangle_{34}\}$, we have CZ$_{34}$ = diag(1, 1, 1, -1).

Note that this is a stabilizer protocol, because all the unitary operators applied after the measurement are Cliffords. In particular, the operator $(\sigma^x - \sigma^y)/\sqrt{2}$ is a Clifford operator, as it can be seen from the decomposition into Hadamard ($H$) and $S$ gates,

$$(\sigma^x - \sigma^y)/\sqrt{2} = \text{HSHSH}\sigma^x, \qquad (20)$$

where $S = \text{diag}(1, i)$.

Let us apply the protocol $\mathcal{E}$ to the input state $|\varphi^*\rangle$. A priori, $\mathcal{E}(|\varphi^*\rangle\langle\varphi^*|)$ 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 $|\varphi^*\rangle$), 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 $\mathcal{E}(|\varphi^*\rangle\langle\varphi^*|) = |\psi^*\rangle\langle\psi^*|$ where

$$|\psi^*\rangle = |1\rangle_1 \otimes |\chi^*\rangle_{234}, \qquad (21)$$

and

$$|\chi^*\rangle = \frac{1}{\sqrt{6}} (|000\rangle + i|001\rangle - |010\rangle - |011\rangle + (1+i)|111\rangle). \qquad (22)$$

Note in particular that

$$|\psi^*\rangle = \frac{1}{\sqrt{p_1}} \Pi_1 |\varphi^*\rangle, \qquad (23)$$

where $p_1 = \langle\varphi^*|\Pi_1|\varphi^*\rangle$, and $\Pi_1 = |1\rangle\langle1| \otimes \mathbb{1}_{234}$.

Finally, we define the difference between the SEs of the input and output state:

$$\begin{aligned} \Delta M(n) := & M_n(|\varphi^*\rangle) - M_n(|\psi^*\rangle) \\ = & M_n(|\varphi^*\rangle) - M_n(|\chi^*\rangle). \end{aligned} \qquad (24)$$

In the second line, we have used the factorization (21), the fact that the SE is additive and that $M_n(|1\rangle) = 0$. The function $\Delta M(n)$ can now be straightforwardly evaluated numerically. We plot it in Fig. 1 for $n \in [0, 15]$, from which we see that $\Delta M(n) < 0$ for $0 \le 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, $|\varphi^*\rangle$ is an eigenstate of the Clifford operation $U$ in (16), i.e.

$$|\varphi^*\rangle = U |\varphi^*\rangle. \qquad (25)$$

Ultimately, this makes it possible to find a Clifford feedback operation which makes the state-transformation deterministic, because

$$\begin{aligned} U |0\rangle\langle0|_{1} |\varphi^*\rangle &= U \frac{(1 + \sigma_1^z)}{2} |\varphi^*\rangle \\ &= \frac{(1 - \sigma_1^z)}{2} U |\varphi^*\rangle = |1\rangle\langle1| |\varphi^*\rangle, \end{aligned} \qquad (26)$$

where we used $U\sigma_1^z = -\sigma_1^z U$.

For $n \ge 2$, we were not able to find examples with $\Delta M(n) < 0$ for $N = 4$ qubits. As we explain in Appendix A, the state $|\varphi^*\rangle$ was found by maximizing numerically the violation of the strong monotonicity condition (3). For $N = 4$, we find violations for $0 \le n < 2$. Increasing $N$, we find violations even for $n \ge 2$, cf. also Sec. 3.2. For instance, for $N = 5$ we find that the strong monotonicity condition can be violated up to $n \simeq 3.2$. The latter, however, does not imply violation of (1). In particular, none of the non-trivial 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 monones for Rényi index $n \ge 2$ (at least when restricted to pure states) remains open.

## 3.2 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 \le n < 2$. To this end, it is enough to consider the counterexample (13). Indeed, let us define

$$|\varphi_0^*\rangle = \frac{1}{\sqrt{p_0}} \Pi_0 |\varphi^*\rangle, \quad |\varphi_1^*\rangle = \frac{1}{\sqrt{p_1}} \Pi_1 |\varphi^*\rangle, \qquad (27)$$

where $\Pi_a = |a\rangle\langle a| \otimes \mathbb{1}_{234}$ and $p_a = \langle\varphi^*|\Pi_a|\varphi^*\rangle$. It is straightforward to show

$$M_n(|\varphi_0^*\rangle) = M_n(|\varphi_1^*\rangle). \qquad (28)$$

Therefore, because of (23), we have

$$p_0 M_n(|\varphi_0^*\rangle) + p_1 M_n(|\varphi_1^*\rangle) = M_n(|\psi^*\rangle). \qquad (29)$$
<figure><img src="image_1.png" /><figcaption>Figure 1: Numerical evaluation of the difference ΔM(n) defined in Eq. (24), showing ΔM(n) < 0 for 0 ≤ n < 2. Inset: same plot for n ∈ [0, 15].</figcaption></figure>

Using the results of Sec. 3.1, we have $M_n(|\psi^*\rangle) > M_n(|\varphi^*\rangle)$ for $0 \le n < 2$. This proves that strong monotonicity is violated for $0 \le n < 2$.

Next, let us show that the condition (3) is also violated for $n \ge 2$. To this end, let us consider a system of N qubits, and define the state

$$|\psi^{\epsilon}\rangle = \frac{1}{\sqrt{N_{\epsilon}}} \left[ |0\rangle^{\otimes N} + \epsilon |\chi\rangle^{\otimes N} \right], \quad (30)$$

where $|\chi\rangle$ is the magic state [4]

$$|\chi\rangle = e^{-i(\pi/4)} \cos \beta |0\rangle + \sin \beta |1\rangle, \quad (31)$$

with $\cos(2\beta) = 1/\sqrt{3}$, and

$$N_{\epsilon} = 1 + \epsilon^2 + 2\epsilon \cos(\beta)^N [\cos(N\pi/4)]. \quad (32)$$

Note that

$$|\chi\rangle \langle\chi| = \frac{1}{2} \left[ 1_1 + \frac{1}{\sqrt{3}} (\sigma^x + \sigma^y + \sigma^z) \right]. \quad (33)$$

We claim that, for finite $\epsilon > 0$ and $n > 1$, we have

$$M_n(|\psi^{\epsilon}\rangle) < c_n(\epsilon), \quad (n > 1) \quad (34)$$

where $c_n(\epsilon)$ is a constant which depends on n and $\epsilon$, 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 $|\psi_0^\epsilon\rangle$, $|\psi_1^\epsilon\rangle$ 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

$$|\psi_1^\epsilon\rangle = |1\rangle_1 \otimes |\chi\rangle^{\otimes(N-1)}. \quad (35)$$

Because of additivity and the fact that $M_n(|1\rangle) = 0$, we have

$$M_n(|\psi_1^\epsilon\rangle) = (N-1)M_n(|\chi\rangle), \quad (36)$$

so that

$$\begin{aligned} p_0 M_n(|\psi_0^\epsilon\rangle) + p_1 M_n(|\psi_1^\epsilon\rangle) &\ge p_1 M_n(|\psi_1^\epsilon\rangle) \\ &= p_1(N-1)M_n(|\chi\rangle). \end{aligned} \quad (37)$$

Finally, we note that $p_1$ remains finite in the thermodynamic limit $N \to \infty$,

$$\lim_{N \to \infty} p_1 = \lim_{N \to \infty} \frac{1}{N_{\epsilon}^{\epsilon^2} \sin^2 \beta} = \frac{\epsilon^2}{1 + \epsilon^2} \sin^2 \beta. \quad (38)$$

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 \ge 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].

## 4 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]

$$D_{\min}(|\psi\rangle) \le \mathrm{LR}(|\psi\rangle). \quad (39)$$

First, we recall the following inequality [28, 32]

$$M_n(|\psi\rangle) \le 2\mathrm{LR}(|\psi\rangle) \quad (n \ge 1/2). \quad (40)$$

For $n > 1$ we now derive an upper bound also in terms of the min-relative entropy, namely

$$M_n(|\psi\rangle) \le \frac{2n}{n-1} D_{\min}(|\psi\rangle). \quad (41)$$

To see this, we note that for any state $|\psi\rangle$ with stabilizer fidelity $F_{\mathrm{STAB}}(|\psi\rangle)$, we can find a Clifford unitary $U_C$ s.t.

$$|\phi\rangle = U_C |\psi\rangle = \sum_k a_k |k\rangle \quad (42)$$

where

$$|a_0|^2 = F_{\mathrm{STAB}}(|\psi\rangle). \quad (43)$$
This follows immediately from $F_{\text{STAB}} = |\langle\psi|\phi_{\text{max}}\rangle|^2 = |\langle\psi|U_{\text{C}}^{\text{max}}|0\rangle|^2 = |a_0|^2$. In Eq. (42), we have denoted by $\{|n\rangle\}$ the set of states in the computational basis. Next, we have

$$2^{-N} \sum_{P \in \mathcal{P}_N} |\langle \phi | P | \phi \rangle|^{2n} \ge 2^{-N} \sum_{P \in \mathcal{P}_z} |\langle \phi | P | \phi \rangle|^{2n} \\ \ge 2^{-2nN} \left( \sum_{P \in \mathcal{P}_z} |\langle \phi | P | \phi \rangle| \right)^{2n} . \quad (44)$$

Here $\mathcal{P}$ is the set of Pauli strings while $\mathcal{P}_z$ is the set of Pauli strings containing 1 and $\sigma^z$ only. In the second line we used the convexity inequality

$$\sum_{i=1}^{m} |a_i|^k \ge \frac{1}{m^{k-1}} \left( \sum_{i=1}^{m} |a_i| \right)^k . \quad (45)$$

Now, using

$$\left( \sum_{P \in \mathcal{P}_z} |\langle \phi | P | \phi \rangle| \right)^{2n} \ge \left| \sum_{P \in \mathcal{P}_z} \langle \phi | P | \phi \rangle \right|^{2n} \quad (46)$$

and

$$\sum_{P \in \mathcal{P}_z} P = 2^N |0\rangle \langle 0|, \quad (47)$$

Eq. (44) yields

$$2^{-N} \sum_{P \in \mathcal{P}_N} |\langle \phi | P | \phi \rangle|^{2n} \ge |a_0|^{4n} = F_{\text{STAB}}(|\psi\rangle)^{2n} . \quad (48)$$

Finally, using that the SE is invariant under unitary Clifford operations, we have $M_n(|\psi\rangle) = M_n(|\phi\rangle)$ 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

$$|\Omega(s)\rangle = \cos(s) |0\rangle + \sin(s) |1\rangle . \quad (49)$$

For small $s$, we can compute

$$M_n(|\Omega(s)\rangle) = a_n s^2 + o(s^2), \quad n > 1, \quad (50)$$

$$M_n(|\Omega(s)\rangle) = b_n s^{2n} + o(s^{2n}), \quad n < 1, \quad (51)$$

where $a_n, b_n$ are $n$-dependent constants. Now, introducing the $N$-qubit state

$$|\Lambda_N(s_0)\rangle = |\Omega(s_0/\sqrt{N})\rangle^{\otimes N}, \quad (52)$$

and using additivity of the SE, we have

$$M_n(|\Lambda_N(s_0)\rangle) = a_n s_0^2 + o(1), \quad n > 1, \quad (53)$$

$$M_n(|\Lambda_N(s_0)\rangle) = b_n s_0^{2n} N^{1-n} + o(N^{1-n}), \quad n < 1. \quad (54)$$

Therefore, using (40), we have that, for $n > 1/2$

$$\frac{\text{LR}(|\Lambda_N(s_0)\rangle)}{M_n(|\Lambda_N(s_0)\rangle)} \ge \frac{1}{2} \frac{M_{1/2}(|\Lambda_N(s_0)\rangle)}{M_n(|\Lambda_N(s_0)\rangle)} \sim N^{\min(n-\frac{1}{2}, \frac{1}{2})}. \quad (55)$$

This shows that it is not possible to derive an $N$-independent lower bound for $n > 1/2$.

On the other hand, since $D_{\text{min}}(|\Omega(s)\rangle) \propto 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_{\text{min}}$. Since $M_n$ vanishes for $n \to \infty$, 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(|\psi\rangle)/D_{\text{min}}(|\psi\rangle)$ 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(|\psi\rangle) \gtrsim 1.7 D_{\text{min}}(|\psi\rangle)$ for all states $|\psi\rangle$ of up to $N = 4$ qubits and $1 \le n \le 2$. For higher $n$, the minimized value for $N \le 4$ appears to scale as $O(1/n)$, which matches the overall scaling of $M_n \sim O(1/n)$ for $n \gg 1$. Overall, our results hint at the possibility that such a lower bound exists. We leave this as an open question for future research.

## 5 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 monotonites. 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 min-relative 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 resource-theory 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) $|\Psi_N\rangle$ [41]

$$|\Psi_N\rangle = \sum_{\{s_k\}} A_1^{s_1} \dots A_N^{s_N} |s_1, \dots, s_N\rangle, \quad (56)$$
where $A_k^s$ are $\chi_k \times \chi_{k+1}$ matrices, with $\chi_1 = \chi_{N+1} = 1$. We call $\chi = \max_j \chi_j$ the bond-dimension. The method developed in Ref. [25] starts from the simple identity

$$ \sum_{P \in \mathcal{P}_N} \frac{\langle \Psi_N | P | \Psi_N \rangle^{2n}}{2^N} = (\langle \Psi_N | \otimes \langle \overline{\Psi}_N | \rangle^{\otimes n} \Lambda_1^{(n)} \otimes \Lambda_2^{(n)} \otimes \cdots \otimes \Lambda_N^{(n)} (|\Psi_N\rangle \otimes |\overline{\Psi}_N\rangle))^{\otimes n}, \quad (57) $$

where $\Lambda_j^{(n)} = (1/2) \sum_{\alpha=0}^3 (\sigma_j^\alpha \otimes \overline{\sigma}_j^\alpha)^{\otimes n}$, while $(\cdot)$ 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 bond-dimension $\chi^{2n}$, which can be carried out efficiently in $N$, at a cost $O(N\chi^{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

$$ M_1(|\psi\rangle) = - \sum_{P \in \mathcal{P}_N} \Xi_P(|\psi\rangle) \log \Xi_P(|\psi\rangle) - N \log 2, \quad (58) $$

and interpret $\Xi_P(|\psi\rangle)$ as a probability distribution. Since each term in the sum is positive, we can evaluate $M_1(|\psi\rangle)$ by sampling the probability distribution $\Xi_P(|\psi\rangle)$. Crucially, for an MPS $|\Psi_N\rangle$, we do not need a Monte Carlo approach, but we can sample from $\Xi_P(|\Psi_N\rangle)$ exactly, slightly generalizing the algorithm of perfect MPS sampling of Ref. [18].

Let us define $\alpha = (\alpha_1, \dots, \alpha_N)$, with $\alpha_j = 0, 1, 2, 3$. The key observation is that the probability of the Pauli string associated with $\alpha$, $P_\alpha = \sigma_1^{\alpha_1} \cdots \sigma_N^{\alpha_N}$, can be written as

$$ \Xi_{P_\alpha} = p(\alpha_1)p(\alpha_2|\alpha_1)p(\alpha_3|\alpha_1, \alpha_2)\cdots, \quad (59) $$

where $p(\alpha_j|\alpha_1, \dots, \alpha_{j-1})$ is the conditional probability that the Pauli matrix $j$ is $\sigma_j^{\alpha_j}$, if the Pauli matrices $k$ with $k < j$ are $\sigma_1^{\alpha_1}, \dots, \sigma_2^{\alpha_2}$, etc. Importantly, these probabilities can be efficiently computed for MPSs. To see this, we write explicitly

$$ p(\alpha_j|\alpha_1, \dots, \alpha_{j-1}) = \mathrm{Tr} \left[ \rho_{1, \dots, j} \sigma_j^{\alpha_1} \cdots \sigma_j^{\alpha_j} \times \rho_{1, \dots, j} \sigma_j^{\alpha_1} \cdots \sigma_j^{\alpha_j} \right], \quad (60) $$

where $\rho_{1, \dots, j}$ is the reduced density matrix over the first $j$ qubits. This formula can be verified making repeated use of the identity $(1/2) \sum_\alpha \sigma^\alpha(\cdot)\sigma^\alpha = 1 \mathrm{Tr}(\cdot)$. For an MPS with bond dimension $\chi$, we can use the standard prescriptions [17, 48] to evaluate the optimal contraction cost of the tensor network associated

with (60), which yields $O(\chi^3 N)$. Therefore, the conditional probabilities can be computed efficiently in $N$. Note that the dependence with $\chi$ 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 $\alpha$, according to a probability which is exactly $\Xi_{P_\alpha}$. Then, we can evaluate $M_1(|\Psi_N\rangle)$ by averaging the value $\log \Xi_{P_\alpha}(|\Psi_N\rangle)$ 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 $\varepsilon_S$ the difference between the average value, obtained by sampling the distribution $\Xi_P$, and the actual value of $M_1$, standard arguments yield $\varepsilon_S = O(1/\sqrt{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 ground-state of the interacting XXZ Heisenberg Hamiltonian

$$ H_{XXZ} = - \sum_{k=1}^{N-1} (\sigma_k^x \sigma_{k+1}^x + \sigma_k^y \sigma_{k+1}^y + \Delta \sigma_k^z \sigma_{k+1}^z), \quad (61) $$

where $\Delta$ is the anisotropy parameter, while we choose open boundary conditions. Note that the total magnetization $C = \sum_{k=1}^N \sigma_k^z$ is conserved, i.e. $[H_{XXZ}, C] = 0$.

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|\psi\rangle = 0$. In Fig. 2a, we plot the SE density $m_n = M_n/N$ against $\Delta$ 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 $\Delta$, where $m$ is first increasing with $\Delta$, 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 $\Delta \to 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 $\Delta = 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 $\chi$ 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(|\psi\rangle)$ up to fixed accuracy could scale exponentially in the number of qubits $N$.
<figure><img src="image_2.png" /><figcaption>Figure 2: a) Von Neumann SE density $m_1$ and Rényi-2 SE density against anisotropy $\Delta$ for the ground state of the XXZ model, at half-filling. For $m_1$, the bond dimension is $\chi = 40$, while the number of samples is $S = 10^5$. For $m_2$, we used the method of Ref. [25] with $\chi = 12$. b) Von Neumann SE density $m_1$ against anisotropy $\Delta$ for the ground state of the XXZ model at half-filling. Here, $\chi = 40$ and $S = 10^5$.</figcaption></figure>

# 6 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 \le 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 min-relative 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 \ge 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 \le 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.

# Acknowledgements

We thank Hyukjoon Kwon, Giacomo de Palma, and especially Ludovico Lami for inspiring discussions.

# References

[1] Scott Aaronson and Daniel Gottesman. Improved simulation of stabilizer circuits. *Phys. Rev. A*, 70:052328, Nov 2004. DOI: 10.1103/PhysRevA.70.052328. URL https://link.aps.org/doi/10.1103/PhysRevA.70.052328.

[2] Hussain Anwar, Benjamin J Brown, Earl T Campbell, and Dan E Browne. Fast decoders for qudit topological codes. *New J. Phys.*, 16(6):063038, 2014. DOI: 10.1088/1367-2630/16/6/063038. URL https://iopscience.iop.org/article/10.1088/1367-2630/16/6/063038/meta.

[3] Michael Beverland, Earl Campbell, Mark Howard, and Vadym Kliuchnikov. Lower bounds on the non-clifford resources for quantum computations. *Quantum Science Tech.*, 5(3):035009, 2020. DOI: 10.1088/2058-9565/ab8963. URL https://doi.org/10.1088/2058-9565/ab8963.

[4] Sergey Bravyi and Alexei Kitaev. Universal quantum computation with ideal clifford gates and noisy ancillas. *Phys. Rev. A*, 71:022316, Feb 2005. DOI: 10.1103/PhysRevA.71.022316. URL https://link.aps.org/doi/10.1103/PhysRevA.71.022316.
[5] Sergey Bravyi, Dan Browne, Padraic Calpin, Earl Campbell, David Gosset, and Mark Howard. Simulation of quantum circuits by low-rank stabilizer decompositions. *Quantum*, 3:181, 2019. DOI: 10.22331/q-2019-09-02-181. URL https://quantum-journal.org/papers/q-2019-09-02-181/.

[6] Kaifeng Bu and Dax Enshan Koh. Efficient classical simulation of clifford circuits with nonstabilizer input states. *Phys. Rev. Lett.*, 123:170502, Oct 2019. DOI: 10.1103/PhysRevLett.123.170502. URL https://link.aps.org/doi/10.1103/PhysRevLett.123.170502.

[7] Kaifeng Bu, Roy J Garcia, Arthur Jaffe, Dax Enshan Koh, and Lu Li. Complexity of quantum circuits via sensitivity, magic, and coherence. *arXiv:2204.12051*, 2022. URL https://doi.org/10.48550/arXiv.2204.12051.

[8] Earl T. Campbell. Catalysis and activation of magic states in fault-tolerant architectures. *Phys. Rev. A*, 83:032317, Mar 2011. DOI: 10.1103/PhysRevA.83.032317. URL https://link.aps.org/doi/10.1103/PhysRevA.83.032317.

[9] Earl T. Campbell. Enhanced fault-tolerant quantum computing in d-level systems. *Phys. Rev. Lett.*, 113:230501, Dec 2014. DOI: 10.1103/PhysRevLett.113.230501. URL https://link.aps.org/doi/10.1103/PhysRevLett.113.230501.

[10] Earl T. Campbell, Hussain Anwar, and Dan E. Browne. Magic-state distillation in all prime dimensions using quantum reed-muller codes. *Phys. Rev. X*, 2:041021, Dec 2012. DOI: 10.1103/PhysRevX.2.041021. URL https://link.aps.org/doi/10.1103/PhysRevX.2.041021.

[11] Earl T Campbell, Barbara M Terhal, and Christophe Vuillot. Roads towards fault-tolerant universal quantum computation. *Nature*, 549(7671):172-179, 2017. DOI: 10.1038/nature23460. URL https://www.nature.com/articles/nature23460.

[12] Liyuan Chen, Roy J Garcia, Kaifeng Bu, and Arthur Jaffe. Magic of random matrix product states. *arXiv:2211.10350*, 2022. URL https://doi.org/10.48550/arXiv.2211.10350.

[13] Eric Chitambar and Gilad Gour. Quantum resource theories. *Rev. Mod. Phys.*, 91: 025001, Apr 2019. DOI: 10.1103/RevModPhys.91.025001. URL https://link.aps.org/doi/10.1103/RevModPhys.91.025001.

[14] J Ignacio Cirac, David Perez-Garcia, Norbert Schuch, and Frank Verstraete. Matrix product density operators: Renormalization fixed points and boundary theories. *Ann. Phys.*, 378:100-149, 2017. DOI: 10.1016/j.aop.2016.12.030. URL http://www.sciencedirect.com/science/article/pii/S0003491616303013.

[15] J Ignacio Cirac, David Perez-Garcia, Norbert Schuch, and Frank Verstraete. Matrix product states and projected entangled pair states: Concepts, symmetries, theorems. *Rev. Mod. Phys.*, 93(4):045003, 2021. DOI: 10.1103/RevModPhys.93.045003. URL https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.93.045003.

[16] Bryan Eastin and Emanuel Knill. Restrictions on transversal encoded quantum gate sets. *Phys. Rev. Lett.*, 102:110502, Mar 2009. DOI: 10.1103/PhysRevLett.102.110502. URL https://link.aps.org/doi/10.1103/PhysRevLett.102.110502.

[17] Glen Evenbly. A practical guide to the numerical implementation of tensor networks: Contractions, decompositions, and gauge freedom. *Frontiers in Applied Mathematics and Statistics*, 8:806549, 2022. DOI: 10.3389/fams.2022.806549. URL https://doi.org/10.3389/fams.2022.806549.

[18] Andrew J. Ferris and Guifre Vidal. Perfect sampling with unitary tensor networks. *Phys. Rev. B*, 85:165146, Apr 2012. DOI: 10.1103/PhysRevB.85.165146. URL https://link.aps.org/doi/10.1103/PhysRevB.85.165146.

[19] Daniel Gottesman. Stabilizer codes and quantum error correction. *Caltech Ph. D. PhD thesis*, Thesis, eprint: quant-ph/9705052, 1997. URL https://doi.org/10.48550/arXiv.quant-ph/9705052.

[20] Daniel Gottesman. The heisenberg representation of quantum computers. *arXiv quant-ph/9807006*, 1998. URL https://doi.org/10.48550/arXiv.quant-ph/9807006.

[21] Daniel Gottesman. Theory of fault-tolerant quantum computation. *Phys. Rev. A*, 57: 127-137, Jan 1998. DOI: 10.1103/PhysRevA.57.127. URL https://link.aps.org/doi/10.1103/PhysRevA.57.127.

[22] Oliver Hahn, Alessandro Ferraro, Lina Hultquist, Giulia Ferrini, and Laura García-Álvarez. Quantifying qubit magic resource with gottesman-kitaev-preskill encoding. *Phys. Rev. Lett.*, 128:210502, May 2022. DOI: 10.1103/PhysRevLett.128.210502. URL https://link.aps.org/doi/10.1103/PhysRevLett.128.210502.

[23] Oliver Hahn, Alessandro Ferraro, Lina Hultquist, Giulia Ferrini, and Laura García-Álvarez. Eratrum: Quantifying qubit magic resource with gottesman-kitaev-preskill encoding [phys. rev. lett. 128, 210502 (2022)]. *Phys. Rev. Lett.*, 131:049901, Jul 2023. DOI: 10.1103/PhysRevLett.131.049901. URL https://link.aps.org/doi/10.1103/PhysRevLett.131.049901.

[24] Tobias Haug and M.S. Kim. Scalable measures of magic resource for quantum computers. *PRX Quantum*, 4:010301, Jan
2023. DOI: 10.1103/PRXQuantum.4.010301.
URL https://link.aps.org/doi/10.1103/PRXQuantum.4.010301.

[25] Tobias Haug and Lorenzo Piroli. Quantifying nonstabilizerness of matrix product states. Phys. Rev. B, 107(3):035148, 2023. DOI: 10.1103/Phys-RevB.107.035148. URL https://link.aps.org/doi/10.1103/PhysRevB.107.035148.

[26] Tobias Haug, Soovin Lee, and MS Kim. Efficient stabilizer entropies for quantum computers. arXiv:2305.19152, 2023. URL https://doi.org/10.48550/arXiv.2305.19152.

[27] Arne Heimendahl, Markus Heinrich, and David Gross. The axiomatic and the operational approaches to resource theories of magic do not coincide. J. Math. Phys., 63 (11):112201, 2022. DOI: 10.1063/5.0085774. URL https://aip.scitation.org/doi/full/10.1063/5.0085774.

[28] Mark Howard and Earl Campbell. Application of a resource theory for magic states to fault-tolerant quantum computing. Phys. Rev. Lett., 118:090501, Mar 2017. DOI: 10.1103/Phys-RevLett.118.090501. URL https://link.aps.org/doi/10.1103/PhysRevLett.118.090501.

[29] Jiaqing Jiang and Xin Wang. Lower bound for the t count via unitary stabilizer nullity. Physical Review Applied, 19(3):034052, 2023. DOI: 10.1103/Phys-RevApplied.19.034052. URL https://journals.aps.org/prapplied/abstract/10.1103/PhysRevApplied.19.034052.

[30] A Yu Kitaev. Fault-tolerant quantum computation by anyons. Ann. Phys., 303(1):2-30, 2003. DOI: 10.1016/S0003-4916(02)00018-0. URL http://www.sciencedirect.com/science/article/pii/S0003491602000180.

[31] Guglielmo Lami and Mario Collura. Quantum magic via perfect sampling of matrix product states. arXiv:2303.05536, 2023. URL https://doi.org/10.48550/arXiv.2303.05536.

[32] Lorenzo Leone, Salvatore F. E. Oliviero, and Alioscia Hamma. Stabilizer rényi entropy. Phys. Rev. Lett., 128:050402, Feb 2022. DOI: 10.1103/PhysRevLett.128.050402. URL https://link.aps.org/doi/10.1103/PhysRevLett.128.050402.

[33] Lorenzo Leone, Salvatore F. E. Oliviero, Gianluca Esposito, and Alioscia Hamma. Phase transition in stabilizer entropy and efficient purity estimation. arXiv:2302.07895, 2023. URL https://doi.org/10.48550/arXiv.2302.07895.

[34] Lorenzo Leone, Salvatore F. E. Oliviero, and Alioscia Hamma. Nonstabilizerness determining the hardness of direct fidelity estimation. Phys. Rev. A, 107:022429, Feb 2023. DOI: 10.1103/PhysRevA.107.022429.

URL https://link.aps.org/doi/10.1103/PhysRevA.107.022429.

[35] Zi-Wen Liu and Andreas Winter. Many-body quantum magic. PRX Quantum, 3: 020333, May 2022. DOI: 10.1103/PRXQuantum.3.020333. URL https://link.aps.org/doi/10.1103/PRXQuantum.3.020333.

[36] Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge University Press, 2011. ISBN 9781107002173. DOI: doi:10.1017/CBO9780511976667.

[37] J Odavić, T Haug, G Torre, A Hamma, F Franchini, and SM Giampaolo. Complexity of frustration: a new source of non-local non-stabilizerness. arXiv:2209.10541, 2022. URL https://doi.org/10.48550/arXiv.2209.10541.

[38] Salvatore F. E. Oliviero, Lorenzo Leone, and Alioscia Hamma. Magic-state resource theory for the ground state of the transverse-field ising model. Phys. Rev. A, 106: 042426, Oct 2022. DOI: 10.1103/Phys-RevA.106.042426. URL https://link.aps.org/doi/10.1103/PhysRevA.106.042426.

[39] Salvatore F. E. Oliviero, Lorenzo Leone, Alioscia Hamma, and Seth Lloyd. Measuring magic on a quantum processor. npj Quantum Information, 8(1):148, 2022. DOI: 10.1038/s41534-022-00666-5. URL https://www.nature.com/articles/s41534-022-00666-5.

[40] Salvatore FE Oliviero, Lorenzo Leone, Seth Lloyd, and Alioscia Hamma. Black hole complexity, unscrambling, and stabilizer thermal machines. arXiv:2212.11337, 2022. URL https://doi.org/10.48550/arXiv.2212.11337.

[41] D Perez-Garcia, F Verstraete, MM Wolf, and JI Cirac. Matrix product state representations. Quantum Inf. Comp., 7(5):401-430, 2007. DOI: 10.26421/QIC7.5-6-1. URL https://arxiv.org/abs/quant-ph/0608197.

[42] John Preskill. Fault-tolerant quantum computation. In Introduction to quantum computation and information, pages 213-269. World Scientific, 1998. DOI: 10.1142/9789812385253'0008. URL https://www.worldscientific.com/doi/10.1142/9789812385253_0008.

[43] Bartosz Regula. Convex geometry of quantum resource quantification. J. Phys. A: Math. Theor., 51(4):045303, 2017. DOI: 10.1088/1751-8121/aa9100. URL https://iopscience.iop.org/article/10.1088/1751-8121/aa9100.

[44] Saubhik Sarkar, Chiranjib Mukhopadhyay, and Abolfazl Bayat. Characterization of an operational quantum resource in a critical many-body system. New J. Phys., 22(8): 083077, 2020. DOI: 10.1088/1367-2630/aba919.
URL https://iopscience.iop.org/article/10.1088/1367-2630/aba919.

[45] James R. Seddon, Bartosz Regula, Hakop Pashayan, Yingkai Ouyang, and Earl T. Campbell. Quantifying quantum speedups: Improved classical simulation from tighter magic monotones. PRX Quantum, 2:010345, Mar 2021. DOI: 10.1103/PRXQuantum.2.010345. URL https://link.aps.org/doi/10.1103/PRXQuantum.2.010345.

[46] Troy J Sewell and Christopher David White. Mana and thermalization: Probing the feasibility of near-clifford hamiltonian simulation. Phys. Rev. B, 106(12):125130, 2022. DOI: 10.1103/PhysRevB.106.125130. URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.106.125130.

[47] Peter W Shor. Fault-tolerant quantum computation. In Proceedings of 37th conference on foundations of computer science, pages 56–65. IEEE, 1996. DOI: 10.1109/SFCS.1996.548464. URL https://ieeexplore.ieee.org/document/548464.

[48] Pietro Silvi, Ferdinand Tschirsich, Matthias Gerster, Johannes Jünemann, Daniel Jaschke, Matteo Rizzi, and Simone Montangero. The tensor networks anthology: Simulation techniques for many-body quantum lattice systems. SciPost Phys. Lect. Notes, page 8, 2019. DOI: 10.21468/SciPostPhysLectNotes.8. URL https://scipost.org/10.21468/SciPostPhysLectNotes.8.

[49] Emanuele Tirrito, Poetri Sonya Tarabunga, Gugliemo Lami, Titas Chanda, Lorenzo Leone, Salvatore FE Oliviero, Marcello Dalmonte, Mario Collura, and Alioscia Hamma. Quantifying non-stabilizerness through entanglement spectrum flatness. arXiv:2304.01175, 2023. URL https://doi.org/10.48550/arXiv.2304.01175.

[50] Xhek Turkeshi, Marco Schirò, and Piotr Sierant. Measuring magic via multifractal flatness. arXiv:2305.11797, 2023. URL https://doi.org/10.48550/arXiv.2305.11797.

[51] Victor Veitch, SA Hamed Mousavian, Daniel Gottesman, and Joseph Emerson. The resource theory of stabilizer quantum computation. New J. Phys., 16(1):013009, 2014. DOI: 10.1088/1367-2630/16/1/013009. URL https://doi.org/10.1088/1367-2630/16/1/013009.

[52] Xin Wang, Mark M Wilde, and Yuan Su. Quantifying the magic of quantum channels. New J. Phys., 21(10):103002, 2019. DOI: 10.1088/1367-2630/ab451d. URL https://doi.org/10.1088/1367-2630/ab451d.

[53] Christopher David White, ChunJun Cao, and Brian Swingle. Conformal field theories are magical. Phys. Rev. B, 103: 075145, Feb 2021. DOI: 10.1103/Phys-

RevB.103.075145. URL https://link.aps.org/doi/10.1103/PhysRevB.103.075145.

## A Details on the violation of monotonicity

In this appendix, we provide further details on the counterexamples presented in Sec. 3. First, we explain how we arrived at the state (13).

We initially looked for violations of the strong monotonicity condition (3). To this end, we defined the functional

$$ \Delta_n(|\psi\rangle) = M_n(|\psi\rangle) - \sum_{a=0}^{1} p_a M_n(|\psi_a\rangle), \quad (62) $$

where

$$ |\psi_a\rangle = \frac{1}{\sqrt{p_a}} \Pi_a |\psi\rangle, \quad (63) $$

and $\Pi_a = |a\rangle\langle a| \otimes 1_{2,3,...,N}$, while $p_a = \langle\psi|\Pi_a|\psi\rangle$. Namely, $|\psi_a\rangle$ are the states obtained after a computational-basis measurement of qubit 1, and $p_a$ are the corresponding probabilities. Note that a negative value of $\Delta_n(|\psi\rangle)$ implies a violation of the strong monotonicity condition.

Next, using a gradient-descent method, we have looked numerically for the minimum of (62) for different values of $N$. For $N \le 3$, we found $\Delta_n(|\psi\rangle) \ge 0$, while for $N = 4$ we found that there are states for which $\Delta_n(|\psi\rangle) < 0$ for $0 \le n < 2$.

The state (13) was obtained as the result of this minimization procedure for $N = 4$ (decimal numbers were rounded as fractions, e.g. $0.408248 \approx 1/\sqrt{6}$, by guesswork). Remarkably, by analytical inspection, we found that this state also has additional properties. Most notably, it is an eigenstate of the Clifford operator (16). This made it possible to devise the Clifford protocol presented in Sec. 3.1, which was done by analytic inspection and guesswork. Therefore, although the state was found by looking for violations of the strong-monotonicity condition, it turned out to also provide a counterexample to the monotonicity condition (1).

The same method applied to $N = 5$ qubits becomes cumbersome. While we were able to carry out the minimization procedure numerically, the states which were obtained in this way did not have the desirable properties of (13). In particular, using this method we did not find states which were eigenstates of Clifford operators.

Next, we provide further details on the state $|\psi^\epsilon\rangle$ of (30). We show in particular that its SE for $n > 1$ is bounded, as anticipated in the main text.

We start by defining

$$ G_n(|\psi^\epsilon\rangle) = \sum_{P \in \mathcal{P}_N} \frac{\langle\psi^\epsilon|P|\psi^\epsilon\rangle^{2n}}{2^N}, \quad (64) $$
so that $M_n(|\psi^\epsilon\rangle) = \log G_n(|\psi^\epsilon\rangle)/(1-n)$. Next, we

note that the state $|\psi^\epsilon\rangle$ is invariant under permutation of qubits. Therefore, we can simplify

$$G_n(|\psi^\epsilon\rangle) = \frac{1}{2^N} \sum_{N_z=0}^{N} \sum_{N_x=0}^{N-N_z} \sum_{N_y=0}^{N-N_z-N_x} \left( N_z, N_x, N_y, N-N_z-N_x-N_y \right) \langle \psi^\epsilon | \prod_{i=1}^{N_z} \sigma_i^z \prod_{j=N_z+1}^{N_z+N_x} \sigma_j^x \prod_{k=N_z+N_x+1}^{N_z+N_x+N_y} \sigma_k^y | \psi^\epsilon \rangle^{2n} \quad (65)$$

In addition, the expectation value of Pauli matrices can be evaluated explicitly, yielding

$$\langle \psi^\epsilon | \prod_{i=1}^{N_z} \sigma_i^z \prod_{j=N_z+1}^{N_z+N_x} \sigma_j^x \prod_{k=N_z+N_x+1}^{N_z+N_x+N_y} \sigma_k^y | \psi^\epsilon \rangle = \frac{\delta_{N_x,0} \delta_{N_y,0} + \epsilon^2 3^{-(N_x+N_y+N_z)/2}}{N_\epsilon} + \frac{2\epsilon \sin(\beta)^{N_x+N_y} \cos(\beta)^{N-N_x-N_y} \cos[\pi(N-N_x+N_y)/4]}{N_\epsilon} \quad (66)$$

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 \sim 50$ qubits. We have done this for different values of $n > 1$, and confirmed that $M_n(|\psi^\epsilon\rangle)$ 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(|\psi^\epsilon\rangle) \approx (I)/N_\epsilon^{2n} + (II)/N_\epsilon^{2n}$, where

$$(I) = \frac{1}{2^N} \sum_{N_z \ge 0} \binom{N}{N_z} \left[ 1 + \epsilon^2 3^{-N_z/2} + 2\epsilon \cos(\beta)^N \cos(\pi N/4) \right]^{2n}, \quad (67)$$

and

$$(II) = \frac{1}{2^N} \sum_{\substack{N_z \ge 0 \\ N_x \parallel N_y \ne 0}} \binom{N}{N_z, N_x, N_y, N-N_x-N_y-N_z} \times (\epsilon^2 3^{-(N_x+N_y+N_z)/2} + 2\epsilon \sin(\beta)^{N_x+N_y} \times \cos(\beta)^{N-N_x-N_y} \cos[\pi(N-N_x+N_y)/4])^{2n}. \quad (68)$$

First, we note that

$$(I) \ge \frac{1}{2^N} \sum_{N_z \ge 0} \binom{N}{N_z} (1 - \delta_N)^{2n} = (1 - \delta_N)^{2n}, \quad (69)$$

where $\delta_N = |2\epsilon \cos(\beta)^N \cos(\pi 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 ($|A_1 + A_2|^{2n} \le 2^{2n-1}(|A_1|^{2n} + |A_2|^{2n}|$), we obtain

$$(II) \le \frac{2^{2n-1}}{2^N} \sum_{N_x, N_y, N_z} \binom{N}{N_z, N_x, N_y, N-N_x-N_y-N_z} \times \left\{ (\epsilon^2 3^{-(N_x+N_y+N_z)/2})^{2n} + [2\epsilon \sin(\beta)^{N_x+N_y} \cos(\beta)^{N-N_x-N_y}]^{2n} \right\}, \quad (70)$$

where we have extended the sum to $N_x, N_y, N_z \ge 0$

<figure><img src="image_3.png" /><figcaption>Figure 3: Numerical test of the asymptotic prediction (72) (for $n=2$). Symbols correspond to the exact evaluation of (65), while the dashed line is the prediction (72). Here we chose $\epsilon = 0.5$.</figcaption></figure>

such that $N_x + N_y + N_z = N$. Applying the multivariate theorem, we obtain

$$|(II)| \le 2^{2n-1} \frac{\epsilon^{4n}}{2^N} \left(1 + 3 \left(\frac{1}{\sqrt{3}}\right)^{2n}\right)^N + 2^{2n-1} \frac{(2\epsilon)^{2n}}{2^N} \left(2 \cos(\beta)^{2n} + 2 \sin(\beta)^{2n}\right)^N. \quad (71)$$

For $n > 1$ we see that (II) is exponentially small in $N$. Therefore, using $\log N_\epsilon \simeq \log(1 + \epsilon^2)$ (up to exponentially small terms in $N$) we have

$$M_n(|\psi^\epsilon\rangle) \le \frac{2n}{n-1} \log(1 + \epsilon^2) + \tilde{\delta}_N, \quad (72)$$

where $\tilde{\delta}_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 re-
sults discussed in Sec. 5. First, we present addi-
tional numerical data about the perfect-MPS sam-
pling 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
|*m*̂₁ - *m*₁|, where *m*̂₁ is an estimate for *m*₁, com-
puted by averaging over a given number of samples
*S*, while *m*₁ is our best prediction, computed by aver-
aging over a very large value of samples, *S*<sub>max</sub> (in this
case, *S*<sub>max</sub> = 10⁵). The plot clearly shows the scaling
|*m*̂₁ - *m*₁| ~ *S*⁻¹/², as expected.

<figure><img src="image_4.png" /><figcaption>Figure 4: Average error in estimation of the von Neumann SE density |m̂₁ − m₁| against number of samples *S*. m̂₁ is the estimated SE density, while m₁ is a highly accurate estimation from *S* = 10⁵ 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*.</figcaption></figure>

Next, we study the dependence of the estimated val-
ues of m₁ and m₂ 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₂/N = m₂ as function of
bond dimension χ, which was computed using the
method of Ref. [25]. We find that the value of m₂
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 approxima-
tion with bond dimension χ, |ΨN(χ)⟩, and the true
ground state |GS⟩. The true ground-state |GS⟩ is esti-
mated by taking a very large bond dimension (in this
case, χ ≃ 200).
a

<figure><img src="image_5.png" /></figure>

b

<figure><img src="image_6.png" /></figure>

c

<figure><img src="image_7.png" /><figcaption>Figure 5: a) Rényi-2 SE density $m_2$ against bond dimension $\chi$. The plot corresponds to the MPS representing the ground state of the Heisenberg model at half-filling. Here, we have used the method of Ref. [25]. b) Estimated value of $m_1$ as a function of $\chi$, using the method presented in Sec. 5. c) Fidelity $F = |\langle\Psi_N(\chi)|GS\rangle|^2$ between the MPS approximation for the ground state, $|\Psi_N(\chi)\rangle$, and the true ground state $|GS\rangle$ as function of $\chi$. The true ground state $|GS\rangle$ is estimated by taking a very large bond dimension (in this case, $\chi \approx 200$). In all plots, we set $\Delta = 0.95$.</figcaption></figure>