Quantifying magic for multi-qubit operations

TL;DR

This paper develops a framework to quantify the non-stabiliserness of quantum operations, introducing new monotones and algorithms that improve classical simulation efficiency for multi-qubit quantum circuits.

Contribution

It extends the concept of magic monotones from states to channels, introduces channel robustness and magic capacity, and provides more efficient simulation algorithms.

Findings

New monotones for quantum channels are defined.

Simulation algorithms outperform previous methods in certain cases.

Techniques for easier calculation of monotones are proposed.

Abstract

The development of a framework for quantifying "non-stabiliserness" of quantum operations is motivated by the magic state model of fault-tolerant quantum computation, and by the need to estimate classical simulation cost for noisy intermediate-scale quantum (NISQ) devices. The robustness of magic was recently proposed as a well-behaved magic monotone for multi-qubit states and quantifies the simulation overhead of circuits composed of Clifford+T gates, or circuits using other gates from the Clifford hierarchy. Here we present a general theory of the "non-stabiliserness" of quantum operations rather than states, which are useful for classical simulation of more general circuits. We introduce two magic monotones, called channel robustness and magic capacity, which are well-defined for general n-qubit channels and treat all stabiliser-preserving CPTP maps as free operations. We present two complementary Monte Carlo-type classical simulation algorithms with sample complexity given by these quantities and provide examples of channels where the complexity of our algorithms is exponentially better than previous known simulators. We present additional techniques that ease the difficulty of calculating our monotones for special classes of channels.

Tables

Table 1. Table 1: Tightness of bound given by Theorem 3.2 , as determined by numerical estimation of diagonal gates, where ℛ Φ subscript ℛ Φ \mathcal{R}_{\Phi} is the robustness of the Choi state, 𝒞 𝒞 \mathcal{C} is the magic capacity, ℛ ∗ subscript ℛ \mathcal{R}_{*} is the trace-preserving variant of ℛ Φ subscript ℛ Φ \mathcal{R}_{\Phi} . Here an equality indicates that in all cases investigated, values calculated were equal up to the precision of the solver. Multicontrol phase gates are taken to be those represented by unitaries of the form diag ​ ( 1 , … , 1 , exp ⁡ [ i ​ π / 2 t ] ) diag 1 … 1 𝑖 𝜋 superscript 2 𝑡 \mathrm{diag}(1,\ldots,1,\exp[i\pi/2^{t}]) .

$n$	2	3	4	5
Multicontrol gates, $t=0$	${ℛ}_{\Phi }=𝒞={ℛ}_{\ast }$	${ℛ}_{\Phi }=𝒞={ℛ}_{\ast }$	${ℛ}_{\Phi }=𝒞={ℛ}_{\ast }$
Multicontrol gates, $t\ge 1$	${ℛ}_{\Phi }=𝒞={ℛ}_{\ast }$	${ℛ}_{\Phi }=𝒞={ℛ}_{\ast }$
Random phase gates	${ℛ}_{\Phi }=𝒞={ℛ}_{\ast }$	${ℛ}_{\Phi }=𝒞\le {ℛ}_{\ast }$	${ℛ}_{\Phi }=𝒞\le {ℛ}_{\ast }$	-

Table 2. Table 2: Number of pure stabiliser states N S subscript 𝑁 𝑆 N_{S} for number of qubits n 𝑛 n .

$n$	$N_S$
1	6
2	60
3	1,080
4	36,720
5	2,423,520
6	315,057,600

Table 3. Table 3: Number of n 𝑛 n -qubit stabiliser states compared with number of affine spaces. By trivial affine spaces we mean those comprised of a single element, which correspond to computational basis states. Diagonal CPTP channels act as the identity on such states.

$n$	Stabiliser states	Total affine spaces	Non-trivial affine spaces
2	60	11	7
3	1,080	51	43
4	36,720	307	291
5	2,423,520	2451	2419

Table 4. Table 4: Final robustness after multicontrol- T 𝑇 T gate applied to input stabiliser states | ℒ ⟩ ket ℒ \ket{\mathcal{L}} with k = dim ( ℒ ) 𝑘 dimension ℒ k=\dim(\mathcal{L}) . In each column, the maximum robustness (i.e. the capacity) is highlighted red.

Linear subspace	Number of qubits, $n$
dimension, $k$	2	3	4	5
1	$1.414$	$1.414$	$1.414$	$1.414$
2	$1.849$	$1.849$	$1.849$	$1.849$
3	-	$2.195$	$2.195$	$2.195$
4	-	-	$2.264$	$2.264$
5	-	-	-	$2.195$

Equations

$$\mathcal{R}\text{quantity}(\rho )=\underset{\vec{q}}{\min }\text{quantity}\Vert \vec{q}{\Vert }_1:\sum _j{q}_j\text{outerproduct}{\varphi }_j{\varphi }_j=\rho ,|{\varphi }_j⟩\in {\mathrm{STAB}}_n.$$

$$\mathcal{R}\text{quantity}(\rho )=\underset{{\rho }_{\pm }\in {\mathrm{STAB}}_n}{\min }\text{quantity}1+2p:(1+p){\rho }_{+}-p{\rho }_{-}=\rho ,p\ge 0.$$

$${\Phi }_{\mathcal{E}}^{AB}=\text{quantity}({\mathcal{E}}^A\otimes \text{mathds}{1}_n)\text{outerproduct}{\Omega }_n{{\Omega }_n}^{A|B},\text{where}{|{\Omega }_n⟩}^{A|B}=\frac{1}{\sqrt{2^n}}\sum _{j=0}^{2^n-1}{|j⟩}^A\otimes {|j⟩}^B,$$

$${\mathcal{R}}_{\ast }(\mathcal{E})=\underset{{\Lambda }_{\pm }\in {\mathrm{SP}}_{n,n}\cap \mathrm{CPTP}}{\min }\text{quantity}2p+1:(1+p){\Lambda }_{+}-p{\Lambda }_{-}=\mathcal{E},p\ge 0,$$

$${\mathcal{R}}_{\ast }(\mathcal{E})=\underset{{\rho }_{\pm }\in {\mathrm{STAB}}_{2n}}{\min }\text{quantity}2p+1:(1+p){\rho }_{+}-p{\rho }_{-}={\Phi }_{\mathcal{E}}^{AB},p\ge 0,{\text{Tr}}_A({\rho }_{\pm })=\frac{\text{mathds}{1}_n}{2^n}.$$

$${\mathcal{R}}_{\ast }(\Lambda \circ \mathcal{E})\le {\mathcal{R}}_{\ast }(\mathcal{E}){\mathcal{R}}_{\ast }(\Lambda )={\mathcal{R}}_{\ast }(\mathcal{E}),$$

$$\mathcal{C}(\mathcal{E})=\underset{|\varphi ⟩\in {\mathrm{STAB}}_{2n}}{\max }\mathcal{R}\text{quantity}[\text{quantity}(\mathcal{E}\otimes \text{mathds}{1}_n)\text{outerproduct}\varphi \varphi ].$$

$$\mathcal{R}\text{quantity}({\Phi }_{\mathcal{E}})\le \mathcal{C}\text{quantity}(\mathcal{E})\le {\mathcal{R}}_{\ast }\text{quantity}(\mathcal{E}).$$

$$\mathcal{R}\text{quantity}({\Phi }_{\mathcal{U}})=\mathcal{C}\text{quantity}(\mathcal{U})={\mathcal{R}}_{\ast }\text{quantity}(\mathcal{U}).$$

$$\text{quantity}({\mathcal{E}}_T\otimes \text{mathds}1)\text{outerproduct}{\Phi }_{+}{\Phi }_{+}=\frac{1}{2}\text{quantity}(\text{outerproduct}0T^{\ast }0T^{\ast }+\text{outerproduct}1T_{\perp }^{\ast }1T_{\perp }^{\ast }),$$

$$\mathcal{R}\text{quantity}[\text{quantity}({\mathcal{E}}^A\otimes \text{mathds}{1}_{n+m})\text{outerproduct}\varphi {\varphi }^{AB}]=\mathcal{R}\text{quantity}[\text{quantity}({\mathcal{E}}^A\otimes \text{mathds}{1}_n)\text{outerproduct}\psi {\psi }^{A{B}^{\prime }}].$$

$${|\varphi ⟩}^{AB}=\text{quantity}(\text{mathds}{1}_n\otimes U^B){|\psi ⟩}^{A{B}^{\prime }}{|{\psi }^{\prime }⟩}^{B^{\prime \prime }},$$

$$\mathcal{R}\text{quantity}[\text{quantity}({\mathcal{E}}^A\otimes \text{mathds}{1}_{n+m})\text{outerproduct}\varphi {\varphi }^{AB}]$$

$$=\mathcal{R}\text{quantity}[\text{quantity}({\mathcal{E}}^A\otimes \text{mathds}{1}_n)\text{outerproduct}\psi {\psi }^{A{B}^{\prime }}],$$

$${\Phi }_{\mathcal{E}}^{AB}=\text{quantity}({\mathcal{E}}^A\otimes \text{mathds}{1}^B)\text{outerproduct}{\Omega }_n{{\Omega }_n}^{A|B},\text{where}{|{\Omega }_n⟩}^{A|B}=\frac{1}{\sqrt{2^n}}\sum _{j=0}^{2^n-1}{|j⟩}^A\otimes {|j⟩}^B.$$

$$\text{Tr}[A\mathcal{E}\text{quantity}(\rho )]=2^n\text{Tr}[{\Phi }_{\mathcal{E}}(A\otimes {\rho }^T)],\forall \rho ,A.$$

$$\text{Tr}(W_n\sigma )\le 0,\forall \sigma \in {\mathrm{STAB}}_n.$$

$$\text{Tr}[(W_n\otimes \text{outerproduct}\varphi \varphi )\sigma ]$$

$$\text{Tr}[\text{quantity}(W_n\otimes \text{mathds}{1}_m)\tilde{\sigma }]$$

$${|{\Omega }_{2n}⟩}^{A{A}^{\prime }|B{B}^{\prime }}=\frac{1}{2^{2n}}\sum _{j,k}{|j⟩}^A{|k⟩}^{A^{\prime }}\otimes {|j⟩}^B{|k⟩}^{B^{\prime }}={|{\Omega }_n⟩}^{A|B}\otimes {|{\Omega }_n⟩}^{A^{\prime }|B^{\prime }},$$

$${\Phi }_{\mathcal{E}\otimes \text{mathds}{1}_n}^{A{A}^{\prime }|B{B}^{\prime }}=({\mathcal{E}}^A\otimes \text{mathds}{1}_n^{A^{\prime }})\otimes \text{mathds}{1}_{2n}^{B{B}^{\prime }}(\text{outerproduct}{\Omega }_{2n}{{\Omega }_{2n}}^{A{A}^{\prime }|B{B}^{\prime }})={\Phi }_{\mathcal{E}}^{AB}\otimes \text{outerproduct}{\Omega }_n{{\Omega }_n}^{A^{\prime }|B^{\prime }}.$$

$$0<\frac{1}{2^n}\text{Tr}[W_{{\rho }^{\prime }}{\rho }^{\prime }]$$

$${\mathcal{R}}_{\ast }(\mathcal{E})=\underset{{\Lambda }_{\pm }\in {\mathrm{SP}}_{n,n}}{\min }\text{quantity}2p+1:(1+p){\Lambda }_{+}-p{\Lambda }_{-}=\mathcal{E},p\ge 0.$$

$$\mathcal{E}(\text{outerproduct}+{+}^{\otimes n})=(1+p){\rho }_{+}-p{\rho }_{-},\text{where}{\rho }_{\pm }\in {\mathrm{STAB}}_n,\text{and}⟨x|{\rho }_{\pm }|x⟩=\frac{1}{2^n},\forall x.$$

$${\mathcal{E}}_j=(1+p_j){\Lambda }_{j,+}-p_j{\Lambda }_{j,-},$$

$${\mathcal{E}}_2\circ {\mathcal{E}}_1=(1+q){\Lambda }_{+}^{\prime }-q{\Lambda }_{-}^{\prime },$$

$${\Lambda }_{+}^{\prime }$$

$${\Lambda }_{-}^{\prime }$$

$$q$$

$$1+2q$$

$${\mathcal{R}}_{\ast }({\mathcal{E}}^A\otimes {\mathcal{E}}^{\prime B})\le {\mathcal{R}}_{\ast }({\mathcal{E}}^A\otimes \text{mathds}{1}_m^B){\mathcal{R}}_{\ast }(\text{mathds}{1}_n^A\otimes {\mathcal{E}}^{\prime B}).$$

Full text

1 Introduction

The Gottesman-Knill theorem showed that circuits comprised of stabiliser state preparations, Clifford gates, Pauli measurements, classical randomness and conditioning can be efficiently simulated by a traditional computer Gottesman (1997); Aaronson and Gottesman (2004). If a circuit involves a relatively small proportion of non-Clifford operations, simulation may be within the reach of a classical computer, albeit with a runtime overhead that is expected to scale exponentially with the amount of resource required.

An important class of devices comprises so-called near-Clifford circuits where simulation may be feasible Bennink et al. (2017); Yoganathan et al. (2019). There are two scenarios where near-Clifford circuits are relevant. As we enter the era of Noisy Intermediate Scale Quantum (NISQ) devices Preskill (2018), many experiments proposed as demonstrators of quantum advantage may be near-Clifford so it is important to rigorously understand when a classical simulation is available. Furthermore, in the NISQ regime the need for classical simulation tools for benchmarking and verification becomes more pressing. The quantification of non-stabiliser resource is also of interest in the context of the magic state model of fault-tolerant quantum computation Knill (2005); Bravyi and Kitaev (2005); Campbell et al. (2017), the second scenario. Any device intended to provide quantum advantage must involve non-stabiliser operations. In circuits employing error-correcting codes, however, it is often not possible for the code to ‘natively’ implement non-Clifford gates fault-tolerantly Campbell et al. (2017). Instead, these gates are implemented indirectly by injection of so-called magic states. These are non-stabiliser states that must be prepared using the experimentally costly process of magic state distillation Knill (2005); Bravyi and Kitaev (2005); Reichardt (2005); Bravyi and Haah (2012); Campbell et al. (2017); Jones (2013); Trout and Brown (2015); Haah et al. (2017); Hastings and Haah (2018); Krishna and Tillich (2018); Campbell and Howard (2018); Wang et al. (2018), which is comprised of Clifford-dominated circuits.

Both of these scenarios motivate the development of a resource theory Vidal and Tarrach (1999); Grudka et al. (2014); Horodecki and Oppenheim (2013); Brandão and Gour (2015); Coecke et al. (2016); Napoli et al. (2016); Stahlke (2014); Wang et al. (2018); Takagi and Zhuang (2018); Regula (2018) where the class of free operations is generated by stabiliser state preparations and rounds of stabiliser operations as described above. For the case of odd $d$-dimensional qudits this problem is largely solved by the discrete phase space formalism Gross (2006); Mari and Eisert (2012); Veitch et al. (2012, 2014); Ahmadi et al. (2018); Delfosse et al. (2017); odd dimension qudit stabiliser states are characterised by a positive discrete Wigner function. In Ref. Pashayan et al. (2015), the discrete Wigner function was cast as a quasiprobability distribution, making a direct connection between the negativity of the distribution, and the complexity of calculating expectation values via a Monte Carlo-type simulation algorithm. These techniques have recently been extended to quantify the magic of odd-dimension qudit quantum channels Wang et al. (2019). However, the discrete phase space approach cannot be applied cleanly to qubits without excluding some Clifford operations from the free operations Delfosse et al. (2015); Raussendorf et al. (2017), or losing the ability to compose representations under tensor product Raussendorf et al. (2019). To retain all multi-qubit stabiliser channels as free operations, then, we must seek alternative approaches.

Howard and Campbell Howard and Campbell (2017) introduced a scheme where density matrices are decomposed as real linear combinations of pure stabiliser state projectors. Non-stabiliser states sit outside the convex hull of the pure stabiliser states, so their decompositions necessarily contain negative terms and can again be viewed as quasiprobability distributions, with $\ell_{1}$-norm strictly larger than 1. The robustness of magic for a state, defined as the minimum $\ell_{1}$-norm over all valid decompositions, is a monotone under stabiliser operations and has several useful resource-theoretic properties. Alternative approaches include stabiliser rank methods, where the state vector is decomposed as a superposition of stabiliser states García et al. (2014); Bravyi and Gosset (2016); Bravyi et al. (2016, 2018). Exact and approximate stabiliser rank, and the associated quantity extent, are measures of magic for pure states. Here we are interested in measures naturally suited for applications to mixed states or general, noisy quantum channels. A stabiliser-based method to simulate noisy circuits by decomposition of states into Pauli operators was recently proposed in Ref. Rall et al. (2019). In this work we characterise the cost of quantum operations with respect to the resource theory of magic. Robustness of magic naturally quantifies the cost for a subclass of non-Clifford operations, namely gates from the third level of the Clifford hierarchy. It is less clear how the framework can be extended to more general quantum operations, and formalising this is one of our main aims.

In Ref. Bennink et al. (2017), Bennink et al. presented an algorithm in which completely positive trace-preserving (CPTP) maps are decomposed as quasiprobability distributions over a subset of stabiliser-preserving operations that we will call $\mathrm{CPR}$. This subset supplements the Clifford unitaries with Pauli reset channels, in which measurement of some Pauli observable is followed by a conditional Clifford correction, so as to reset a state to a particular +1 Pauli eigenstate. While Bennink et al. showed that $\mathrm{CPR}$ spans the set of CPTP maps, there is no guarantee that all stabiliser-preserving CPTP maps can be found within its convex hull. Indeed, we will see in Section 7 there exist channels that are stabiliser-preserving, but are nevertheless assigned a non-trivial cost by the algorithm of Ref. Bennink et al. (2017). The implication is that decomposition in terms of elements of $\mathrm{CPR}$ is not the best strategy for simulating general non-stabiliser operations. An obvious extension of Ref. Bennink et al. (2017) is to replace $\mathrm{CPR}$ by the full set of stabiliser-preserving CPTP maps. The technical question to be answered is then how to correctly and concisely represent this set; how can we be sure that we have captured all possible stabiliser-preserving channels? This issue is addressed in Sections 4 and 5.

In this paper we introduce two magic monotones for channels: the channel robustness $\mathcal{R}_{*}$ and the magic capacity $\mathcal{C}$. Both are closely related to the robustness of magic for states. They are well-defined for general $n$-qubit channels and treat all stabiliser-preserving CPTP maps as free operations. We will see that these monotones give the sample complexity of two classical simulation algorithms. Other magic monotones have been proposed Veitch et al. (2014); Ahmadi et al. (2018); Wang et al. (2018) but without known connections to classical simulation algorithms. Furthermore, we give several examples of channels where the simulation complexities of our approaches are exponentially faster (as a function of gate count) than other quasiprobability simulators such as the Bennink et al. simulator Bennink et al. (2017). To our knowledge, our algorithms are the first that are known to efficiently simulate all stabiliser-preserving CPTP maps, as opposed to the set of stabiliser operations generated by Clifford gates and Pauli measurements.

The paper is structured as follows. In Section 2 we review the properties of robustness of magic and give some definitions. Next, we summarise our main results in Section 3, before pinning down what we mean by stabiliser-preserving operations in Section 4. Sections 5 and 6 are chiefly concerned with proving important properties of our monotones. Two classical simulation algorithms, each related to one of our monotones, are described in Section 7. Finally, in Section 8 we calculate the numerical values of our monotones for operations on up to five qubits, using techniques developed in Appendix E. MATLAB code to calculate each of our measures is provided at the public repository Ref. Seddon (2019).

2 Preliminaries

Let $\mathrm{STAB}_{n}$ be the set of $n$-qubit stabiliser states. In an abuse of notation we will use $\ket{\phi}\in\mathrm{STAB}_{n}$ to mean a pure state from this set, and $\rho\in\mathrm{STAB}_{n}$ to mean the density matrix of a state taken from the stabiliser polytope, the convex hull of pure stabiliser states. The pure states in $\mathrm{STAB}_{n}$ form an overcomplete basis for the set of $2^{n}$-dimensional density matrices $\mathcal{D}_{n}$. We can therefore write the density matrix for any state as an affine combination of pure stabiliser state projectors $\rho=\sum_{j}q_{j}\outerproduct{\phi_{j}}{\phi_{j}}$ where $\ket{\phi_{j}}\in\mathrm{STAB}_{n}$, and $\sum_{j}q_{j}=1$. In general, $q_{j}$ can be negative. The robustness of magic is defined as the minimal $\ell_{1}$-norm $\|\vec{q}\|_{1}=\sum_{j}\absolutevalue{q_{j}}$ over all possible decompositions Howard and Campbell (2017):

[TABLE]

In the definition above, the state of interest is expressed as a decomposition over pure stabiliser states. By collecting together all terms of the same sign, any state can instead be expressed in terms of a pair of mixed stabiliser states (Figure 1). An equivalent definition is then:

[TABLE]

The robustness of magic is a well-behaved magic monotone, having the following properties:

Convexity: $\mathcal{R}\quantity(\sum_{j}q_{j}\rho_{j})\leq\sum_{j}\absolutevalue{q_{j}}\mathcal{R}\quantity(\rho_{j})$; 2. 2.

Faithfulness: If $\rho\in\mathrm{STAB}_{n}$, then $\mathcal{R}\quantity(\rho)=1$. Otherwise $\mathcal{R}\quantity(\rho)>1$; 3. 3.

Monotonicity under stabiliser operations: If $\Lambda$ is a CPTP stabiliser-preserving operation, then $\mathcal{R}\quantity(\Lambda\quantity(\rho))\leq\mathcal{R}\quantity(\rho)$; 4. 4.

Submultiplicativity under tensor product: $\mathcal{R}\quantity(\rho^{A}\otimes\rho^{\prime B})\leq\mathcal{R}\quantity(\rho^{A})\mathcal{R}\quantity(\rho^{\prime B})$.

The quantity $\mathcal{R}$ also has a clear operational meaning, quantifying the classical simulation cost in a Monte Carlo-type scheme that samples from a quasiprobability distribution over stabiliser states Pashayan et al. (2015); Howard and Campbell (2017); Bennink et al. (2017). These algorithms estimate the expectation value of a Pauli observable after a stabiliser channel is applied to a non-stabiliser input state. The minimum number of samples required to achieve some stated accuracy scales with $\mathcal{R}^{2}$.

The robustness of magic can be calculated using standard linear programming techniques Boyd and Vandenberghe (2004) (for example using the MATLAB package CVX Grant and Boyd (2014)). The naive formulation of the linear program is practical on a desktop computer for up to five qubits (the number of stabiliser states increases super-exponentially with $n$). It was recently shown by Heinrich and Gross Heinrich and Gross (2019) that when states possess certain symmetries, the original optimisation problem can be mapped to a more tractable one, so that the robustness of magic can be calculated for up to 10 copies of a state.

The framework naturally extends to a subclass of non-stabiliser circuits: those that may be implemented by deterministic state injection Howard and Campbell (2017), including all gates from the third level of the Clifford hierarchy (Figure 2). The canonical example is the T-gate, $T=\mathrm{diag}\quantity(1,e^{i\pi/4})$, which can be implemented by consuming so-called magic states as a resource Bravyi and Kitaev (2005). The classical simulation overhead for implementing a gate is then the robustness of magic for the consumed resource state. Not all non-stabiliser operations can be implemented in this way, however.

Informally we say that an operation is stabiliser-preserving if it always maps stabiliser states to stabiliser states. To make this precise, define $\mathrm{SP}_{n,m}$ to be the set of $n$-qubit operations $\mathcal{E}$ such that $\quantity(\mathcal{E}\otimes\mathds{1}_{m})\sigma\in\mathrm{STAB}_{n+m}$ for all $\sigma\in\mathrm{STAB}_{n+m}$, where $\mathds{1}_{m}$ is the identity map for an $m$-qubit Hilbert space. The set $\mathrm{SP}_{n,0}$ is then the set of channels that map $n$-qubit stabiliser states to $n$-qubit stabiliser states. We say a channel is “completely” stabiliser-preserving if $\mathcal{E}\in\mathrm{SP}_{n,m}$ for all $m$.

3 Overview of main results

Our first result is a characterisation of the class of completely stabiliser-preserving operations, making use of the well-known Choi-Jamiołkowski isomorphism Jamiolkowski (1972); Choi (1975); Jiang et al. (2013). Throughout this paper we will consider Hilbert spaces comprised of $n$-qubit subspaces, eg. $\mathcal{H}=\mathcal{H}^{A}\otimes\mathcal{H}^{B}$. Where it is unclear from the context, we use superscripts to indicate which subspace each operator belongs to. For the maximally entangled state $\ket{\Omega_{n}}$, as defined in Theorem 3.1 below, we use the superscript $A|B$ to indicate the partition across which the state is entangled.

Theorem 3.1** (Completely stabiliser-preserving operations).**

Given an $n$-qubit CPTP channel $\mathcal{E}$, for all $m>0$, $\mathcal{E}\in\mathrm{SP}_{n,n+m}$ if and only if $\mathcal{E}\in\mathrm{SP}_{n,n}$. Furthermore, $\mathcal{E}\in\mathrm{SP}_{n,n}$ if and only if the Choi state

[TABLE]

is a stabiliser state. Here, $\ket{j}$ are the $n$-qubit computational basis states.

We prove this in section 4 . We take this to be the set of free operations in our resource theories.

Our first new monotone is channel robustness $\mathcal{R}_{*}$. For an $n$-qubit CPTP channel $\mathcal{E}$ this is defined as:

[TABLE]

where $\Lambda_{\pm}$ are completely stabiliser-preserving and CPTP maps. To fully enumerate this class of maps, we notice that the associated Choi state must satisfy two conditions: (i) $\Phi_{\mathcal{E}}^{AB}$ is a stabiliser state, and (ii) $\Phi_{\mathcal{E}}^{AB}$ satisfies the trace-preservation condition $\Tr_{A}(\Phi_{\mathcal{E}})=\frac{\mathds{1}_{n}}{2^{n}}$. We can therefore write:

[TABLE]

This can now be calculated by linear program given access to a list of all stabiliser states (see Appendix C and the code repository Ref. Seddon (2019)). Channel robustness satisfies the following:

Faithfulness: If $\mathcal{E}$ is a CPTP channel, then $\mathcal{R}_{*}\quantity(\mathcal{E})=1$ if $\mathcal{E}$ is completely stabiliser-preserving and strictly larger than $1$ otherwise; 2. 2.

Convexity: $\mathcal{R}_{*}\quantity(\sum_{j}q_{j}\mathcal{E}_{j})\leq\sum_{j}\absolutevalue{q_{j}}\mathcal{R}_{*}\quantity(\mathcal{E}_{j})$; 3. 3.

Submultiplicativity under composition: $\mathcal{R}_{*}\quantity(\mathcal{E}_{2}\circ\mathcal{E}_{1})\leq\mathcal{R}_{*}\quantity(\mathcal{E}_{2})\mathcal{R}_{*}\quantity(\mathcal{E}_{1})$; 4. 4.

Submultiplicativity under tensor product: $\mathcal{R}_{*}\quantity(\mathcal{E}^{A}\otimes\mathcal{E}^{\prime B})\leq\mathcal{R}_{*}\quantity(\mathcal{E}^{A})\mathcal{R}_{*}\quantity(\mathcal{E}^{\prime B})$.

As a special case, if $\Lambda$ is a CPTP stabiliser channel, then

[TABLE]

and similarly $\mathcal{R}_{*}(\mathcal{E}\circ\Lambda)\leq\mathcal{R}_{*}(\mathcal{E})$. This combines submultiplicativity under composition and faithfulness, to show that $\mathcal{R}_{*}$ is suitably monotonically non-increasing under compositions with stabiliser channels. This is the sense in which channel robustness is a magic monotone for channels. We prove submultiplicativity in Section 5. For completeness we prove convexity and faithfulness in Appendix B.

The approach above is very close to the stabiliser decomposition of channels employed by Bennink et al. in Ref. Bennink et al. (2017). The main difference is that Bennink et al. optimise their decomposition with respect to $\mathrm{CPR}$, the set of Clifford unitaries supplemented by Pauli reset channels, rather than $\mathrm{SP}_{n,n}$. The set $\mathrm{CPR}$ turns out to be a strict subset of the stabiliser-preserving CPTP maps, so $\mathcal{R}_{*}$ is a lower bound to the $\ell_{1}$-norm of any $\mathrm{CPR}$ decomposition (though the bound is tight in many cases). Just as the $\ell_{1}$-norm in Bennink et al. quantifies the sample complexity of a classical simulation algorithm, we can construct a related algorithm where the runtime depends on $\mathcal{R}_{*}$ in a similar way. We give the details of this algorithm in Section 7.1.

Before proceeding, let us reflect on the condition $\Tr_{A}(\rho_{\pm})=\mathds{1}_{n}/2^{n}$ that enforces that the corresponding channels $\mathcal{E}_{\pm}$ are trace-preserving. Dropping this condition would instead lead to $\mathcal{R}(\Phi_{\mathcal{E}})$, the robustness of the Choi state. For gates from the third level of the Clifford hierarchy, deterministic state injection is always possible, and hence the resource cost of the gate $U$ can be equated with the robustness of magic of the corresponding resource state. These resource states can always (by Clifford-equivalence) be taken to have the form $\outerproduct{U}{U}=\quantity(U_{A}\otimes\mathds{1}_{B})\outerproduct{\Omega}{\Omega}$. This is precisely the Choi state, so it is natural to ask if $\mathcal{R}(\Phi_{\mathcal{E}})$ also quantifies non-stabiliserness for more general channels. We find that $\mathcal{R}(\Phi_{\mathcal{E}})$ exhibits faithfulness, convexity and submultiplicativity under tensor product, but lacks submultiplicativity under composition. This arises from the fact that the decomposition of the Choi state corresponds to a decomposition of the channel into maps that are not necessarily trace-preserving. See Appendix A for a counterexample. Despite this shortcoming, we will see that $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$ is a useful quantity to compare to more well-behaved measures. Moreover, given that $\mathcal{R}(\Phi_{\mathcal{E}})$ does give the resource cost for third level Clifford hierarchy gates, for consistency of the framework it is reasonable to require that our new monotones should be equal to $\mathcal{R}(\Phi_{\mathcal{E}})$ for this restricted class of gates. We will show later that our monotones do have this property.

Our second new monotone is the magic capacity. Given an $n$-qubit channel $\mathcal{E}$, it is natural to consider the largest possible increase in robustness of magic, over any possible input state. By analogy with the resource theories of entanglement Campbell (2010) and coherence Stahlke (2014), we define the magic capacity as:

[TABLE]

Note that the definition of capacity involves forming a tensor product of an $n$-qubit channel $\mathcal{E}$ with the $n$-qubit identity. This is necessary because there exist $n$-qubit channels that generate their maximum robustness when applied to part of an $m$-qubit state, where $m>n$. Nevertheless, the $n$-qubit identity suffices for our definition; This is a consequence of Lemma 4.1 in Section 4. The capacity has the following useful properties:

Faithfulness: If $\mathcal{E}$ is a completely positive, trace-preserving (CPTP) channel, then $\mathcal{C}=1$ if $\mathcal{E}$ is stabiliser-preserving (SP), and strictly larger than $1$ otherwise; 2. 2.

Convexity: $\mathcal{C}\quantity(\sum_{j}q_{j}\mathcal{E}_{j})\leq\sum_{j}|q_{j}|\mathcal{C}\quantity(\mathcal{E}_{j})$; 3. 3.

Submultiplicativity under composition: $\mathcal{C}\quantity(\mathcal{E}_{1}\circ\mathcal{E}_{2})\leq\mathcal{C}\quantity(\mathcal{E}_{1})\mathcal{C}\quantity(\mathcal{E}_{2})$; 4. 4.

Submultiplicativity under tensor product: $\mathcal{C}\quantity(\mathcal{E}^{A}\otimes\mathcal{E}^{\prime B})\leq\mathcal{C}\quantity(\mathcal{E}^{A})\mathcal{C}\quantity(\mathcal{E}^{\prime B})$; 5. 5.

Maximum increase in robustness: $\frac{\mathcal{R}[\quantity(\mathcal{E}\otimes\mathds{1})\rho]}{\mathcal{R}\quantity(\rho)}\leq\mathcal{C}\quantity(\mathcal{E}),\,\forall\rho$.

In Section 6 we prove properties 3-5, (convexity and faithfulness are shown in Appendix D). We will also prove the following theorem relating magic capacity to channel robustness and $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$.

Theorem 3.2** (Sandwich Theorem).**

For any CPTP map $\mathcal{E}$, the following inequalities hold:

[TABLE]

Moreover, if the unitary operation $\mathcal{U}$ is in the third level of the Clifford hierarchy, then we have equality:

[TABLE]

We are interested in whether or not these inequalities are tight for more general operations. In Table 1 we summarise numerical results for a selection of diagonal gates. The results for these gates are presented in full in Section 8.

The magic capacity also quantifies the sample complexity for a Monte Carlo-type classical simulation algorithm, presented in Section 7.2. This differs from previous algorithms such as Bennink et al. Bennink et al. (2017) in that a convex optimisation must be solved at each step. While this results in an increase in runtime per sample, it can be the case that $\mathcal{C}(\mathcal{E})\ll\mathcal{R}_{*}(\mathcal{E})$, which can lead to an improvement in sample complexity over the algorithm of Section 7.1.

4 Completely stabiliser-preserving operations

In this section, we justify setting $\mathrm{SP}_{n,n}$ as the class of free operations. We begin with an example channel $\mathcal{E}\in\mathrm{SP}_{n,0}$ that fails to be stabiliser-preserving when acting on part of a larger system. Consider the single-qubit channel $\mathcal{E}_{T}$ defined by the Kraus operators $\quantity{\outerproduct{0}{T},\outerproduct{1}{T_{\perp}}}$, where $\ket{T}=T\ket{+}$ and $\ket{T_{\perp}}=T\ket{-}$. Clearly, applied to any single-qubit state, the output will be some probabilistic mixture of $\ket{0}$ and $\ket{1}$, and so must have $\mathcal{R}=1$, so $\mathcal{E}_{T}\in\mathrm{SP}_{n,0}$. But if $\mathcal{E}_{T}$ is applied to one qubit in a Bell pair, we obtain:

[TABLE]

where $\ket{T^{*}}=T^{\dagger}\ket{+}$, $\ket{T^{*}_{\perp}}=T^{\dagger}\ket{-}$. From this output state, we can deterministically recover a pure magic state on qubit 2 using only stabiliser operations, by making a $Z$-measurement on qubit 1 and then performing a rotation on qubit 2 conditioned on the outcome. The output state has robustness $\mathcal{R}\quantity(\quantity(\mathcal{E}_{T}\otimes\mathds{1})\outerproduct{\Phi_{+}}{\Phi_{+}})=\mathcal{R}\quantity(\ket{T})=\sqrt{2}$.

So, there exist channels where $\mathcal{E}\in\mathrm{SP}_{n,m}$ but $\mathcal{E}\notin\mathrm{SP}_{n,m+1}$. To call a channel completely stabiliser-preserving, then, we need to be sure $\mathcal{E}\otimes\mathds{1}_{m}$ remains stabiliser-preserving for all $m>0$. We now show we only need tensor with the identity of the same dimension as the original channel.

Lemma 4.1** (Maximum robustness achieved on $2n$ qubits).**

Let $\mathcal{E}$ be an $n$-qubit quantum channel. Then for $m>0$, for any $\ket{\phi}\in\mathrm{STAB}_{2n+m}$, there exists some state $\ket{\psi}\in\mathrm{STAB}_{2n}$ such that:

[TABLE]

Proof.

Consider a $(2n+m)$-qubit stabiliser state $\ket{\phi}$, with partition $A|B$ between the first $n$ and last $n+m$ qubits. Ref. Fattal et al. (2004) shows that the state $\ket{\phi}^{AB}$ is local Clifford-equivalent to $p$ independent Bell pairs entangled across the partition $A|B$ (here “local” means with respect to the bipartition rather than per qubit). Since there are $n$ qubits in partition $A$, $p$ is at most $n$. Let $B^{\prime}|B^{\prime\prime}$ be a partition of $B$ into $n$ and $m$ qubits. Then by local permutation of qubits within B, we can take these $p\leq n$ Bell pairs to be entangled across $A|B^{\prime}$. So we have:

[TABLE]

where $U^{B}$ is a Clifford operation, $\ket{\psi}^{AB^{\prime}}\in\mathrm{STAB}_{2n}$ and $\ket{\psi^{\prime}}^{B^{\prime\prime}}\in\mathrm{STAB}_{m}$. So writing the channel corresponding to $U^{B}$ as $\mathcal{U}^{B}$, for any $\mathcal{E}$ on $n$ qubits, we know that:

[TABLE]

Since $\mathds{1}_{n}\otimes\mathcal{U}^{B}$ represents a (reversible) Clifford gate, by monotonicity of robustness of magic:

[TABLE]

where in the last line we used the fact that $\ket{\psi^{\prime}}^{B^{\prime\prime}}$ is a stabiliser state, and hence does not contribute to the robustness. The state $\outerproduct{\psi}{\psi}^{AB^{\prime}}$ is a $2n$-qubit state, so this proves the result. ∎

This lemma allows us to prove the first claim of Theorem 3.1, which says that $\mathcal{E}$ is completely stabiliser-preserving if and only if $\mathcal{E}\in\mathrm{SP}_{n,n}$. The inclusion $\mathrm{SP}_{n,n}\subseteq\mathrm{SP}_{n,n+m}$ is immediate since the stabiliser states are preserved under tracing out of auxiliary systems. The interesting inclusion is $\mathrm{SP}_{n,n+m}\subseteq\mathrm{SP}_{n,n}$. Suppose that $\mathcal{E}\in\mathrm{SP}_{n,n}$ and consider any $\sigma\in\mathrm{STAB}_{2n+m}$. By Lemma 4.1 there exists some stabiliser state $\sigma^{\prime}\in\mathrm{STAB}_{2n}$ such that $\mathcal{R}\quantity(\quantity(\mathcal{E}\otimes\mathds{1}_{n+m})\sigma)=\mathcal{R}\quantity(\quantity(\mathcal{E}\otimes\mathds{1}_{n})\sigma^{\prime})$. But if $\mathcal{E}\in\mathrm{SP}_{n,n}$, then $\quantity(\mathcal{E}\otimes\mathds{1}_{n})\sigma^{\prime}$ is a stabiliser state, so the robustness is equal to 1. By the faithfulness of robustness of magic, $\quantity(\mathcal{E}\otimes\mathds{1}_{n+m})\sigma$ is a stabiliser state. Therefore, $\mathcal{E}\in\mathrm{SP}_{n,n}$ implies $\mathcal{E}\in\mathrm{SP}_{n,n+m}$. Next, we discuss a straightforward test for membership of this set, which does not require mechanically checking all possible input stabiliser states.

We can associate every $n$-qubit channel $\mathcal{E}$ with a unique density operator on $2n$ qubits Jiang et al. (2013); Jamiolkowski (1972); Choi (1975):

[TABLE]

Here $\ket{j}$ label the computational basis states. We will also use the following property:

[TABLE]

Consider the robustness of magic of the Choi state, $\mathcal{R}(\Phi_{\mathcal{E}})$. We mentioned earlier that $\mathcal{R}(\Phi_{\mathcal{E}})$ quantifies simulation cost for gates from the third level of the Clifford hierarchy. This motivates us to consider its properties for more general operations, and it turns out that $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$ gives us our first criterion for completely stabiliser-preserving channels.

Lemma 4.2** (Faithfulness of robustness of the Choi state).**

Consider the $n$-qubit CPTP channel $\mathcal{E}$. If $\mathcal{E}\in\mathrm{SP}_{n,n}$, then $\mathcal{R}\quantity(\Phi_{\mathcal{E}})=1$. Otherwise, $\mathcal{R}\quantity(\Phi_{\mathcal{E}})>1$.

Proof.

The fact that $\mathcal{E}\in\mathrm{SP}_{n,n}$ implies $\mathcal{R}\quantity(\Phi_{\mathcal{E}})=1$ is easy to see. Since $\outerproduct{\Omega_{n}}{\Omega_{n}}$ is itself a $2n$-qubit stabiliser state, $\mathcal{E}\in\mathrm{SP}_{n,n}$ guarantees that $\Phi_{\mathcal{E}}$ is a stabiliser state, so must have robustness 1. The implication in the other direction is less obvious; one might imagine there perhaps exist maps that send $\outerproduct{\Omega_{n}}{\Omega_{n}}$ in particular to a stabiliser state, but are not stabiliser-preserving in general. We show that this is not the case using an argument based on witnesses for non-stabiliser states, in part inspired by the conditions for free operations (SPO) given by Ahmadi et al. Ahmadi et al. (2018) for odd prime dimension qudits and for the single qubit case. The criteria for SPO were based on a class of witness defined by phase point operators. Here we instead consider the following family of witnesses for $n$ qubits. We say that $W_{n}$ is a good witness for $n$-qubit non-stabiliser states if:

[TABLE]

The hyperplane separation theorem Boyd and Vandenberghe (2004) guarantees that such witnesses exist and can be constructed for any non-stabiliser state $\rho$. That is, for any $\rho\notin\mathrm{STAB}_{n}$, there always exists an operator $W_{\rho}$ such that $\Tr(W_{\rho}\rho)>0$ and yet is a good $n$-qubit witness as defined above.

We first show that for any good $n$-qubit witness $W_{n}$, the operator $W_{n}\otimes\outerproduct{\phi}{\phi}$, where $\ket{\phi}\in\mathrm{STAB}_{m}$, is a good witness for $(n+m)$-qubit non-stabiliser states. For any $\sigma\in\mathrm{STAB}_{n+m}$:

[TABLE]

so that ${\widetilde{\sigma}=\quantity(\mathds{1}_{n}\otimes\outerproduct{\phi}{\phi})\sigma\quantity(\mathds{1}_{n}\otimes\outerproduct{\phi}{\phi})}$, where we used cyclicity of the trace and the fact that $\outerproduct{\phi}{\phi}$ is a projector. If $\widetilde{\sigma}=0$ then the inequality (17) is trivially always satisfied by $W_{n}\otimes\outerproduct{\phi}{\phi}$. Otherwise, $\widetilde{\sigma}$ is a stabiliser state (non-normalised) and so too is $\Tr_{B}[\widetilde{\sigma}]$. Then:

[TABLE]

The inequality follows because $W_{n}$ is a good witness and $\Tr_{B}[\widetilde{\sigma}]$ is a stabiliser state. Therefore, $W_{n}\otimes\outerproduct{\phi}{\phi}$ is also a valid witness.

Now suppose $\mathcal{E}\notin\mathrm{SP}_{n,n}$. Then there is some stabiliser state $\ket{\phi^{\prime}}\in\mathrm{STAB}_{2n}$, such that $\rho^{\prime}=\quantity(\mathcal{E}\otimes\mathds{1}_{n})\outerproduct{\phi^{\prime}}{\phi^{\prime}}\notin\mathrm{STAB}_{2n}$. By the hyperplane separation theorem, there exists a good $2n$-qubit witness $W_{\rho^{\prime}}$ such that $\Tr[W_{\rho^{\prime}}\rho^{\prime}]>0$. Consider that the $4n$-qubit state $\ket{\Omega_{2n}}^{AA^{\prime}|BB^{\prime}}$ is unentangled across the partition $AB|A^{\prime}B^{\prime}$, so that we can write:

[TABLE]

taking care to note the permutation of subspaces. Therefore the Choi state for $\quantity(\mathcal{E}\otimes\mathds{1}_{n})$ is:

[TABLE]

We then use equation (16) to obtain:

[TABLE]

But $(W_{\rho^{\prime}}\otimes\outerproduct{\phi^{\prime}}{\phi^{\prime}}^{T})$ is a good witness, so $\Phi_{\mathcal{E}}\otimes\outerproduct{\Omega_{n}}{\Omega_{n}}\notin\mathrm{STAB}_{4n}$ and therefore $\Phi_{\mathcal{E}}$ is a non-stabiliser state. So, by faithfulness of robustness of magic, if $\mathcal{E}\notin\mathrm{SP}_{n,n}$ then $\mathcal{R}\quantity(\Phi_{\mathcal{E}})>1$. ∎

Combined the above two lemmas provide a proof of both claims given in Theorem 3.1. This does not mean the robustness of the Choi state is a reliable monotone, since despite being faithful it fails to be submultiplicative under composition (see Appendix A). Rather, we use the faithfulness of the Choi state as a tool to give an alternative definition of the channel robustness as captured by Eq. (5).

5 Channel robustness

A natural extension of the algorithm of Bennink et al. Bennink et al. (2017) is to replace $\mathrm{CPR}$ (the set of Clifford gates and Pauli reset channels) with $\mathrm{SP}_{n,n}$. We therefore define the channel robustness as:

[TABLE]

To calculate this in practice, we decompose the Choi state $\Phi_{\mathcal{E}}$ as per equation (5), adapting the robustness of magic optimisation problem from Ref. Howard and Campbell (2017). The details are given in Appendix C. We also note that for diagonal channels $\mathcal{E}$, the problem is equivalent to a decomposition of the state $\mathcal{E}(\outerproduct{+}{+}^{\otimes n})$ as:

[TABLE]

Here, the condition on the partial trace for the general case is replaced by the requirement that all diagonal elements of the states $\rho_{\pm}$ are equal to $1/2^{n}$. In the special case where $\rho_{\pm}$ are pure, this implies they are diagonal Clifford-equivalent to graph states, though more generally they may be mixtures of states that individually do not have full support in the standard basis. In practical terms, this reduction to an $n$-qubit problem confers a significant advantage, since the number of stabiliser states (which form the extreme points of the linear programming problem) grows super-exponentially with $n$. Full technical details of this simplification are given in Appendix E. We now return to consider the properties of channel robustness.

The channel robustness is convex and faithful with these properties inherited from the robustness of magic (see Appendix B for details). Here we discuss additional properties.

Submultiplicativity under composition: $\mathcal{R}_{*}(\mathcal{E}_{2}\circ\mathcal{E}_{1})\leq\mathcal{R}_{*}(\mathcal{E}_{1})\mathcal{R}_{*}(\mathcal{E}_{2})$. The channels $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ will have an optimal decomposition:

[TABLE]

where $\mathcal{R}_{*}(\mathcal{E}_{j})=1+2p_{j}$ and $\Lambda_{j,\pm}$ are CPTP maps and completely stabiliser preserving. Using these decompositions, we obtain that

[TABLE]

where

[TABLE]

The set of CPTP completely stabiliser preserving channels is closed under composition and convex, so both $\Lambda_{\pm}^{\prime}$ are in this set. Therefore, we have a valid decomposition for $\mathcal{E}_{2}\circ\mathcal{E}_{1}$ that entails $\mathcal{R}_{*}(\mathcal{E}_{2}\circ\mathcal{E}_{1})\leq 1+2q$. One finds

[TABLE]

which completes the proof.

Submultiplicativity under tensor product: $\mathcal{R}_{*}\quantity(\mathcal{E}^{A}\otimes\mathcal{E}^{\prime B})\leq\mathcal{R}_{*}\quantity(\mathcal{E}^{A})\mathcal{R}_{*}\quantity(\mathcal{E}^{\prime B})$. We treat tensor product as a special case of composition. For $n$-qubit $\mathcal{E}_{A}$ and $m$-qubit $\mathcal{E}^{\prime}_{B}$:

[TABLE]

To complete the proof we will confirm that

[TABLE]

As noted earlier, we can write $\ket{\Omega_{n+m}}^{AA^{\prime}|BB^{\prime}}=\ket{\Omega_{n}}^{A|B}\otimes\ket{\Omega_{m}}^{A^{\prime}|B^{\prime}}$, so that the Choi state for $\mathcal{E}^{A}\otimes\mathds{1}_{m}$ is given by:

[TABLE]

The state $\Phi_{\mathcal{E}}$ will have some optimal decomposition $\Phi_{\mathcal{E}}={(1+p)\rho_{+}-p\rho_{-}}$, with channel robustness $\mathcal{R}_{*}(\mathcal{E}_{A})=1+2p$, so that:

[TABLE]

This is a valid, not necessarily optimal, stabiliser decomposition satisfying the trace condition, so we have

[TABLE]

This is enough to show submultiplicativity; for completeness, in Appendix B we will also show ${\mathcal{R}_{*}(\mathcal{E})\leq\mathcal{R}_{*}(\mathcal{E}\otimes\mathds{1})}$ so that in fact we have equality.

6 Magic capacity

6.1 Properties

We now turn to our second monotone, which quantifies the capacity of a channel to generate magic. Recall:

[TABLE]

where $\mathcal{R}$ is the robustness of magic. Notice that we only need optimise over the pure stabiliser states. For mixed states or even non-stabiliser states, the capacity still captures the possible increase in robustness of magic by virtue of the maximum increase in robustness property:

[TABLE]

Here we prove this property, using similar arguments to those deployed in Campbell (2010). Consider an $n$-qubit channel $\mathcal{E}$. Any $2n$-qubit input state $\rho$ will have an optimal stabiliser state decomposition $\rho=\sum_{j}q_{j}\outerproduct{\phi_{j}}{\phi_{j}}$, where $\sum_{j}q_{j}=1$, and such that $\mathcal{R}\quantity(\rho)=\sum_{j}|q_{j}|$. By linearity we have:

[TABLE]

By convexity of robustness of magic, we then have:

[TABLE]

The optimal pure stabiliser state $\ket{\phi_{*}}$, satisfies:

[TABLE]

for any $j$. So we have:

[TABLE]

Rearranging we obtain inequality (37).

Submultiplicativity under composition: $\mathcal{C}\quantity(\mathcal{E}_{1}\circ\mathcal{E}_{2})\leq\mathcal{C}\quantity(\mathcal{E}_{1})\mathcal{C}\quantity(\mathcal{E}_{2})$. Take the composition of two linear maps $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$. There exists some stabiliser state $\rho_{*}=\outerproduct{\phi_{*}}{\phi_{*}}$ that achieves the optimal robustness:

[TABLE]

The operator $(\mathcal{E}_{1}\otimes\mathds{1}_{n})\rho_{*}$ will have some optimal decomposition $(\mathcal{E}_{1}\otimes\mathds{1}_{n})\rho_{*}=\sum_{k}q_{1k}\outerproduct{\phi_{k}}{\phi_{k}}$ such that $\mathcal{R}\quantity((\mathcal{E}_{1}\otimes\mathds{1}_{n})\rho_{*})=\sum_{k}\absolutevalue{q_{1k}}$. So by linearity:

[TABLE]

Then by convexity of robustness of magic:

[TABLE]

where to go from (44) to (45) we used the fact that since $\ket{\phi_{k}}$ are stabiliser states, $\mathcal{R}\quantity((\mathcal{E}_{2}\otimes\mathds{1})\outerproduct{\phi_{k}}{\phi_{k}})$ can be no larger than $\mathcal{C}\quantity(\mathcal{E}_{2})$. Finally, using the fact that $\mathcal{R}\quantity((\mathcal{E}_{1}\otimes\mathds{1})\rho_{*})\leq\mathcal{C}\quantity(\mathcal{E}_{1})$, we have $\mathcal{C}\quantity(\mathcal{E}_{2}\circ\mathcal{E}_{1})\leq\mathcal{C}\quantity(\mathcal{E}_{2})\mathcal{C}\quantity(\mathcal{E}_{1})$, completing the proof.

Submultiplicativity under tensor product: $\mathcal{C}\quantity(\mathcal{E}^{A}\otimes\mathcal{E}^{\prime B})\leq\mathcal{C}\quantity(\mathcal{E}^{A})\mathcal{C}\quantity(\mathcal{E}^{\prime B})$. This follows directly from submultiplicativity under composition, since

[TABLE]

We saw in Section 4 that any gains in robustness achievable by tensoring $\mathcal{E}_{A}$ with the identity and acting on a larger Hilbert space are already taken care of by the $\otimes\mathds{1}_{n}$ in the definition (36), so that $\mathcal{C}\quantity(\mathcal{E}^{A}\otimes\mathds{1}_{m}^{B})=\mathcal{C}\quantity(\mathcal{E}^{A})$ and $\mathcal{C}\quantity(\mathds{1}_{n}^{A}\otimes\mathcal{E}^{\prime B})=\mathcal{C}\quantity(\mathcal{E}^{\prime B})$. Substituting this into inequality (47) gives the desired result.

6.2 Sandwich theorem

We will now prove Theorem 3.2, which stated that $\mathcal{R}\quantity(\Phi_{\mathcal{E}})\leq\mathcal{C}\quantity(\mathcal{E})\leq\mathcal{R}_{*}\quantity(\mathcal{E})$, for any CPTP channel $\mathcal{E}$, and that $\mathcal{R}\quantity(\Phi_{\mathcal{U}})=\mathcal{C}\quantity(\mathcal{U})=\mathcal{R}_{*}\quantity(\mathcal{U})$ for any unitary operation $\mathcal{U}$ from the third level of the Clifford hierarchy.

Proof.

By definition $\Phi_{\mathcal{E}}=\quantity(\mathcal{E}\otimes\mathds{1}_{n})\outerproduct{\Omega_{n}}{\Omega_{n}}$. But $\ket{\Omega_{n}}$ is a stabiliser state, so $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$ can be no larger than $\mathcal{R}\quantity(\quantity(\mathcal{E}\otimes\mathds{1}_{n})\outerproduct{\phi_{*}}{\phi_{*}})=\mathcal{C}(\mathcal{E})$, where $\ket{\phi_{*}}$ is the stabiliser state that achieves the capacity, and so:

[TABLE]

Now suppose $\mathcal{E}=(1+p)\Lambda_{+}-p\Lambda_{-}$ is the optimal decomposition of $\mathcal{E}$ into CPTP stabiliser-preserving maps, $\Lambda_{\pm}\in\mathrm{SP}_{n,n}$, so that $\mathcal{R}_{*}\quantity(\mathcal{E})=1+2p$. Then for any input stabiliser state $\sigma\in\mathrm{STAB}_{2n}$, we can write down a valid stabiliser decomposition of the output state:

[TABLE]

In particular this is true for the stabiliser state $\sigma_{*}=\outerproduct{\phi_{*}}{\phi_{*}}$ that is optimal with respect to the capacity. But equation (49) could be a non-optimal decomposition, so its $\ell_{1}$-norm $1+2p$ is at least as large as $\mathcal{R}\quantity(\quantity(\mathcal{E}\otimes\mathds{1}_{n})\sigma_{*})$. So:

[TABLE]

completing the proof of the first statement. Having done so, to prove the second statement it suffices to show that $\mathcal{R}\quantity(\Phi_{\mathcal{U}})=\mathcal{R}_{*}\quantity(\mathcal{U})$.

For any $n$-qubit gate $U$ from the third level of the Clifford hierarchy, corresponding to the channel $\mathcal{U}$, deterministic state injection is possible Bravyi and Kitaev (2005); Gottesman and Chuang (1999). That is, given a Hilbert space $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}\otimes\mathcal{H}_{C}\otimes\mathcal{H}_{D}$, where each subspace is comprised of $n$ qubits, there exists a completely-stabiliser-preserving circuit $\Lambda$ such that, for any $2n$-qubit input state $\rho$:

[TABLE]

Where $\Phi_{\mathcal{U}}$ is the Choi state for the channel $\mathcal{U}$. The circuit $\Lambda$ is comprised of a complete Bell measurement on $BC$, followed by a Clifford correction on subspace $A$ conditioned on the outcome of the Bell measurement. It can be represented by Kraus operators:

[TABLE]

where $M_{j}=\mathds{1}^{A}\otimes\outerproduct{\Phi_{j}}{\Phi_{j}}^{BC}\otimes\mathds{1}^{D}$ are the Kraus operators corresponding to elements of the Bell basis $\ket{\Phi_{j}}$, and $C_{j}$ is some unitary Clifford correction.

Now consider an optimal decomposition of the Choi state:

[TABLE]

We now show that by substitution into equation (51) we can obtain a decomposition of the channel that satisfies the trace-preservation condition required for channel robustness. We have:

[TABLE]

where $\widetilde{\rho}_{\pm}^{AD}=\Tr_{BC}[\Lambda(\rho_{\pm}^{AB}\otimes\outerproduct{\Omega}{\Omega}^{CD})]\in\mathrm{STAB}_{2n}$, since $\Lambda\in\mathrm{SP}_{n,n}$. If we can show that $\Tr_{A}(\widetilde{\rho}_{\pm})=\mathds{1}_{n}/2^{n}$, then we have satisfied the required condition. First, note that $\Tr_{A}(\widetilde{\rho}_{\pm})$ is independent of the Clifford corrections $C_{j}$, since the partial trace depends only on the outcome probabilities of the Bell measurement, $p_{j}=\Tr[M_{j}(\rho_{\pm}^{AB}\otimes\outerproduct{\Omega}{\Omega}^{CD})M_{j}^{\dagger}]$. Therefore $\Tr_{A}(\widetilde{\rho}_{\pm})=\Tr_{ABC}(\rho^{\prime}_{\pm})$, where $\rho^{\prime}_{\pm}=\sum_{j}M_{j}(\rho_{\pm}^{AB}\otimes\outerproduct{\Omega}{\Omega}^{CD})M_{j}^{\dagger}$ is the state following the Bell measurement. We then have:

[TABLE]

In going to the second line, we used the fact that the partial trace over $BC$ is cyclic with respect to operators that act non-trivially only on $\mathcal{H}_{B}\otimes\mathcal{H}_{C}$. In going from the second to the third line, we used the fact that $\{M_{j}\}$ is a complete set of Kraus operators, so $\sum_{j}M_{j}^{\dagger}M_{j}=\mathds{1}$. We have shown that the decomposition (55) satisfies the trace preservation criterion. Since the decomposition may not be optimal, we have that $\mathcal{R}_{*}\quantity(\mathcal{U})\leq 1+2p=\mathcal{R}\quantity(\Phi_{\mathcal{U}})$. But from the proof of the first statement $\mathcal{R}\quantity(\Phi_{U})\leq\mathcal{C}\quantity(\mathcal{U})\leq\mathcal{R}_{*}\quantity(\mathcal{U})$, so it must be that equality holds. ∎

We note that the result that $\mathcal{R}(\Phi_{\mathcal{U}})=\mathcal{R}_{*}\quantity(\mathcal{U})$ for third-level gates carries over to the case of decompositions of $\mathcal{U}(\outerproduct{+}{+}^{\otimes n})$ for diagonal third-level gates. That is, there always exists a decomposition satisfying the constraints of equation (24) that is optimal with respect to $\mathcal{R}(\Phi_{\mathcal{U}})=\mathcal{R}(\mathcal{U}(\outerproduct{+}{+}^{\otimes n}))$. This can be seen by following the argument of Theorem 3.2, but replacing the full $4n$-qubit teleportation circuit with a $2n$-qubit state injection circuit (Figure 2).

7 Classical simulation algorithms

Here we propose two classical simulation algorithms. The channel robustness $\mathcal{R}_{*}$ relates to the runtime of our first simulator, which we call the static simulator. The magic capacity $\mathcal{C}$ relates to the runtime of the second simulator, called the dynamic simulator. In both cases, we consider a circuit composed from a sequence of channels with $\{\mathcal{E}_{1},\mathcal{E}_{2},\ldots,\mathcal{E}_{L}\}$ acting on an initial stabiliser state, which we take to be $\ket{0^{n}}$. The circuit ends with some final state $\rho=\mathcal{E}_{L}\ldots\circ\mathcal{E}_{2}\circ\mathcal{E}_{1}(\ket{0^{n}}\bra{0^{n}})$ and measurement of some Pauli observable $Z$. We assume that each channel acts non-trivially on a bounded number of qubits (e.g. 2 or 3) so we can evaluate the relevant monotones. Our goal is to estimate the expectation value $\mathrm{Tr}[Z\rho]$ to within additive error. In the language of Ref. Pashayan et al. (2017) our simulators will be poly-boxes.

Both our algorithms are inspired by previous methods that collect a large number of Monte Carlo samples that scales quadratically with the negativity of some quasiprobability distribution Pashayan et al. (2015); Bennink et al. (2017); Howard and Campbell (2017). The static Monte Carlo simulator uses a precomputed, and therefore static, quasiprobability distribution. The dynamic Monte Carlo simulator recomputes optimal quasiprobability distributions at each step, which can lead to fewer samples required but with a higher runtime per sample. As such, there are subtle trade-offs in the runtime complexities.

Both our algorithms use that completely stabiliser preserving operations $\mathrm{SP}_{n,n}$ acting on a stabiliser state can be classically efficiently simulated. This follows from the fact that given the Choi state $\Phi_{\mathcal{E}}$ for an $n$-qubit channel, the channel may be implemented by performing a Bell measurement on $\Phi_{\mathcal{E}}\otimes\sigma$, postselecting on the $\Omega$ outcome to obtain $(\mathds{1}_{n}\otimes\outerproduct{\Omega}{\Omega})\Phi_{\mathcal{E}}\otimes\sigma(\mathds{1}_{n}\otimes\outerproduct{\Omega}{\Omega})$ and then tracing out the last $2n$ qubits Gottesman and Chuang (1999). This can be simulated using Gottesman-Knill when $\sigma$ is a stabiliser state and $\mathcal{E}\in\mathrm{SP}_{n,n}$. Curiously, it is unclear whether $\mathrm{SP}_{n,n}$ can be physically realised using Clifford unitaries and Pauli measurements but without the use of postselection.

7.1 Static Monte Carlo

In Algorithm 1 we give pseudocode for our first simulator, which we call static Monte Carlo, and which can be viewed as a generalisation of the algorithm of Bennink et al. Bennink et al. (2017), differing in two important ways. First, whereas their algorithm involved an optimisation over the set CPR, the channel robustness is optimised with respect to $\mathrm{SP}_{m,m}$, so that all completely stabiliser-preserving maps are represented non-negatively. Second, while the positive and negative parts of the decomposition are each CPTP, the subroutine in step (c) of the algorithm involves updating the state with maps that may not be trace-preserving. This is necessary to ensure that the updated state after step (c)iii. remains pure. We explain below how this subroutine is carried out and show that it does not increase the sample complexity.

We assume a pre-computation stage in which for each circuit element we determine an optimal decomposition as per the definition (5), so that:

[TABLE]

where $\mathcal{E}_{j,k}\in\mathrm{SP}_{m,m}$. The runtime of this pre-computation stage is bounded since each circuit element acts non-trivially on only $m$ qubits. Once these decompositions have been calculated, one can define an overall quasiprobability distribution:

[TABLE]

where $\vec{k}\in\mathbb{Z}_{2}^{L}$ is a vector representing a choice of either $\mathcal{E}_{j,0}$ or $\mathcal{E}_{j,1}$ at each circuit element $\mathcal{E}_{j}$. We can renormalise this to obtain a product probability distribution:

[TABLE]

Hence we can write the output of the circuit as follows:

[TABLE]

where each $\Lambda_{\vec{k}}=\mathcal{E}_{L,k_{L}}\circ\ldots\circ\mathcal{E}_{1,k_{1}}$ gives a trajectory of $\mathrm{SP}_{m,m}$ channels through the circuit, $\lambda_{\vec{k}}=\mathrm{sign}(q_{\vec{k}})$, and $R=\prod_{j=1}^{L}\mathcal{R}_{*}\quantity(\mathcal{E}_{j})$. The expectation value of an observable $Z$ at the end of the circuit is given by:

[TABLE]

This decomposition is a quasi-probability distribution with $\ell_{1}$-norm

[TABLE]

By sampling from the stabiliser-preserving trajectories $\Lambda_{\vec{k}}$ with probability distribution $\{p_{\vec{k}}\}$, and calculating $\lambda_{\vec{k}}R\Tr[Z\Lambda_{\vec{k}}(\outerproduct{\phi_{0}}{\phi_{0}})]$, we have an unbiased estimator for $\Tr[Z\mathcal{E}(\outerproduct{\phi_{0}}{\phi_{0}})]$. Crucially, (62) gives a product distribution, so this first sampling step is efficient. The variance is increased by $R>1$, but by standard arguments Pashayan et al. (2015); Bennink et al. (2017); Howard and Campbell (2017) the Hoeffding inequalities can be used to show that any constant error $\delta>0$ in the mean estimate can be achieved with arbitrary success probability $(1-\epsilon)$ by repeating the sampling procedure $N$ times, where:

[TABLE]

The simulator proceeds by tracking the evolution of a pure stabiliser state through the sampled trajectory $\Lambda_{\vec{k}}$, as described in Algorithm 1. In practice, while each $\mathcal{E}_{j,k_{j}}$ is constrained to be CPTP, the output of the linear program will be a decomposition of $\mathcal{E}_{j,k_{j}}$ into maps that are stabiliser-preserving but not necessarily trace-preserving, corresponding to Choi states that are proportional to pure stabiliser projectors:

[TABLE]

Nevertheless, because $\mathcal{E}_{j,k_{j}}$ is a CPTP channel, it can be simulated by sampling from a proper probability distribution as defined in step (c)i. of Algorithm 1. We stress that while $\Phi_{\mathcal{E}_{j,k_{j}}}$ is subject to the constraint ensuring that each $\mathcal{E}_{j,k_{j}}$ is trace-preserving, this need not apply to the individual terms $\outerproduct{\psi_{l,j,k_{j}}}{\psi_{l,j,k_{j}}}$ in its decomposition. Indeed, any decomposition of a non-unital channel cannot be represented as a linear combination of unitary Clifford operations alone. While Clifford gates are represented in the distribution by maximally entangled states, product states correspond to non-trace-preserving maps involving projections. For example, $\Phi_{\mathcal{T}}=\outerproduct{01}{01}$ corresponds to the single Kraus operator $\mathcal{T}=\outerproduct{0}{1}$ which can be seen as a $Z$-measurement post-selected on the “$\ket{1}$” outcome, followed by an $X$ gate. Under these well-defined stabiliser operations, the state update in step (c)iii. can be carried out efficiently as per the Gottesman-Knill theorem Gottesman (1997); Aaronson and Gottesman (2004). More generally, the state update corresponding to any pure Choi state $\ket{\psi_{l,j,k_{j}}}$ can be modelled as a post-selected Bell measurement Gottesman and Chuang (1999):

[TABLE]

Note that the new state $\ket{\phi_{j}}$ is guaranteed to be pure, since the projection of the system BC onto a Bell state removes any correlation across the partition $A|BC$. We emphasise that the probability for carrying out this update depends not just on the known prefactors $|q_{l,j,k_{j}}|$, but on $\Tr[\mathcal{T}_{l,j,k_{j}}(\outerproduct{\phi_{j-1}}{\phi_{j-1}})]$, which are input state dependent. Computing this trace amounts to evaluating the overlap between $(\mathds{1}_{n}\otimes\outerproduct{\Omega_{n}}{\Omega_{n}})$ and some stabiliser state, as per equation (68), which can be done efficiently using the stabiliser tabulex method Gottesman (1997); Aaronson and Gottesman (2004). However, since $\ket{\phi_{j-1}}$ will have been chosen randomly in the previous step of the algorithm, these traces cannot be calculated ahead of time. This means that our algorithm involves additional per-sample runtime overhead in order to calculate $\{p_{l}^{(j)}\}$. However, there are a finite and bounded number of such calculations since each circuit element acts on a small number of qubits. The key point is that the constraint on $\mathcal{E}_{j,k_{j}}$ ensures that, despite being composed of non-trace-preserving elements, the distribution over the outcomes $\{p_{l}^{(j)}\}$ for any given state forms a proper probability distribution. Consequently, the intermediate sampling in step (c)i. does not increase the variance, so the sample complexity of the simulator depends only on $\prod_{j}\mathcal{R}_{*}(\mathcal{E}_{j})^{2}$.

Suppose we want to compare the runtime of our simulator for a particular circuit with that of the Bennink et al. algorithm Bennink et al. (2017). Their decompositions are in terms of $\mathrm{CPR}$, the set of Cliffords and Pauli reset channels, and an associated cost function is

[TABLE]

The sample complexity for simulating a given circuit element $\mathcal{E}_{j}$ is proportional to $\mathcal{R}_{\mathrm{CPR}}(\mathcal{E}_{j})^{2}$. Since $\mathrm{CPR}\subseteq\mathrm{SP}_{n,n}$, it must be the case that $\mathcal{R}_{*}\leq\mathcal{R}_{\mathrm{CPR}}$, potentially leading to lower simulation sample complexity if there exist channels with $\mathcal{R}_{*}<\mathcal{R}_{\mathrm{CPR}}$. We give here a simple toy example demonstrating a significant advantage to our static simulator.

Consider the single-qubit CPTP map $\Lambda_{H}$ defined by a $Z$-measurement followed by a Hadamard gate conditioned on the “-1” outcome. This has Kraus representation:

[TABLE]

This is clearly a completely stabiliser-preserving map, so has channel robustness $\mathcal{R}_{*}\quantity(\Lambda_{H})=1$. For a single qubit, $\mathrm{CPR}$ consists of the 24 Clifford gates, and 6 Pauli reset channels. Using this set, we calculate $\mathcal{R}_{\mathrm{CPR}}\quantity(\Lambda_{H})=2$. Since $\Lambda_{H}\in\mathrm{SP}_{1,1}$, this confirms that $\mathrm{CPR}$ is a strict subset of the completely stabiliser-preserving channels, and indicates that $\mathcal{R}_{\mathrm{CPR}}$ is not a monotone under stabiliser operations. We also note that the calculated value is larger than the robustness of magic for any single-qubit state, despite $\Lambda_{H}$ being a stabiliser operation. For a circuit containing M uses of the channel $\Lambda_{H}$, the samples required for a CPR simulator would be proportional to $\mathcal{R}_{\mathrm{CPR}}(\Lambda_{H})^{2M}=4^{M}$. But for our simulator, $\Lambda_{H}$ can be simulated efficiently, as $\mathcal{R}_{*}\quantity(\Lambda_{H})^{2M}=1$.

While the above example is quite artificial, a reduction in sample complexity is also achieved for channels where $\mathcal{R}_{*}(\mathcal{E})>1$, but is strictly smaller than $\mathcal{R}_{\mathrm{CPR}}(\mathcal{E})$. Given a circuit decomposed as $\mathcal{E}=\mathcal{E}_{L}\circ\ldots\circ\mathcal{E}_{2}\circ\mathcal{E}_{1}$ , the sample complexity for the CPR simulator would be proportional to $\prod_{j}^{L}\mathcal{R}_{\mathrm{CPR}}(\mathcal{E}_{j})^{2}$. It is always the case that $\mathcal{R}_{*}(\mathcal{E}_{j})\leq\mathcal{R}_{\mathrm{CPR}}(\mathcal{E}_{j})$, so the sample complexity for our simulator will never be greater. Suppose we find that there are $M$ circuit elements such that $\mathcal{R}_{*}(\mathcal{E}_{j})/\mathcal{R}_{\mathrm{CPR}}(\mathcal{E}_{j})\leq k$ for some constant $0<k<1$. Then we would find that using our simulator gives a reduction in sample complexity by a factor of $k^{2M}$. While our simulator sometimes incurs a modest increase in the runtime per sample, this must be weighed against a reduction in sample complexity that is exponential in the number of circuit elements where $\mathcal{R}_{*}(\mathcal{E}_{j})<\mathcal{R}_{\mathrm{CPR}}(\mathcal{E}_{j})$. The obvious next question is whether there are any natural non-trivial examples where this happens. We show in Section 8 that gate sequences subject to amplitude-damping noise provide one such case.

We also note that calculation of optimal $\mathrm{CPR}$ decompositions is only tractable for one- and two-qubit circuit elements, as the three-qubit case already involves a linear program with nearly 93 million variables Bennink et al. (2017). For the most general quantum channels, we encounter a similar problem, as for three-qubit circuit elements, we in principle need to optimise over six-qubit stabiliser states. However, in Appendix E we show that for diagonal channels the problem can be greatly simplified, and the problem becomes tractable for operations on up to five qubits. This allows our algorithm to take advantage of the submultiplicativity of channel robustness; for example for diagonal channels where $\mathcal{R}(\mathcal{E}^{\otimes n})<\mathcal{R}(\mathcal{E})^{n}$, it is advantageous to compose $n$ single-qubit circuit elements together as a single $n$-qubit circuit element, before running the linear program. We will see in Section 8 that this strategy is useful for single-qubit $Z$-rotations.

7.2 Dynamic Monte Carlo

In the previous simulator, all convex optimisations are calculated in advance. However, we have found examples of channels where $\mathcal{C}(\mathcal{E})<\mathcal{R}_{*}(\mathcal{E})$. For such channels, and for any stabiliser state $\rho$, the robustness of the output state $\mathcal{E}(\rho)$ will always be less than the $\ell_{1}$-norm of the decomposition of $\mathcal{E}$ into stabiliser-preserving CPTP channels. Our next simulator takes advantage of this, and we present pseudocode in Algorithm  2.

As in the previous simulator, we can represent the trajectory through the circuit by a vector $\vec{k}$, such that the output of the true quantum circuit would be $\rho=\sum_{\vec{k}}q_{\vec{k}}\sigma_{\vec{k}}$. The major difference is that $q_{\vec{k}}$ cannot be written in the form of equation (62) with $p_{j}$ independent of $\vec{k}$, as the quasiprobabilities for each intermediate decomposition will depend on the stabiliser state sampled in the previous step. Consider what happens for the $j$th circuit element. After step (d)ii. of the algorithm, we have some stabiliser state $\sigma_{\vec{k}_{j}}$, where $\vec{k}_{j}$ labels the trajectory up to the $j$th element. Here we do not assume $\vec{k}_{j}$ is a binary vector; instead the elements of the vector label each pure stabiliser state. After steps (d)iii. and iv. we have a non-stabiliser state $\rho_{\vec{k}_{j}}=(\mathcal{E}_{j}\otimes\mathds{1})(\sigma_{\vec{k}_{j}})$, decomposed as:

[TABLE]

Here the summation is over all $(j+1)$-step trajectories consistent with the previous $j$-step trajectory labelled by $\vec{k}_{j}$. In steps (d)v. and vi. we will choose the stabiliser state $\sigma_{\vec{k}_{j+1}}$ with probability $|q_{\vec{k}_{j+1}}|/\mathcal{R}(\rho_{\vec{k}_{j}})$ and the variable $R$ picks up a factor $\lambda_{\vec{k}_{j+1}}\mathcal{R}(\rho_{\vec{k}_{j}})$, where $\lambda_{\vec{k}_{j+1}}$ is the sign of the corresponding quasiprobability.

The final state after the full sequence of quantum operations may be written:

[TABLE]

The true expectation value for the observable $Z$ is given by:

[TABLE]

A key difference with the previous simulator is that we never explicitly calculate the full distribution $p_{\vec{k}}$. Nevertheless, each time the simulator samples, it produces output $R_{\vec{k}}E_{\vec{k}}$ with probability $p_{\vec{k}}$. These probabilities exactly match the weightings in the above equation, so the simulator is an unbiased estimator. The number of samples required can be again derived using the Hoeffding inequalities, which depend on the maximum possible values of the output of each sample. Each output is bounded by $|R_{\vec{k}}E_{\vec{k}}|\leq|R_{\vec{k}}|\leq\prod_{j}\mathcal{R}(\rho_{\vec{k}_{j}})\leq\prod_{j}\mathcal{C}(\mathcal{E}_{j})$. Therefore, the sample complexity is upper-bounded by order $\prod_{j}\mathcal{C}(\mathcal{E}_{j})^{2}$.

Notice that for every sample, $L$ convex optimisations are performed, as well as $L$ steps involving the algorithm of Fattal et al. Fattal et al. (2004), which has runtime polynomial in the total number of qubits. If $\mathcal{C}(\mathcal{E})=\mathcal{R}_{*}(\mathcal{E})$ then we would simply not use this method so that the only convex optimisations are in the preprocessing. However, $\mathcal{C}(\mathcal{E})$ now determines the sample complexity, so if $\mathcal{C}(\mathcal{E})\ll\mathcal{R}_{*}(\mathcal{E})$ then the dynamic simulator may run much faster than the static simulator; here we have a trade-off of increase in per-sample runtime, versus a possibly exponential reduction in sample complexity. Indeed, since the sample complexity is typically the bottleneck, this approach would lead to significant improvements for some quantum channels. In the following section we investigate which types of channels may lead to an advantage.

8 Numerical results

The numerical results in this section have been produced using code and data files available from the public repository detailed in Ref. Seddon (2019).

8.1 Single-qubit rotation with amplitude damping

Consider the setting discussed in Section 7, where a many-qubit circuit evolution $\mathcal{E}$ is decomposed as a series of few-qubit circuit elements $\mathcal{E}=\mathcal{E}_{L}\ldots\circ\mathcal{E}_{2}\circ\mathcal{E}_{1}$. Many implementations of quantum algorithms can be expected to involve single-qubit rotations about some Pauli axis. In near-term devices, the circuit will be subject to noise. Consider a simple model of a noisy computation where a noise channel $\Lambda$ acts between each unitary gate $\mathcal{U}_{j}$, so the overall channel representing the circuit would be:

[TABLE]

Here we study the simulation cost for a single step in such a computation, comprised of a single-qubit rotation and a noise channel. Note that for intermediate steps in a circuit decomposition such as (74), we have a choice of ordering. We can take the circuit elements to be either $\Lambda\circ\mathcal{U}_{j}$ or $\mathcal{U}_{j}\circ\Lambda$. These choices are equivalent in terms of the output of the simulation, but could lead to different sample complexity depending on the cost function used. We studied circuit elements made up of a single-qubit Pauli $X$-rotation $U(\theta)=\exp(iX\theta)$ composed with an amplitude damping channel $\Lambda_{p}$ with noise parameter $p$, defined by Kraus operators:

[TABLE]

We calculated channel robustness $\mathcal{R}_{*}$, magic capacity $\mathcal{C}$, and the Bennink et al. Bennink et al. (2017) cost function $\mathcal{R}_{\mathrm{CPR}}$, for both $\Lambda_{p}\circ\mathcal{U}(\theta)$ and $\mathcal{U}(\theta)\circ\Lambda_{p}$, and for a range of values of $p$ and $\theta$. For noise $p=0.1$ we see that when the noise channel follows the gate, there is no difference between the three quantities (Figure 3(i)). However, if the noise channel acts before the unitary, both our monotones show a reduced value, whereas $\mathcal{R}_{\mathrm{CPR}}$ increases. This suggests that for this noise model, the better strategy with respect to sample complexity would be to choose the ordering $\mathcal{U}(\theta)\circ\Lambda_{p}$, and use one of our simulators (subject to the caveats mentioned in the previous section). In Figure 3(ii) we show how different levels of noise affect the channel robustness. We do not plot capacity since we find that $\mathcal{R}_{*}(\mathcal{E})=\mathcal{C}(\mathcal{E})$ for this class of operation, up to solver precision. We also compare channel robustness with the Choi state robustness (Figure 3(iii)). We find that $\mathcal{R}\quantity(\Phi_{\mathcal{E}})=\mathcal{R}_{*}(\mathcal{E})$ for $\theta$ up to approximately $\pi/16$, but $\mathcal{R}\quantity(\Phi_{\mathcal{E}})<\mathcal{R}_{*}(\mathcal{E})$ for larger angles.

8.2 Multiqubit phase gates

Recall that from Theorem 3.2 we know $\mathcal{R}\quantity(\Phi_{\mathcal{E}})\leq\mathcal{C}\quantity(\mathcal{E})\leq\mathcal{R}_{*}\quantity(\mathcal{E})$. For the particular class of channel studied above, we saw numerically that $\mathcal{C}\quantity(\mathcal{E})=\mathcal{R}_{*}\quantity(\mathcal{E})$ up to solver precision in all cases investigated, and in the absence of noise the numerical results suggested $\mathcal{R}\quantity(\Phi_{\mathcal{E}})=\mathcal{C}\quantity(\mathcal{E})=\mathcal{R}_{*}\quantity(\mathcal{E})$. Moreover, we know from Theorem 3.2 that all measures are equal for gates from the third level of the Clifford hierarchy. Under what conditions does this equality persist for multi-qubit operations? As explained earlier, to calculate each of our quantities for $n$-qubit channels, in general we must solve an optimisation problem over all $2n$-qubit stabiliser states. Since this problem is only tractable for up to $5$-qubit states, in practice we are limited to studying two-qubit channels, in the most general case. However, it turns out the problem can be greatly simplified for certain types of operation. In particular, Appendix E shows how the problem size can be reduced for channels diagonal in the computational basis, using a representation of stabiliser states in terms of affine spaces over binary vectors due to Dehaene and De Moor Dehaene and De Moor (2003); Gross and Van den Nest (2008). MATLAB code to calculate our measures for this reduced problem is provided in the public repository given in Ref. Seddon (2019). This allows us to calculate values for diagonal operations on up to 5 qubits, which we present here.

As a special case we consider multicontrol phase gates of the form:

[TABLE]

where $n$ denotes the number of qubits. We note that by convention, controlled-phase gates typically apply the phase to the all-one state $\ket{1^{n}}$, where $1^{n}=(1,\ldots,1)^{T}$, but the form given above is Clifford-equivalent to the conventional version, and is more convenient for the techniques used in Appendix E. The family includes familiar gates such as $CZ$ ($t=0$, $n=2)$, $CCZ$ ($t=0$, $n=3$), multicontrol-$S$ ($t=1)$ and multicontrol-$T$ ($t=2$).

The main findings were that the inequalities are tight for the $n=2$ and $n=3$ cases, but that this does not persist for larger system sizes (Figure 4). The $t=0$ case (the family of multicontrol-$Z$ gates) turns out to be a special case (Figure 4, left panel). Here we find equality for all three quantities up to $n=4$. For the $t=0$, $n=5$ case, $\mathcal{R}\quantity(\Phi_{M_{0,5}})=\mathcal{C}\quantity(M_{0,5})$ holds, but $\mathcal{R}_{*}\quantity(M_{0,5})$ is strictly greater than both. Note also that for $t=0$, all three quantities increase with each increment in $n$.

The families of gates with $t>0$ follow a pattern qualitatively similar to each other. The results for the $t=1$ (multicontrol-$S$) and $t=2$ (multicontrol-$T$) cases are shown in the middle and right panels of Figure 4. For $n=4$, $t>0$, the same situation holds as for $n=5$, $t=0$, as we find $\mathcal{R}\quantity(\Phi_{M_{t,4}})=\mathcal{C}\quantity(M_{t,4})<\mathcal{R}_{*}\quantity(M_{t,4})$. At $n=5$, all three quantities separate. In contrast with the multicontrol-$Z$, we see that $\mathcal{R}\quantity(\Phi_{M_{t,n}})$ decreases as we go from four to five qubits, while $\mathcal{C}\quantity(M_{t,n})$ levels off. We see similar behaviour for all non-zero values of $t$ investigated numerically. Our current techniques limit us to five-qubit gates, but we have reason to believe that the capacity will remain level for $n>5$, and we make the following conjecture, which we justify more fully in Appendix E.4.

Conjecture 8.1**.**

For any fixed $t$, the maximum increase in robustness of magic for $M_{t,n}$ is achieved at some finite number of qubits $n=K$ by acting on the state $\ket{+}^{\otimes K}$. Therefore $\mathcal{C}\quantity(M_{t,n})=\mathcal{R}\quantity(M_{t,K}\ket{+}^{\otimes K})$ for all $n\geq K$.

We also numerically investigated the robustness of diagonal unitaries

[TABLE]

with $\theta_{x}$ chosen uniformly at random. We were particularly interested in understanding when $\mathcal{R}(\Phi_{\mathcal{E}})\leq C(\mathcal{E})\leq\mathcal{R}_{*}(\mathcal{E})$ is tight or loose. In Figure  5 we compare the Choi robustness with the channel robustness. For every 2-qubit gate tested we observed that $\mathcal{R}(\Phi_{\mathcal{E}})=\mathcal{R}_{*}(\mathcal{E})$ up to numerical precision. Whereas, for 3 and 4 qubit gates we typically saw that $\mathcal{R}(\Phi_{\mathcal{E}})<\mathcal{R}_{*}(\mathcal{E})$, though the gap is not often large. While the difference is slight for a single gate, these quantities influence the rate of exponential scaling when considering $N$ uses of such a unitary and will lead to a large gap for modest $N$.

We also compared the Choi robustness with the magic capacity but do not plot this data as it was equal for every random instance we observed. This is curious since in Figure  4 we clearly see that there do exist diagonal gates, the multicontrol phase gates, for which there is a gap between the Choi robustness and the magic capacity. While such gates exist, our random sampling of diagonal gates does not tend to provide such examples. We discuss this further in Appendix E.4.

Finally, we are also interested in the normalised channel robustness for single-qubit gates $U$, defined as $\quantity[\mathcal{R}_{*}(U^{\otimes n})]^{1/n}$. This allows us to quantify the per-gate savings in sample complexity that can be achieved by grouping single-qubit rotations in $n$-qubit blocks. In Figure 6 we present results for qubit $Z$-rotations $U=\exp[iZ\theta]$, up to four qubits. We find that strict submultiplicativity is observed for all values of $\theta$, with significant reductions between the $n=2$ and $n=4$ cases for a wide range of angles.

9 Conclusion

We have presented two new magic monotones for general quantum channels: the magic capacity $\mathcal{C}$, which quantifies the ability of a channel to generate magic, and the channel robustness $\mathcal{R}_{*}$, which is related to finding the minimal quasiprobability decomposition of a channel into stabiliser-preserving CPTP maps. Each of these monotones is directly related to the sample complexity for an associated Monte Carlo-type classical simulation algorithm. We found that for certain quantum channels, our static simulator would lead to a exponentially better sample complexity as compared to that for the algorithm due to Bennink et al. Bennink et al. (2017). In particular we found a reduction in sample complexity for the case of a sequence of single-qubit rotations subject to amplitude-damping noise. Since our decompositions can be calculated for up to five qubits in the case of diagonal operations, our static simulator is also able to take advantage of the submultiplicativity of channel robustness under tensor product and composition.

For some channels, further improvements in sample complexity are possible using a different simulator that is related to the capacity. That is, we found the capacity can be strictly less than channel robustness for certain multi-qubit entangling gates. This simulator has to introduce on-the-fly convex optimisation, however, so each sample will be more difficult to obtain.

For simulation of realistic quantum devices involving many qubits, one would need to decompose the circuit into a sequence of operations on smaller number of qubits, as described in Section 7, in order that the associated optimisation problems are tractable. It is a non-trivial problem to decide what is the optimal way to block together the few-qubit operations making up a given many-qubit circuit: we saw in Section 8 that the ordering of operations can make a difference to the sample complexity. We leave this problem for a future work.

Since our monotones are submultiplicative under tensor product and compositions, it is generally preferable to combine circuit elements where possible. In practice, to calculate our monotones, circuit elements can involve up to two qubits for the most general case, or at most five qubits for diagonal operations, where we can make use of the techniques described in Appendix E. In Ref. Heinrich and Gross (2019), Heinrich and Gross show that robustness of magic can be calculated for up to 10 copies of the resource states employed in the standard magic state model of fault-tolerant quantum computation. Their methods rely on both the permutation symmetry due to having multiple copies, and the stabiliser symmetries of the magic states considered. Another direction for future work could therefore be to investigate whether (and under what conditions) similar techniques can be applied to the channel picture to increase the number of qubits that can be involved in each circuit element.

Acknowledgements.

This work was supported by the Engineering and Physical Sciences Research Council [grant numbers EP/P510270/1 (J.R.S.) and EP/M024261/1 (E.T.C.)] . We would like to thank Mark Howard for many useful discussions and assistance with linear optimisation. We would also like to thank Ryuji Takagi for interesting discussions on stabiliser-preserving operations. Thank you to Hakop Pashayan and Dan Browne for helpful discussions on classical simulation and the discrete Wigner function. Finally we would like to thank the anonymous reviewers for their comments and suggestions on the manuscript.

Appendix A Properties of robustness of the Choi state

Here we confirm that the robustness of the Choi state, $\mathcal{R}\quantity(\Phi)$ has the properties convexity and submultiplicativity under tensor product. We then give an example to show that it is not submultiplicative under composition.

Convexity: This follows immediately from convexity of robustness of magic. Consider a real linear combination of $n$-qubit channels: $\mathcal{E}=\sum_{k}q_{k}\mathcal{E}_{k}$. The Choi state for $\mathcal{E}$ is:

[TABLE]

where in the last line we identified $\quantity(\mathcal{E}_{k}\otimes\mathds{1}_{n})\outerproduct{\Omega}{\Omega}$ as the Choi state for $\mathcal{E}_{k}$. Then by convexity of robustness of magic:

[TABLE]

which shows $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$ is convex in $\mathcal{E}$.

Submultiplicativity under tensor product: The maximally entangled state $\ket{\Omega_{n+m}}^{AA^{\prime}|BB^{\prime}}$ as defined by equation (15) in the main text can be factored as $\ket{\Omega_{n+m}}^{AA^{\prime}BB^{\prime}}=\ket{\Omega_{n}}^{A|B}\ket{\Omega_{m}}^{A^{\prime}|B^{\prime}}$. So the Choi state for a channel $\mathcal{E}^{AA^{\prime}}=\mathcal{E}^{A}\otimes\mathcal{E}^{\prime A^{\prime}}$, where $\mathcal{E}^{A}$ and $\mathcal{E}^{\prime A^{\prime}}$ are respectively $n$-qubit and $m$-qubit channels, can be written:

[TABLE]

Then by submultiplicativity of robustness of magic for states, we have:

[TABLE]

which is the desired property.

Failure of submultiplicativity under composition: Let $\mathcal{E}_{1}$ be the single-qubit $Z$-reset channel defined by Kraus operators $\quantity{\outerproduct{0}{0},\outerproduct{0}{1}}$, and let $\mathcal{E}_{2}$ be the conditional channel defined by $\quantity{\outerproduct{T}{0},\outerproduct{1}{1}}$, where $\ket{T}=T\ket{+}$. These channels respectively have Choi states $\Phi_{\mathcal{E}_{1}}=\outerproduct{0}{0}\otimes\frac{\mathds{1}}{2}$, and $\Phi_{\mathcal{E}_{2}}=\frac{1}{2}\quantity(\outerproduct{T0}{T0}+\outerproduct{11}{11})$, with robustness of magic $\mathcal{R}\quantity(\Phi_{\mathcal{E}_{1}})=1$ and $\mathcal{R}\quantity(\Phi_{\mathcal{E}_{2}})\approx 1.207$.

The composed channel $\mathcal{E}_{2}\circ\mathcal{E}_{1}$ has a Kraus representation $\quantity{\outerproduct{T}{0},\outerproduct{T}{1}}$, and so has a Choi state $\Phi_{\mathcal{E}_{2}\circ\mathcal{E}_{1}}=\outerproduct{T}{T}\otimes\frac{\mathds{1}}{2}$, with $\mathcal{R}\quantity(\Phi_{\mathcal{E}_{2}\circ\mathcal{E}_{1}})\approx 1.414>\mathcal{R}\quantity(\Phi_{\mathcal{E}_{2}})\mathcal{R}\quantity(\Phi_{\mathcal{E}_{1}})$. So it is not the case that the robustness of the Choi state is submultiplicative under composition.

More intuitively, such counterexamples arise for channels $\mathcal{E}$ where the stabiliser state $\ket{\phi_{*}}$ that results in maximal final robustness $\mathcal{R}\quantity[\quantity(\mathcal{E}\otimes\mathds{1}_{n})\outerproduct{\phi_{*}}{\phi_{*}}]$ is not the maximally entangled state $\ket{\Omega_{n}}$, as then we can always boost the output robustness by using a stabiliser-preserving operation to prepare $\ket{\phi_{*}}$ before applying $\mathcal{E}$.

Appendix B Properties of channel robustness

Faithfulness: Suppose $\mathcal{E}$ is an $n$-qubit CPTP map. There are two cases:

(i) $\mathcal{E}\in\mathrm{SP}_{n,n}$. In this case, $\Phi_{\mathcal{E}}$ is itself a mixed stabiliser state, and since $\mathcal{E}$ is trace-preserving it satisfies $\Tr(\Phi_{\mathcal{E}})=\mathds{1}_{n}/2^{n}$. So $\Phi_{\mathcal{E}}$ is already trivially a decomposition of the correct form, with $p=0$, so that $R_{*}(\mathcal{E})=1+2p=1$.

(ii) $\mathcal{E}\notin\mathrm{SP}_{n,n}$. Then by faithfulness of robustness of magic (Lemma 4.2 in the main text), $\Phi_{\mathcal{E}}$ has $\mathcal{R}\quantity(\Phi_{\mathcal{E}})>1$. Since the definition of $\mathcal{R}_{*}$ is a restriction of $\mathcal{R}\quantity(\Phi)$, it must be the case that $\mathcal{R}\quantity(\Phi_{\mathcal{E}})\leq\mathcal{R}_{*}\quantity(\mathcal{E})$. Therefore $\mathcal{R}_{*}(\mathcal{E})>1$.

Convexity:

Suppose we have a set of Choi states $\Phi_{\mathcal{E}_{j}}$ corresponding to channels $\mathcal{E}_{j}$, with optimal decompositions:

[TABLE]

where each $\rho_{j_{\pm}}$ separately satisfies the condition $\Tr_{A}(\rho_{\pm})=\frac{\mathds{1}_{n}}{2^{n}}$, so that $\mathcal{R}_{*}(\mathcal{E}_{j})=1+2p_{j}$. Now take a real linear combination of such channels:

[TABLE]

where $P$ is the set of indices such that $q_{j}\geq 0$, and $N$ is the set such that $q_{k}<0$. We assume that $\sum_{i}q_{i}=1$ so that the trace of $\Tr(\Phi_{\mathcal{E}})=1$. Then the corresponding Choi state for channel $\mathcal{E}$ is:

[TABLE]

Note that the terms inside the brackets all have positive coefficients, hence we can interpret as non-normalised mixtures over stabiliser states. To normalise them we can define:

[TABLE]

and

[TABLE]

Then writing:

[TABLE]

one can check that:

[TABLE]

This allows us to rewrite the Choi state as: $\Phi_{\mathcal{E}}=(1+\widetilde{p})\widetilde{\rho}_{+}-\widetilde{p}\widetilde{\rho}_{-}.$ Since $\widetilde{\rho}_{\pm}$ are convex mixtures over stabiliser states satisfying $\Tr_{A}(\rho_{j\pm})=\frac{\mathds{1}_{n}}{2^{n}}$, they must satisfy the same condition. We also know that $\widetilde{p}\geq 0$, so it is clear that the decomposition is in the form required for the definition of $\mathcal{R}_{*}$, except that it is not necessarily optimised to minimise $1+2\widetilde{p}$. So we have:

[TABLE]

which gives us the required result.

Invariance under tensor with identity: In Section 5 of the main text we saw that $\mathcal{R}_{*}\quantity(\mathcal{E}^{A}\otimes\mathds{1})\leq\mathcal{R}_{*}\quantity(\mathcal{E}^{A})$. We now complete the proof that $\mathcal{R}_{*}\quantity(\mathcal{E}^{A}\otimes\mathds{1})=\mathcal{R}_{*}\quantity(\mathcal{E}^{A})$ by showing that $\mathcal{R}_{*}\quantity(\mathcal{E}^{A})\leq\mathcal{R}_{*}\quantity(\mathcal{E}^{A}\otimes\mathds{1})$.

Consider an optimal decomposition for $\Phi_{\mathcal{E}^{A}\otimes\mathds{1}_{m}}=(1+p^{\prime})\rho^{\prime}_{+}-p^{\prime}\rho^{\prime}_{-}$, such that $\mathcal{R}_{*}(\mathcal{E}^{A}\otimes\mathds{1}_{m})=1+2p^{\prime}$, where $\Tr_{AA^{\prime}}(\rho_{\pm})=\mathds{1}_{n+m}/2^{n+m}$. Here we do not assume that $\rho^{\prime}_{\pm}$ are products across the partition $AB|A^{\prime}B^{\prime}$, as was the case in equation (34) in the main text. However, we have just seen that $\Phi_{\mathcal{E}^{A}\otimes\mathds{1}_{m}}$ can be written as a product, so that by tracing out systems $A^{\prime}B^{\prime}$ we obtain:

[TABLE]

Partial trace of a stabiliser state remains a stabiliser state, so this is a stabiliser decomposition. We just need to check that the partial trace condition holds, so we want to show:

[TABLE]

but this is clearly the case from the fact that $\rho^{\prime}_{\pm}$ were constrained such that $\Tr_{AA^{\prime}}(\rho_{\pm})=\mathds{1}_{n+m}/2^{n+m}$. Hence again we have a valid, not necessarily optimal decomposition and:

[TABLE]

Combining with the inequality $\mathcal{R}_{*}(\mathcal{E}^{A}\otimes\mathds{1}^{B})\leq\mathcal{R}_{*}(\mathcal{E}^{A})$, shown in the main text, we obtain the equality:

[TABLE]

Appendix C Optimisation problem for channel robustness

In Howard and Campbell Howard and Campbell (2017), the optimisation problem for calculating robustness of magic for states was cast as follows:

[TABLE]

where $\vec{q}$ is a vector of coefficients, $\vec{b}$ is the vector of Pauli expectation values for the target state $\Phi_{\mathcal{E}}$, and $A$ is a matrix whose columns are the Pauli vectors for the stabiliser states. For $n$-qubit channels, we have $2n$-qubit Choi states, so the number of generalised Paulis is $N_{P}=4^{2n}$, and the number of stabiliser states is $N_{S}=2^{2n}\prod_{j=1}^{2n}(2^{j}+1)$ Howard and Campbell (2017). Then $\vec{b}$ has $N_{P}$ entries, $\vec{q}$ has $N_{S}$ entries, and the dimension of $A$ is $(N_{P}\times N_{S})$. From this construction we can recover optimal decompositions of the form: $\Phi_{\mathcal{E}}=\sum_{j}q_{j}\outerproduct{\phi_{j}}{\phi_{j}}$, where $\sum_{j}q_{j}=1$ and $\ket{\phi_{j}}$ are the pure stabiliser states.

We want to restrict the problem to decompositions of the form:

[TABLE]

where $p\geq 0$ and $\rho_{\pm}$ correspond to trace-preserving channels, and can in general be mixed. Rather than enumerating all the extreme points of the set of stabiliser states corresponding to maps in $\mathrm{SP}_{n,n}$, it is more convenient to retain the same $A$ matrix and modify the constraints. We still need to start from a finite set of extreme points, i.e. pure stabiliser states, so first rewrite as:

[TABLE]

where $q_{j_{+}}$ are the positive quasiprobabilities, $q_{j_{-}}$ are the negative quasiprobabilities, and $p_{j_{\pm}}=|q_{j_{\pm}}|$. In the Pauli vector picture we can write this as $\vec{b}=A\vec{p}_{+}-A\vec{p}_{-}$, where all the entries of $\vec{p}_{\pm}$ are non-negative. We define a new variable vector $\vec{p}$ which will have twice the length of the previous $\vec{q}$, i.e. $2N_{S}$ entries:

[TABLE]

and define a new $(N_{P}\times 2N_{S})$ matrix $A^{\prime}$ in block form, $A^{\prime}=\begin{pmatrix}A&-A\end{pmatrix}$. Then we have:

[TABLE]

So now we need to minimise $\|\vec{p}\|_{1}=\sum_{j}p_{j}$ subject to $A^{\prime}\vec{p}=\vec{b}$ and $\vec{p}\geq 0$.

Next, we need the trace-preserving condition. Provided $\mathcal{E}$ is CPTP, if one part of the decomposition is trace-preserving, then the other will be as well, so we only need enforce the constraint on one of $\rho_{+}$ or $\rho_{-}$. Assume that we check $\rho_{+}$. The condition for a Choi state $\Phi^{AB}=\mathcal{E}^{A}\otimes\mathds{1}^{B}(\ket{\Omega}\bra{\Omega}^{AB})$ to be trace-preserving is:

[TABLE]

where $d$ is the dimension of the subsystem. We need to convert this to a constraint on the vector $\vec{b}_{+}$ corresponding to $\phi_{+}$, which is given by $\vec{b}_{+}=A\vec{p}_{+}$. First, note that all Paulis are traceless except for the identity $P_{0}=\mathds{1}$, so for the maximally mixed state:

[TABLE]

so if the first entry in a Pauli vector is always $\langle\mathds{1}\rangle$, the maximally mixed state has Pauli vector:

[TABLE]

where $\vec{0}$ is the zero vector. However, we need this to hold just for the reduced state on $B$ rather than the full Pauli vector. Consider that if the whole state is written $\Phi^{AB}=\sum_{j,k}r_{j,k}P_{j}\otimes P_{k}$ for some set of coefficients $r_{j,k}$, then the expectation values are given by:

[TABLE]

The reduced state is:

[TABLE]

and the entries of the reduced Pauli vector will be:

[TABLE]

So for condition (A103) to hold for the reduced state on $B$, we combine equations (A104) and (A108) to get:

[TABLE]

That is, we just need to look at the entries of $\vec{b}_{+}$ corresponding to Paulis of the form $\mathds{1}\otimes P_{j}$. These should all be zero except the first entry, which corresponds to $\langle\mathds{1}\otimes\mathds{1}\rangle$. Note that $\vec{b}_{+}=A\vec{p}_{+}$ will in general not be normalised, but this does not matter, since we are only interested in whether or not entries are zero. We can use a binary matrix $M$ to pick out the values of interest. As an example we consider the two-qubit case, and assume that the entries are ordered as:

[TABLE]

Here, we are only interested in the 2nd, 3rd and 4th entries. We form a new vector $\vec{c}$ by left multiplying with $M$:

[TABLE]

Then the condition we need is just $\vec{c}=0$. To convert this to a condition on the $2N_{S}$-entry variable $\vec{p}=\begin{pmatrix}\vec{p}_{+}\\ \vec{p}_{-}\end{pmatrix}$, we first pad $A$ with zeroes: $A_{+}=\begin{pmatrix}A&\overline{0}\end{pmatrix}$, where $\overline{0}$ is the $(N_{P}\times N_{S})$ zero matrix. We then have:

[TABLE]

so that $\vec{c}=M\vec{b}_{+}=MA_{+}\vec{p}$. Therefore, our condition for trace-preserving $\rho_{+}$ is $MA_{+}\vec{p}=0$. We can therefore specify the new optimisation problem as:

[TABLE]

where $A^{\prime}=\begin{pmatrix}A&-A\end{pmatrix}$, and $A_{+}=\begin{pmatrix}A&\overline{0}\end{pmatrix}$, with $A$ and $\vec{b}$ having the same definitions as previously, $\overline{0}$ is the zero matrix with dimension the same as $A$, and with $M$ being the binary matrix that picks out the $\langle\mathds{1}\otimes P_{j}\rangle$ entries from the vector $A_{+}\vec{p}$. Most of this is straightforward to implement. The step that requires some care is in correctly constructing the matrix $M$, as it will depend on the choice of ordering of Pauli operators in the construction of $A$ and $\vec{b}$. If the $B$ subsystem has $n$ qubits, then we will need to constrain $4^{n}-1$ non-trivial $\langle\mathds{1}\otimes P_{j}\rangle$ expectation values to zero, so $M$ should have dimension $((4^{n}-1)\times N_{P})$. If the Paulis are ordered as in the example given above for $2$-qubit Choi states, then the construction is just $M=\begin{pmatrix}\vec{0}&\mathds{1}^{\prime}&\vec{0}&\cdots&\vec{0}\end{pmatrix}$, where $\mathds{1}^{\prime}$ is the $((4^{n}-1)\times(4^{n}-1))$ identity, and $\vec{0}$ denotes a column of zeroes. We have implemented this linear program in MATLAB, using the convex optimisation package CVX Grant and Boyd (2014), and have made the code available from the repository Ref. Seddon (2019).

Appendix D Properties of magic capacity

Faithfulness: For any $n$-qubit stabiliser-preserving CPTP channel $\Lambda$, if $\rho\in\mathrm{STAB}_{2n}$ is a stabiliser state, then $(\Lambda\otimes\mathds{1}_{n})\rho$ is also a stabiliser state. So by the faithfulness of robustness of magic, $\mathcal{R}\quantity((\Lambda\otimes\mathds{1}_{n})\rho)=1$ for any input stabiliser state $\rho\in\mathrm{STAB}_{2n}$, and $\mathcal{C}\quantity(\Lambda)=1$.

Suppose instead that $\mathcal{E}$ is non-stabiliser-preserving, but still CPTP. Then there exists at least one stabiliser state $\rho\in\mathrm{STAB}_{2n}$ such that $(\mathcal{E}\otimes\mathds{1})\rho$ is a normalised state, but not a stabiliser state. Then by faithfulness of $\mathcal{R}$ when applied to states, $\mathcal{R}\quantity((\mathcal{E}\otimes\mathds{1})\rho)>1$, and so $\mathcal{C}\quantity(\mathcal{E})>1$.

Convexity: Suppose we have a real linear combination of $n$-qubit CPTP maps $\mathcal{E}_{k}$:

[TABLE]

There exists some optimal stabiliser state $\rho_{*}$ that achieves $\mathcal{C}\quantity(\mathcal{E})=\mathcal{R}\quantity(\mathcal{E}\otimes\mathds{1}\quantity(\rho_{*}))$. Then

[TABLE]

where the last line follows by convexity of the robustness of magic. But each robustness $\mathcal{R}\quantity(\quantity(\mathcal{E}_{k}\otimes\mathds{1}_{n})\rho_{*})$ can be no larger than $\mathcal{C}\quantity(\mathcal{E}_{k})$. So we have:

[TABLE]

Appendix E Calculating monotones for diagonal channels

E.1 Reducing the problem size

As mentioned earlier, the size of the optimisation problem for calculating our monotones (as well as $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$) quickly becomes prohibitively large for $n$-qubit states, since the number of stabiliser states increases super-exponentially with $n$ (Table 2).

The issue is even worse than it first appears, since for an $n$-qubit channel we must in general consider $2n$-qubit stabiliser states. Direct calculation of either monotone is impractical for $n$-qubit channels with $n>2$. This difficulty is aggravated when calculating the capacity as in principle we have to repeat the optimisation for every $\quantity(\mathcal{E}\otimes\mathds{1}_{n})\outerproduct{\phi}{\phi}$ such that $\ket{\phi}\in\mathrm{STAB}_{2n}$. In some cases we can ameliorate these problems by looking for Clifford gates that commute with the channel of interest. Here we consider the case where $\mathcal{E}$ is a diagonal channel, meaning it has a Kraus representation where each Kraus operator is diagonal in the computational basis. This of course includes diagonal unitaries as a special case. One could likely reduce the problem size further by exploiting symmetries of channels using techniques similar to those used in Ref. Heinrich and Gross (2019), but we will not consider this strategy here.

It is straightforward to see how the problem can be simplified for calculating $\mathcal{R}\quantity(\Phi_{\mathcal{E}})$ and $\mathcal{R}_{*}(\mathcal{E})$. If $\mathcal{E}$ is diagonal, the operation $\mathcal{E}\otimes\mathds{1}_{n}$ commutes with any sequence of CNOTs targeted on the last $n$ qubits. But the maximally entangled state $\ket{\Omega_{n}}$ can be written:

[TABLE]

Here $U_{C}=\otimes_{j=1}^{n}U_{j}$, where $U_{j}$ is the CNOT controlled on qubit $j$ and targeted on qubit $n+j$. By the monotonicity of robustness of magic, we immediately see that:

[TABLE]

For the channel robustness we would like to decompose $\mathcal{E}\quantity(\outerproduct{+}{+}^{\otimes n})$ in terms of states $\rho_{\pm}\in\mathrm{STAB}_{n}$, but need to take care that the trace condition $\Tr_{A}(\rho^{\prime}_{\pm})=\mathds{1}_{n}/2^{n}$ is satisfied for the equivalent $2n$-qubit Choi states $\rho^{\prime}_{\pm}$. In Appendix E.2 we show that the criterion is satisfied provided all diagonal elements of $\rho_{\pm}$ are equal to $1/2^{n}$. So for diagonal channels we can write:

[TABLE]

So calculation of $R_{*}(\mathcal{E})$ and $\mathcal{R}(\Phi_{\mathcal{E}})$ is tractable up to five qubits provided $\mathcal{E}$ is diagonal. We will see below in Section E.3 that this is also true for the magic capacity.

E.2 Trace condition for diagonal channels

Consider that the Choi state for a diagonal channel has a decomposition

[TABLE]

where $U_{C}=\otimes_{j=1}^{n}U_{j}$ is the tensor product of CNOTs $U_{j}$ that are controlled on the $j$th qubit and targeted on the $n+j$th. Then

[TABLE]

where $\rho^{\prime}_{\pm}$ are still stabiliser states since $U_{C}$ is Clifford. Now consider the stabiliser-preserving channel $\mathds{1}_{n}\otimes\Lambda$ that resets the last $n$ qubits to $\outerproduct{0}{0}^{\otimes n}$. Applying this to both sides of equation (A121) we get a new decomposition

[TABLE]

Then referring back to equation (A120), we obtain $\rho_{\pm}=U_{C}\quantity(\rho^{\prime\prime}_{+}\otimes\outerproduct{0}{0}^{\otimes n})U_{C}^{\dagger}$. So, the trace-preserving condition becomes:

[TABLE]

where $\ket{x}$ are the computational basis states on subsystem $A$. Recalling that $U_{C}$ can be written as a tensor product of CNOTs $U_{C}=\otimes_{j=1}^{n}U_{j}$ one can check that this equation can be written:

[TABLE]

Therefore, the decomposition corresponds to a pair of trace-preserving channels provided that all diagonal elements of $\rho^{\prime\prime}_{\pm}$ are equal to $1/2^{n}$.

For a given diagonal channel, there always exists a decomposition that satisfies these conditions and has $\ell_{1}$-norm equal to the channel robustness as defined for the full Choi state. We do not give the full proof here, but sketch the argument. Given any diagonal channel $\mathcal{E}$ decomposition of the full Choi state $\Phi_{\mathcal{E}}=(1+p)\rho_{+}-p\rho_{-}$ satisfying the trace condition, one can always find a new decomposition $\Phi_{\mathcal{E}}=(1+p)\Lambda(\rho_{+})-p\Lambda(\rho_{-})$ where $\Lambda(\rho_{\pm})$ still satisfy $\Tr_{A}(\Lambda(\rho_{\pm}))$, but are now the Choi states for diagonal channels. The map $\Lambda$ used to obtain this decomposition is in effect an error correction circuit that takes general stabiliser Choi states to the subspace corresponding to the diagonal channels. Specifically, we note that the Choi states for diagonal maps $\mathcal{T}$ have the general form:

[TABLE]

In general $c_{j,k}$ can be complex or zero, but terms on the diagonal are constrained. In particular, trace-preserving diagonal channels cannot change the weight of particular computational basis states, so the probability distribution for computational basis states will be the same as for $\ket{\Omega}$:

[TABLE]

The circuit $\Lambda$ is defined by the following steps. For each $j$ from $1$ to $n$:

Perform a parity measurement ($Z\otimes Z$) between qubits $j$ and $n+j$. 2. 2.

If even parity (+1 outcome), do nothing. If odd parity (-1 outcome), perform an $X$ gate on qubit $j$.

This stabiliser-preserving channel leaves Choi states for diagonal maps (and crucially, the target Choi state $\Phi_{\mathcal{E}}$) invariant, but updates general Choi states to have the form (A126). One can check that the circuit preserves the property $\Tr_{A}(\rho_{\pm})=\mathds{1}_{n}/2^{n}$. We then obtain a decomposition in the desired form:

[TABLE]

Where $\Lambda(\rho_{\pm})=(\mathcal{E}_{\pm}\otimes\mathds{1})\outerproduct{\Omega}{\Omega}$ are Choi states for $n$-qubit diagonal channels $\mathcal{E}_{\pm}$. But as described above, the CNOT sequence $U_{C}$ commutes with diagonal channels acting on the first $n$ qubits, so we can obtain $n$-qubit representatives of these channels:

[TABLE]

Discarding the last $n$ qubits we obtain the desired $n$-qubit decomposition:

[TABLE]

E.3 Magic capacity in the affine space picture

In this section we will make use of the formalism due to Dehaene and De Moor, in which stabiliser states are cast in terms of affine spaces and quadratic forms over binary vectors Dehaene and De Moor (2003); Gross and Van den Nest (2008), to prove the following theorem:

Theorem A1** (Capacity for diagonal operations).**

Suppose the $n$-qubit channel $\mathcal{E}_{D}$ is diagonal. Let:

[TABLE]

where $x\in\mathbb{F}_{2}^{n}$ are binary vectors and $\mathcal{K}\subseteq\mathbb{F}_{2}^{n}$ is an affine space. Then:

[TABLE]

That is, given an $n$-qubit channel $\mathcal{E}$, provided the channel is diagonal, the capacity $\mathcal{C}\quantity(\mathcal{E})$ may be calculated by optimisation over only the $n$-qubit states $\ket{\mathcal{K}}$ as defined in equation (A131), rather than over all $2n$-qubit stabiliser states.

We first review the formalism of Ref. Dehaene and De Moor (2003). Computational basis states $\ket{x}$ can be labelled by binary column vectors $x=\quantity(x_{1},\ldots,x_{n})^{T}\in\mathbb{F}_{2}^{n}$, so that $x_{j}\in\quantity{0,1}$ relates to the $j$th qubit. Any pure $n$-qubit stabiliser state may be written:

[TABLE]

where $\mathcal{K}\subseteq\mathbb{F}_{2}^{n}$ is an affine space, $d$ is some fixed binary vector, and $q(x)$ has the form:

[TABLE]

Here $Q$ is a binary, strictly upper triangular matrix, $\lambda$ is a vector, and addition is modulo 2. Conversely, any state that can be written in this way is a stabiliser state.

An affine space $\mathcal{K}$ is a linear subspace $\mathcal{L}$ shifted by some constant binary vector $h$, modulo 2: $\mathcal{K}=\mathcal{L}+h$. Every affine space is related in this way to exactly one linear subspace, and the dimension $k=\dim(\mathcal{K})$ of an affine space means the dimension of the corresponding subspace. Instead of enumerating all elements of an affine space, we can specify it by a shift vector $h$ and an $n\times k$ matrix where each column is one of the generators of the corresponding linear space:

[TABLE]

We have freedom in our choice of $k$ independent generators, and we can transform between equivalent generating sets by adding any two columns of $G$. We are also free to swap any two columns. A general transform between generating sets can therefore be represented by an invertible matrix $S$ of dimension $k\times k$, multiplying on the right $G\longrightarrow GS$.

Any non-trivial linear transformation of the affine space can be fully specified by the transformation of the generators and the shift vector. In particular, we can represent the action of a single CNOT by multiplication on the left by a matrix $C$. If the CNOT has control qubit $j$ and target qubit $k$, then $C$ has 1s on the diagonal, a 1 in the $j$th element of the $k$th row, and zeroes everywhere else. A sequence for a $2n$-qubit system, in which CNOTs are always controlled on the first $n$ qubits, and targeted on the last $n$ qubits can be represented in block form:

[TABLE]

where each block has dimension $n\times n$, and $M$ can be any binary matrix. We use this formalism to prove the following lemma, which leads directly to Theorem A1:

Lemma A1** (Equivalences for diagonal channels).**

Suppose $\mathcal{E}_{D}$ is a diagonal CPTP channel. Then:

All input stabiliser states with the same affine space $\mathcal{K}$ result in the same final robustness:

[TABLE] 2. 2.

Given a $2n$-qubit state $\ket{\phi}\in\mathrm{STAB}_{2n}$ , there exists some $n$-qubit $\ket{\phi^{\prime}}\in\mathrm{STAB}_{n}$ such that:

[TABLE]

Proof.

We first prove statement 1. Since robustness of magic is invariant under Clifford unitaries, we need to show that there exists a Clifford unitary $U$ that converts $(\mathcal{E}_{D}\otimes\mathds{1})\outerproduct{\mathcal{K},q,d}{\mathcal{K},q,d}$ to $(\mathcal{E}_{D}\otimes\mathds{1})\outerproduct{\mathcal{K},q^{\prime},d^{\prime}}{\mathcal{K},q^{\prime},d^{\prime}}$. A suitable choice for $U$ is one such that $U\ket{\phi_{\mathcal{K},q,d}}=\ket{\phi_{\mathcal{K},q^{\prime},d^{\prime}}}$, and, crucially, that commutes with the channel $\mathcal{E}_{D}$. Since $\mathcal{E}_{D}$ is given to be diagonal, any diagonal Clifford $U$ will suffice. The affine space $\mathcal{K}$ remains unchanged, so we only need show there is always a diagonal Clifford that maps $q\to q^{\prime}$ and $d\to d^{\prime}$ for any $q$, $q^{\prime}$, $d$ and $d^{\prime}$. That this is always possible is perhaps already evident from Ref. Dehaene and De Moor (2003), but for completeness we give the argument here.

We can convert $d$ to $d^{\prime}$ using appropriately chosen $S_{j}$ gates, meaning the gate $\mathrm{diag}(1,i)$ acting on the $j$th qubit. Consider the action of $S_{j}$ on a basis vector:

[TABLE]

If we define basis vector $e_{j}$ so that it has $1$ in the $j$th position and zeroes elsewhere, we can write the action of $S_{j}$ as:

[TABLE]

Note that the form of this equation is independent of the value of $x$, so we can write:

[TABLE]

So, we can flip any bit of $d$ by applying the correct $S$ gate. The quadratic form $q(x)$ is left unchanged.

Now consider $q(x)=x^{T}Qx+\lambda^{T}x$, which we must convert to some other ${q^{\prime}(x)=x^{T}Q^{\prime}x+\lambda^{\prime T}x}$. We can use the same trick as above to convert any $\lambda$ to any other $\lambda^{\prime}$, by replacing $S_{j}$ with the $Z_{j}$ gate, i.e. $\mathrm{diag}(1,-1)$ acting on the $j$th qubit. For $Q$ we can use the controlled-$Z$ gate between the $j$th and $k$th qubit, which we denote $CZ_{jk}$. This has the following effect on a basis state:

[TABLE]

where $M_{jk}$ is the $n\times n$ matrix with a $1$ in position $(j,k)$ and zeroes everywhere else. The set of all $\quantity{M_{jk}}$ form a basis for $n\times n$ binary matrices, hence we can convert any $Q$ to any other $Q^{\prime}$ by an appropriately chosen sequence of $CZ$ gates, leaving $d$ and $\lambda$ untouched. This completes the proof of statement 1.

Now to prove statement 2. From statement 1 any stabiliser state $\ket{\phi}$ is equivalent to:

[TABLE]

up to some diagonal Clifford, for some $\mathcal{K}$. The strategy is to find a Clifford unitary $U$ that commutes with $\mathcal{E}_{D}$, and converts the $2n$-qubit stabiliser state $\ket{\mathcal{K}}$ to some product of two $n$-qubit states $\ket{\mathcal{K}^{\prime}}=\ket{\mathcal{K}^{\prime}_{A}}\otimes\ket{\mathcal{K}^{\prime}_{B}}$. Then we have:

[TABLE]

where the last step follows as $\ket{\mathcal{K}^{\prime}_{B}}$ is a stabiliser state so makes no contribution to the robustness. The final state $\ket{\mathcal{K}^{\prime}}$ can be factored as $\ket{\mathcal{K}^{\prime}_{A}}\otimes\ket{\mathcal{K}^{\prime}_{B}}$ provided its generator $G^{\prime}$ can be written in block matrix form as:

[TABLE]

where $G^{\prime}_{A}$ and $G^{\prime}_{B}$ have $n$ rows, and represent the generators for affine spaces $\mathcal{K}^{\prime}_{A}$ and $\mathcal{K}^{\prime}_{B}$.

We now show that we can always reach this form by a Clifford $U_{C}$ comprised of a sequence of CNOTs targeted on the last $n$ qubits. Such a sequence always commutes with $\mathcal{E}_{D}\otimes\mathds{1}_{n}$. Suppose we have some $2n\times k$ generator $G$ for an affine space $\mathcal{K}$ with $k=\dim(\mathcal{K})$:

[TABLE]

where $G_{A}$ and $G_{B}$ are each $n\times k$ submatrices. The full matrix $G$ will have rank $k$, and $G_{A}$ will have some rank $m\leq k$. Either $G_{A}$ is already full rank ($m=k$), or it can be reduced to the following form by elementary column operations, which is equivalent to multiplication on the right by a $k\times k$ matrix $S$:

[TABLE]

where $G_{A}^{\prime}$ is $n\times m$ (and hence full column rank), and [math] is $n\times(k-m)$. Multiplying $G$ on the right by $S$, we interpret as a change in the choice of generating set:

[TABLE]

Now, apply the Clifford $U_{C}$ described by the matrix $C$ in equation (A136). This transforms the generator to:

[TABLE]

Note that if $G_{A}$ was already full rank, the change of generating set is not necessary. If we can set the bottom-left submatrix to zero, then $U_{C}\ket{\mathcal{K}}$ can be factored as described above. This is possible if there exists a binary matrix $M$ such that $MG^{\prime}_{A}=G^{\prime\prime}_{B}$. But $G^{\prime}_{A}$ has full column rank $m$, so there exists an $m\times n$ left-inverse $G^{\prime-1}_{A,\mathrm{left}}$ such that $G^{\prime-1}_{A,\mathrm{left}}G^{\prime}_{A}=\mathds{1},$ where $\mathds{1}$ is $m\times m$. Then we can set $M=G^{\prime\prime}_{B}G^{\prime-1}_{A,\mathrm{left}}$, so that:

[TABLE]

Then $G^{\prime}=CGS$ is in the form (A147), so $U_{C}\ket{\mathcal{K}}=\ket{\mathcal{K}^{\prime}_{A}}\otimes\ket{\mathcal{K}^{\prime}_{B}}$, as required. ∎

Lemma A1 shows that if $\mathcal{E}_{D}$ is diagonal then for any $2n$-qubit stabiliser state $\ket{\phi}$ we have that $\mathcal{R}\quantity(\quantity(\mathcal{E}_{D}\otimes\mathds{1}_{n})\outerproduct{\phi}{\phi})=\mathcal{R}\quantity(\mathcal{E}_{D}\quantity(\outerproduct{\mathcal{K}}{\mathcal{K}}))$ for some $n$-qubit affine space $\mathcal{K}$. This shows that the capacity can be calculated by maximising over just the representative states $\ket{\mathcal{K}}$, proving Theorem A1. Table 3 illustrates the reduction in problem size. For example, whereas naively for a two-qubit channel we would need to calculate robustness for all $36,720$ four-qubit stabiliser states, using the result above we only need check one stabiliser state for each of the $7$ non-trivial affine spaces. Cases up to five qubits are now tractable using this method.

E.4 Dimension of affine space

Here we make further observations that will help interpret numerical results from Section 8 of the main text.

Observation A1** (Dimension of affine space limits achievable robustness).**

Suppose $U$ is a diagonal unitary acting on $n$ qubits, and suppose $\ket{\mathcal{K}}$ is a stabiliser state associated with some affine space $\mathcal{K}$, $k=\dim(\mathcal{K})$. Then $\mathcal{R}(U\ket{\mathcal{K}})=\mathcal{R}(U^{\prime}\ket{\phi^{\prime}})$ where $U^{\prime}\ket{\phi^{\prime}}$ is a state on only $k$ qubits, and $U^{\prime}$ is some $k$-qubit unitary. Therefore $\mathcal{R}(U\ket{\mathcal{K}})$ is upper-bounded by the maximum robustness achievable for a $k$-qubit state.

Proof.

We prove the result by showing that there is a sequence of Clifford gates that takes $U\ket{\mathcal{K}}$ to the product of a $k$-qubit state and an $(n-k)$-qubit stabiliser state. We know from Lemma A1 that for diagonal unitaries, all states with same affine space result in the same robustness, so it is enough to consider the state:

[TABLE]

A diagonal unitary will map this to:

[TABLE]

where $\quantity{e^{i\theta_{x}}}$ will be some subset of the diagonal elements of $U$. The affine space $\mathcal{K}$ will have a generator matrix of rank $k$. As we saw in Lemma A1, a sequence of elementary row operations on the generator matrix can be realised by a sequence of CNOT gates. So we can use Clifford gates to transform any rank $k$ generator matrix as:

[TABLE]

where $\mathds{1}$ is the $k\times k$ identity. Each element of $\mathcal{K}$ can be written $x=\sum_{j}g_{j}+h$, where $\sum_{j}g_{j}$ is some combination of columns of $G$, and $h$ is a fixed shift vector. The transformation $A$ corresponds to a sequence of CNOTs that we collect in a single Clifford unitary $U_{A}$, that acts on $n$-qubit computational basis states $\ket{x}$, where $x\in\mathcal{K}$, as follows:

[TABLE]

where $h^{\prime}$ is an $(n-k)$-length vector, and $y(x)$ is a $k$-length vector given by:

[TABLE]

Note that $y(x)$ is only defined for $x\in\mathcal{K}$, and that $h^{\prime}$ is independent of $x$. Elements $x\in\mathbb{F}_{2}^{n}$ that are not in $\mathcal{K}$ could be mapped to a vector where the last $n-k$ bits are not $h^{\prime}$, but these never appear as terms of $U\ket{\mathcal{K}}$. Since $U_{A}$ must preserve orthogonality, each $\ket{x}$, where $x\in\mathcal{K}$, maps to a distinct element of the $k$-qubit basis set $\quantity{\ket{y}}$. In fact, since $y$ are length $k$ and there are $2^{k}$ distinct elements, they must form the $k$-bit linear space $\mathcal{L}^{\prime}=\mathbb{F}_{2}^{k}$. So we can write:

[TABLE]

where $\ket{\mathcal{L}^{\prime}}$ is a $k$-qubit stabiliser state, and $U^{\prime}$ is the $k$-qubit diagonal unitary with $e^{i\theta^{\prime}_{y(x)}}=e^{i\theta_{x}}$ as the non-zero elements. The state $\ket{h^{\prime}}$ is a stabiliser state, so cannot contribute to the robustness of $U_{A}U\ket{\mathcal{K}}$, and therefore $\mathcal{R}(U\ket{\mathcal{L}\quantity(\mathcal{K})})=\mathcal{R}(U_{A}U\ket{\mathcal{L}\quantity(\mathcal{K})})=\mathcal{R}(U^{\prime}\ket{\mathcal{L}^{\prime}})$, where $U^{\prime}\ket{\mathcal{L}^{\prime}}$ is a $k$-qubit state. ∎

Recall that in Section 8 of the main text, we found that highly structured examples of diagonal unitaries $U$ exist where $\mathcal{C}(U)$ is strictly larger than $\mathcal{R}(\Phi_{U})$, whereas for all the random diagonal unitaries sampled, we found them to be exactly equal. We can now explain this by a concentration effect, in conjunction with Observation A1. The $n$-qubit random diagonal gates concentrate (with high probability) within a narrow range of values for the magic capacity, close to the maximum possible magic capacity for an $n$-qubit diagonal gate. If $\mathcal{R}(\Phi_{U})<C(U)$ then by Theorem  A1 we must have that $C(U)=\mathcal{R}(U\outerproduct{\mathcal{K}}{\mathcal{K}}U^{\dagger})$ for some affine space $\mathcal{K}$ of non-maximal dimension. However, $U\outerproduct{\mathcal{K}}{\mathcal{K}}U^{\dagger}$ is Clifford equivalent to an $(n-1)$-qubit stabiliser state acted on by a diagonal unitary. Then $\mathcal{R}(U\outerproduct{\mathcal{K}}{\mathcal{K}}U^{\dagger})$ would be upper bounded by the maximum $C(\mathcal{E})$ for $(n-1)$-qubit diagonal unitaries. But if $C(\mathcal{E})$ is close to the maximum possible for $n$-qubit diagonal unitaries, then it is impossible for $U\outerproduct{\mathcal{K}}{\mathcal{K}}U^{\dagger}$ to achieve the magic capacity.

Finally, we consider the special case of multi-control phase gates $M_{t,n}$, which we defined in the main text as:

[TABLE]

Note that the gate $M_{t,n}$ acts as the identity on states $\ket{\mathcal{K}}$ unless $\mathcal{K}$ contains the zero vector $0^{n}=(0,\ldots,0)^{T}$, so if $0^{n}\notin\mathcal{K}$, we get $\mathcal{R}(M_{t,n}\ket{\mathcal{K}})=1$. But if $0^{n}\in\mathcal{K}$, then $\mathcal{K}$ is a linear subspace. So for this type of gate, to find all possible values of $\mathcal{R}(M_{t,n}\ket{\mathcal{K}})>1$ we need only consider linear subspaces. The following theorem implies that we actually only need solve one optimisation for each possible dimension of linear subspace rather than one for every linear subspace.

Theorem A2**.**

Consider the $n$-qubit gate $M_{t,n}$ defined by equation (A160), and let $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ be linear subspaces such that $\dim(\mathcal{L}_{A})=\dim(\mathcal{L}_{B})=k$. Then:

[TABLE]

Proof.

We largely repeat the arguments of Observation A1, for the special case where the phases are given by:

[TABLE]

Since $\dim(\mathcal{L}_{A})=\dim(\mathcal{L}_{B})$, their generator matrices $G_{A}$ and $G_{B}$ have the same rank. It follows from the arguments of Observation A1 that there exists an invertible $C$, corresponding to a sequence of CNOT gates, such that $G_{B}=CG_{A}$, and $\ket{\mathcal{L}_{A}}=U_{C}\ket{\mathcal{L}_{A}}$, where $U_{C}$ is a unitary Clifford operation.

If we consider instead the state $M_{t,n}\ket{\mathcal{L}_{A}}$, which involves terms in the same basis vectors as $\ket{\mathcal{L}_{A}}$, we just need to track what happens to the phase $\exp(i\theta_{0})$. Clearly, since any CNOT acts as the identity on $\ket{0^{n}}$, we obtain:

[TABLE]

Since $U_{C}$ is a reversible Clifford operation, by monotonicity of robustness of magic, equation (A161) follows. ∎

From Theorem A2, then, to find $\mathcal{C}\quantity(M_{t,n})$, we only need calculate $\mathcal{R}\quantity(M_{t,n}\ket{\mathcal{L}})$ for a single representative subspace for each possible value of $\dim(\mathcal{L})$. Recall that for $n$-qubit stabiliser states $\ket{\mathcal{L}}$, $k=\dim{\mathcal{L}}$ can take integer values from [math] to $n$. The states with $k=0$ correspond to single computational basis states without superposition, so are unaffected by phase gates. That is, for $n$-qubit multicontrol phase gates we only have to calculate $n$ robustnesses. Compare this to the number of optimisation problems we would need to solve without using the above observations (Table 3).

We can go further. From Observation A1 we know that for a subspace with $\dim(\mathcal{L})=k<n$, it must be the case that $M_{t_{n}}\ket{\mathcal{L}}$ is Clifford-equivalent to $(U^{\prime}\ket{\mathcal{L^{\prime}}})\otimes\ket{h^{\prime}}$ for the $k$-qubit state $\ket{\mathcal{L^{\prime}}}$ and $(n-k)$-qubit computational basis state $\ket{h^{\prime}}$, and some diagonal $k$-qubit unitary $U^{\prime}$. By inspection of the phases given by equation (A162), $U^{\prime}$ can only be the $k$-qubit multicontrol gate $M_{t,k}$. This leads to the following statement:

Observation A2** ($n$-qubit multicontrol gates).**

For any fixed $t$ and $n$-qubit state $\ket{\mathcal{L}}$ where $\dim(\mathcal{L})=k<n$, we have:

[TABLE]

where $\ket{\mathcal{L}^{\prime}}$ is the $k$-qubit state with $\mathcal{L}^{\prime}=\mathbb{F}_{2}^{k}$.

Observation A2 partially justifies our Conjecture 8.1 in Section 8 of the main text, that for fixed $t$, the maximum increase in robustness achievable for $M_{t,n}$, over any $n$, is given by $\mathcal{R}\quantity(M_{t,K}\ket{+}^{\otimes K})$, for some finite number of qubits $K$. To unpack this claim further, let us consider the maximisation over input stabiliser states performed to calculate the capacity $\mathcal{C}$. In this Appendix, we have seen that for the family of gates $M_{t,n}$, we only need to calculate robustness for one representative input stabiliser state for each possible dimension of linear subspace; that is, for $M_{t,n}$ there are only $n$ robustnesses to calculate. In Table 4 we present the relevant values for the family of multicontrol-$T$ gates ($t=2$) and make two observations. First, looking across the rows of Table 4, notice that the values for fixed $k$ are constant with $n$, assuming $k\leq n$. Indeed, this is a generic feature of the $M_{t,n}$ gates as formalised by Observation A2. Second, looking down the last column of Table 4, we see that up until $k=4$, $\mathcal{R}\quantity(M_{t,n}\ket{\mathcal{L}})$ increases with $\dim(\mathcal{L})$, but at $k=5$ the value drops. With a little thought we can see that this is necessarily the case if $\mathcal{R}\quantity(\Phi_{M_{t,5}})<\mathcal{C}\quantity(M_{t,5})$; we saw earlier that for diagonal gates $U$ the Choi state robustness is equal to $\mathcal{R}\quantity(U\ket{+}^{\otimes n})$, and $\ket{+}^{\otimes n}$ is a representative state for the $k=n$ case.

Our current techniques limit us to five-qubit operations, so we are unable to confirm whether $\mathcal{R}(M_{t,n}\ket{\mathcal{L}})$ continues to decrease with increasing $\dim(\mathcal{L})$. An intuition for why a decrease is plausible goes as follows. A stabiliser state $\ket{\mathcal{L}}$ with $\dim(\mathcal{L})=k$ will have $2^{k}$ equally weighted terms when written in the computational basis, so will have a normalisation factor of $2^{-k/2}$. The non-stabiliser state $M_{t,n}\ket{\mathcal{L}}$ is identical to $\ket{\mathcal{L}}$ apart from the phase on the all-zero term $\ket{0\ldots 0}$. As $k$ becomes large, the amplitude of the term $\frac{e^{i\pi/2^{t}}}{2^{k/2}}\ket{0\ldots 0}$ becomes very small, so that $M_{t,n}\ket{\mathcal{L}}$ has high fidelity with the stabiliser state $\ket{\mathcal{L}}$. We would therefore expect $M_{t,n}\ket{\mathcal{L}}$ to have a small robustness if $k$ is large.