Scalable randomized benchmarking of non-Clifford gates

Randomized benchmarking is a widely used experimental technique to characterize the average error of quantum operations. Benchmarking procedures that scale to enable characterization of $n$-qubit circuits rely on efficient procedures for manipulating those circuits and, as such, have been limited to subgroups of the Clifford group. However, universal quantum computers require additional, non-Clifford gates to approximate arbitrary unitary transformations. We define a scalable randomized benchmarking procedure over $n$-qubit unitary matrices that correspond to protected non-Clifford gates for a class of stabilizer codes. We present efficient methods for representing and composing group elements, sampling them uniformly, and synthesizing corresponding $\mathrm{poly}(n)$-sized circuits. The procedure provides experimental access to two independent parameters that together characterize the average gate fidelity of a group element.

Randomized benchmarking is a widely used experimental technique to characterize the average error of quantum operations. Benchmarking procedures that scale to enable characterization of n-qubit circuits rely on efficient procedures for manipulating those circuits and, as such, have been limited to subgroups of the Clifford group. However, universal quantum computers require additional, non-Clifford gates to approximate arbitrary unitary transformations. We define a scalable randomized benchmarking procedure over n-qubit unitary matrices that correspond to protected non-Clifford gates for a class of stabilizer codes. We present efficient methods for representing and composing group elements, sampling them uniformly, and synthesizing corresponding poly(n)-sized circuits. The procedure provides experimental access to two independent parameters that together characterize the average gate fidelity of a group element. A key step to realizing a large-scale universal quantum computer is demonstrating that decoherence and other realistic imperfections are small enough to be overcome by fault-tolerant quantum computing protocols [1,2]. Randomized benchmarking (RB) [3][4][5][6] has become a standard experimental technique for characterizing the average error of quantum gates due in part to its insensitivity to state preparation and measurement errors. Benchmarking provides robust estimates of average gate fidelity [6,7] and can characterize specific interleaved gate errors [8,9], addressability errors [10], and leakage errors [11][12][13].
RB techniques that efficiently scale to many qubits have been limited to subgroups of gates in the Clifford group, since computations with this group are tractable [6]. However, the Clifford group is not enough for general quantum computations [14]. Previous work generalizes RB to groups that include non-Clifford gates [15,16], but only on single qubits, a significant limitation.
In this Letter we present a scalable RB procedure that includes important non-Clifford circuits, such as circuits composed from T = 4 √ Z and controlled-NOT (CNOT) gates that naturally occur in fault-tolerant quantum computations. The n-qubit matrix groups we study are a generalization of the standard dihedral group and coincide in some cases with protected gates in stabilizer codes, such as k-dimensional color codes [17]. Circuits built from these gates cannot be universal but do constitute significant portions of magic state distillation protocols [18,19], repeat-until-success circuits [20], and the vital quantum Fourier transform [21]. We show that there are efficient methods for representing and composing group elements, sampling them uniformly, and synthesizing corresponding circuits whose size grows polynomially with the number of qubits n. The benchmarking procedure provides experimental access to two independent noise parameters through exponential decays of average sequence fidelities.
The quantum circuits we consider are products of CNOT gates Λ 12 (X)|u, v := |u, u ⊕ v , bit-flip gates X|u := |u⊕1 , and single qubit m-phase gates Z m |u := ω u m |u where ω m = e i2π/m . More concisely, the circuits of interest are given by the group We call this group a CNOT-dihedral group since it is generated by CNOTs and a single qubit dihedral group [28]. Although we prove certain results for general m, we focus mainly on the case of m = 2 k as this affords efficient benchmarking and contains various non-Clifford gates of interest, such as T , controlled-S, and controlledcontrolled-Z [29]. RB over G 2 k -The benchmarking procedure we present here both generalizes [16] [30] and extends naturally to interleaving gates to estimate individual gate fidelities [8,9]. The procedure is as follows. Choose a sequence of + 1 unitary gates where the first gates are uniformly random elements g j1 , g j2 , . . . , g j of G 2 k and the ( + 1) th gate is g −1 j := g † j1 . . . g † j where j denotes the -tuple (j 1 , . . . , j ) labeling the sequence. We show later that elements of G 2 k can be efficiently sampled and g −1 j can be efficiently computed. For each sequence, we prepare an input state ρ, apply S j := g −1 j g j . . . g j1 , and measure an operator E. The overlap with E is Tr[ES j (ρ)]. Averaging this overlap over K independent sequences of length gives an estimate of the average sequence fidelity F seq ( , E, ρ) := Tr[ES (ρ)] where S (ρ) := 1 K j S j (ρ) is the average quantum channel. The error operator of the final gate g j is attributed to measurement error, perturbing E to a new operator E . We decompose the input state and this final measurement operator in the Pauli basis to give ρ = P x P P/2 n and E = P e P P . The average sequence fidelity is [10] where A Z = P ∈Z\{I} e P x P and A R = P ∈P\Z e P x P . In a spirit similar to simultaneous RB [10], by choosing arXiv:1510.02720v1 [quant-ph] 9 Oct 2015 appropriate input states, namely |0 . . . 0 and |+· · · + := b∈{0,1} n |b , each of the two exponential decays can be observed independently. State preparation errors may lead to deviation from single exponential decay, but this is detectable. The channel parameters α Z and α R can be extracted by fitting the average sequence fidelity to such a decay. The corresponding depolarizing channel parameter is a weighted average α = (α Z + 2 n α R )/(2 n + 1), and the average gate error is given by r = (2 n −1)(1− α)/2 n (see [6]).
Several natural questions arise from this work. First, one might address the asymptotically optimal cost of circuit synthesis for elements of the CNOT-dihedral groups, as well as the practical question of finding optimal circuit decompositions for elements of the smallest groups. We expect optimal circuits are computationally hard to find as n grows, but experimentally it is important to minimize the number of two-qubit gates. Second, unlike the Clifford group, the CNOT-dihedral group is not a 2-design [5]. It would be interesting to find a group (or set) containing a non-Clifford gate and that is a 2design and in which benchmarking can be done efficiently. Third, our results show that we can efficiently perform RB. However, we have not addressed the precise sense in which quantum computations over the CNOT-dihedral group can be efficiently simulated. This may be a subtle problem [22,23]. Lastly, there are generalized stabilizer formalisms, such as [24], and it is natural to wonder if one of these describes how this group acts on some set of states.
The remainder of the Letter is devoted to proving the various results utilized in the benchmarking procedure: canonical decomposition of G m , efficient computation in G m , and twirling over G m , each of which is interesting in its own right. Let m be general and let us briefly set some notation. The matrix representation of G m is set by identifying g ∈ G m to the matrix that maps |0 n := |00 . . . 0 to |b := |b 1 b 2 . . . b n with unit phase. We define the phase-flip gates Z|u := (−1) u |u and controlled-Z (CZ) gates Λ 12 (Z)|u, v := (−1) uv |u, v . The support of a bit string v ∈ {0, 1} n is supp(v) = {j|v j = 1} ⊆ [n] := {1, 2, . . . , n}. We refer to v and its support interchangeably, treating v as a set and vice versa. Let U be a single qubit gate and U (v) denote the gate acting as U only on qubits in the support of v. Given J ⊆ [n] or elements i, j, · · · ∈ [n], we also use the shorthand U (J) and U (i, j, . . . ). P := X(j), Z(j) / i denotes the nqubit Pauli group, and we define X := X(j)|j ∈ [n] , Z := Z(j)|j ∈ [n] , cX := Λ ij (X)|i, j ∈ [n], i = j , and cZ := Λ ij (Z)|i, j ∈ [n], i < j .
Canonical form of G m -Our first goal will be to put G m in a canonical form (the main result is contained in Theorem 1). The rewriting identities shown in Fig. 1 allow us to commute diagonal elements of G m through Λ ij (X) and X(j) gates. The rules for bit-flip gates are a special case of the CNOT rules. The following Lemma follows directly from definitions and formalizes the role of the rewriting identities in understanding the group's structure.
Lemma 1. Let W m denote the subgroup of diagonal matrices of G m and let Π = Λ ij (X), X(j) denote the subgroup of permutation matrices. Then G m is isomorphic to a semi-direct product of groups G m W m Π.
The proof of Lemma 1 is given in the Supplemental Material [25]. Note that by definition, Π = X cX . Since X F n 2 and cX GL n (F 2 ), each element of Π can be associated to an n-bit string c ∈ F n 2 and an n by n invertible 0 − 1 linear transformation B ∈ GL n (F 2 ) such that π|b = |Bb ⊕ c . Here, F 2 denotes the field with two elements. Furthermore, |Π| = 2 n n−1 =0 (2 n − 2 ). It remains to better understand W m (see Lemma 3 for the main result). Let D m denote the group of 2 n by 2 n diagonal unitary matrices D with elements b|D|b = ω Here f : F n 2 → Z m is a function that assigns mth roots of unity to the diagonal and Z m is the ring of integers modulo m. Since G m is generated by permutation matrices and products m-phase gates, W m ⊆ D m .
Let R ⊂ Z m [x 1 , . . . , x n ] denote the polynomial ring whose elements are p(x) := p(x 1 , . . . , is a monomial. The multi-index takes values in {0, 1} n as a convenient notation since we will evaluate p(x) on binary strings, so The degree of a monomial is denoted |α|. We mainly consider R as an additive group. The next Lemma follows from the definition of group isomorphism and the fact that each function f (b) can be expressed as a polynomial in R.
The proof of Lemma 2 is given in the Supplemental Material [25]. The rewriting identities give the action of Π on W m by conjugation. LetW m := Z m (j) . Based on a similar application of the rewriting identities as in Lemma 1, W m = πW m π † | π ∈ Π . Since W m ⊆ D m R, Φ −1 associates a polynomial in R to each element of W m . By our chosen convention, matrices representing elements w ∈ W m are given modulo a global phase factor ω m such that w|0 n = |0 n . Therefore the preimages Φ −1 (w) have zero constant term, i.e. p α = 0 when |α| = 0. Through Φ, the rewriting identities define an action of Π on R that respectively takes Eq. 3 increments the degree of a monomial and multiplies its coefficient by −2, whereas Eq. 4 does not change the degree. Another way to understand iterated applications of Eq. 3 is to observe that This fact relates how single qubit Z m gates acting on mod 2 linear combinations of input bits are equivalent to products of certain controlled-phase gates.
There is an element of W m corresponding to each monomial term of non-zero degree, and the coefficient of this term has the form p α ∈ (−2) |α|−1 Z m , as we will now see [25]. We choose a subset of qubits J, fix any j ∈ J, and define a permutation gate and corresponding polynomial By Eq. 5, Φ(p J ) = π J,j Z m (j)π † J,j ∈ W m , i.e. this circuit has a polynomial with one term of degree |J|. Since Φ(Z m (j)) = x j , scaled monomials of successive degrees |α| and with coefficients in (−2) |α|−1 Z m can be generated inductively by composing these circuits. Take all linear combinations of these over Z m to find We can now directly compute |G m |.
Proof. Let o m (a) = LCM(a, m)/a denote the order of a in Z m . Observe that (−2) |α|−1 Z m Z om(2 |α|−1 ) as additive groups. Therefore W m is isomorphic to a direct product of additive cyclic groups A m := n t=1 Z ( n t ) om(2 t−1 ) . This shows that |G m | = |A m ||Π|.
Putting everything together, we have Theorem 1. Any element of G m can be written in canonical form as the composition of a sequence of phase gates (comprising an element of W m whose form is given in Lemma 3), a sequence of CNOT gates, and a sequence of bit-flip gates.
Efficient computation in G 2 k -Our next goal is to present efficient methods for computing with G m . Suppose we consider a fixed value of m. Any labeling of group elements will have length proportional to s = log 2 |G m |.
Therefore, s = Ω(2 n ) whenever m is odd, and in general we cannot efficiently represent elements of G m as the number of qubits grows. However, s = O(n k ) for the special case m = 2 k , and the story is different. We focus on this special case for the remainder of this Letter.
An element g ∈ G m can be written as a product g = uvw where w ∈ W m is a diagonal matrix, v ∈ cX is a CNOT circuit, and u ∈ X is a tensor product of bitflips. This transforms n-qubit quantum states as g|b = ω p(b) m |Bb + c where p ∈ W, B ∈ GL n (F 2 ), and c ∈ F n 2 . Group elements are in bijective correspondence with the triples (p, B, c). The polynomial p has maximum degree k and at most k t=0 n t = O(n k ) nonzero coefficients, each contained in Z 2 k .
The product of group elements g 1 , g 2 ∈ G m , is given by the triple The products B 2 B 1 and B 2 c 1 ⊕ c 2 can be computed in O(n 3 ) time, and polynomials in W can be added in O(n k ) ring operations. We need to show that p 2 (B 1 x ⊕ c 1 ) can also be computed efficiently. Consider a triple (p, B, c) and let B j denote the jth row of B and J j = supp(B j ). Define x = Bx ⊕ c. Then for any j ∈ [n], using Eqs. 5 and 6, has maximum degree k. When we substitute x = x 1 . . . x n into the degree k polynomial p(x), computations occur with coefficients in Z 2 k . We compute each monomial (x ) α with O(k) multivariate polynomial multiplications, each of which can be done term-by-term in O(n 2k+1 ) ring operations. We compute the term (−2) |α|−1 p α (x ) α with an additional O(n k ) ring operations to multiply each term of (x ) α by a (−2) |α|−1 p α and accumulate the result. There are O(n k ) terms in p(x), so the total number of ring operations to compute p(x ) is O(n 3k+1 ). If c = 0 n , it is possible that p(x ) has a non-zero constant term. With an additional O(n k ) ring operations, p(x ) can be mapped to an equivalent polynomial in W.
Uniformly sampling from G 2 k is equivalent to uniformly and independently sampling from W, GL n (F 2 ), and F n 2 . This can be done efficiently since elements of W have maximum degree k; see also [25,26].
Given a triple (p, B, c), we synthesize a corresponding circuit from products of CNOT gates, bit-flip gates, and a single qubit m-phase gates. Our goal is to efficiently synthesize a circuit whose size (number of gates) is polynomial in n but not to optimize this circuit. We independently synthesize circuits coinciding with p, B, and c. Since c corresponds to X(c), and a CNOT circuit for B can be found by Gaussian elimination [14], the new part of the algorithm synthesizes a circuit for p.
We describe the circuit synthesis for p informally. The algorithm proceeds in k rounds. Begin by initializing a working polynomial q(x) ← p(x), set a round counter t ← k, and set a quantum circuit U ← I. Here "←" denotes assignment. In round t, we synthesize a circuit corresponding to a polynomial p (t) (x) that coincides with q(x) on its degree t terms. For each of the O(n t ) degree-t terms (−2) |α|−1 p α x α of q(x), we apply the constant-sized circuit g α := π J,j (Z 2 k (j)) pα π † J,j setting U ← g α U , where J = supp(α) as in the proof of Lemma 3. The product of the g α corresponds to p (t) (x) := α⊆[n],|α|=t p α p J (x). Therefore we update q(x) ← q(x) − p (t) (x), which now has maximum degree t − 1, decrement the round counter, and proceed to the next round. The algorithm terminates when q(x) = 0 and t = 0. The total algorithm run-time and circuit size of the output U is O(n k ).
Twirling over G 2 k -A quantum channel is a completely-positive trace-preserving map whose operator sum decomposition is E(ρ) = k A k ρA † k where k A † k A k = I. The twirl of E over a finite group G (G-twirl) is given bȳ In what follows, we use several facts about group twirls. If G = AB is a direct product of groups then where the twirl over the factor group G/A is over a set of coset representatives. Twirling any map over the Pauli group produces a Pauli channel [5]. Consider a Pauli channel E(ρ) = Q∈P η Q QρQ. Twirl this channel over any finite group G that has a permutation action on the set P. The orbit of P ∈ P is O P := {V † P V | V ∈ G} and the stabilizer is S P := {V ∈ G | V † P V = P }. The orbits define an equivalence relation P ∼ Q if and only if O P = O Q . This relation partitions P into a disjoint union of orbits. By the orbit-stabilizer theorem and Lagrange's theorem [27], |O P | = |G/S P | = |G|/|S P |. Therefore the twirl, Eq. 12, can be written where C is a set of representative elements, one from each orbit. These facts allow us to compute the twirl over G 2 k when k > 1 by expressing it as a sequence of twirls. We begin by decomposing the group. LetW 2 k := W 2 k \ (W 2 k \{I}) and recall thatW 2 k := Z 2 k (j) . Then W 2 k = W 2 kW 2 k . Since cZ W 2 k and Z W 2 k , we form the corresponding factor groups. Therefore an element w ∈ W 2 k can be written as w =ww =w 1w2w1w2 wherẽ w 1 labels cosetsw 1 cZ,w 2 ∈ cZ,w 1 labels cosetsw 1 Z, andw 2 ∈ Z. Finally, by Lemma 1, any element g ∈ G 2 k factors as g = uvw where u ∈ X , v ∈ cX , and w ∈ W 2 k . Therefore, we have g = uw 2 vw 2w1w1 wherē Our strategy is to use the decomposition to express the G 2 k -twirl as a sequential P-twirl, cX -twirl, cZ-twirl, W 2 k /Z-twirl, andW 2 k /cZ-twirl. Each twirl can be computed in a straightforward manner using the facts we have described and reduces the number of independent parameters describing the channel until we have twirled over the whole of G 2 k [25]. The final twirled map is In the Liouville representation in the Pauli basis which has matrix elements R (E) P Q = Tr(P E(Q))/4 n where P and Q are n-qubit Pauli operators, this map has three diagonal blocks corresponding to I, Z \ {I}, and P \ Z with elements 1, α Z := 1 − 4 n β R , and α R := 1 − 2 n β Z − (4 n − 2 n )β R , respectively.
Conclusion -Our results enable scalable benchmarking of a natural family of non-Clifford circuits related to quantum error-correcting codes. Because our procedure scales, in principle it allows efficient benchmarking of isolated non-Clifford gates as well as large sub-circuits for state distillation [18,19] or repeat-until-success protocols [20]. These sub-circuits can be characterized with our procedure using physical gates or logical gates on protected qubits. Together with standard Clifford benchmarking, our procedures enable characterization of the full range of gates used in the leading fault-tolerant quantum computing protocols. As multi-qubit benchmarking is well within experimental reach, we expect an optimized implementation of our procedure to be quite practical.
Supplemental Materials: Scalable randomized benchmarking of non-Clifford gates Let Λ x1x2...xpj (U ) denote a controlled-U gate targeting qubit j with controls x 1 , x 2 , . . . , x p . Let O m and E m denote odd and even strings of length m respectively. Given a subset J ⊆ [n],J := [n] \ J. For a support J and string u ∈ {0, 1} |J| , u[J] ∈ Z n 2 denotes the string whose restriction to J is u and is zero otherwise. The indicator function δ(u, v) is one if u = v and zero otherwise.
The circuit diagrams in Fig. 1 are an intuitive way to present the rewriting identities, but they can also be written concisely as For completeness, we give the proofs of Lemma 1, Lemma 2, Equation 5, and Lemma 3.
Proof of Lemma 1. By definition of the semi-direct product, we must show that (a) W m ∩ Π = {I}, where I denotes the identity matrix, (b) G m = ΠW m , and (c) W m G m . Every element of Π is a permutation matrix. Therefore the only diagonal matrix in Π is the identity matrix, proving (a). Consider an element g = g N . . . g 2 g 1 ∈ G m expressed as a product of generators. Suppose the th generator is g = Z m (j) ∈ W m . By repeated application of rewriting identities, we can write g = g m . . . g +1 g −1 . . . g 2 g 1 w where w = g † 1 . . . g † −1 g g −1 . . . g 1 ∈ W m . Repeating this argument for all occurrences of Z m (j) in g, we obtain g = πw where π ∈ Π and w ∈ W m ; thus G m < ΠW m . Conversely, since Π and W m are both subgroups of G m , any πw ∈ ΠW m < G m . This proves (b). By (b), each element of G m can be written as g = πw g for some π ∈ Π and w g ∈ W m , so gwg † = πw g ww † g π † = πw π † where w = w g ww † g ∈ W m . Decomposing π into a product of generators X(j) and Λ ij (X) gates and applying rewriting identities to each element in the product, we find πw π † = w ∈ W m and thus prove (c).
The corollary to Lemma 3 allows us to compute the order of G m for any m and any number of qubits n. Both to verify this formula and get a sense of the size of the group, we have tabulated orders of the smallest groups G m in Table I by explicitly constructing them. The group order increases rapidly with n but is comparable in size to the Clifford group.  We provide more details about how to sample a uniformly random element of G m . Assume we have access to a string of uniformly random bits. We take n bits as the element c ∈ Z n 2 . In time M (n) + O(n 2 ) and with n 2 + 3 random bits, where M (n) is the time to multiply two n by n matrices, we can generate a uniformly random non-singular matrix B ∈ GL n (F 2 ) [S1]. Finally, for all α satisfying 0 < |α| ≤ k, we take k random bits to draw a uniformly random element from Z 2 k , which we multiply by (−2) |α|−1 modulo 2 k to obtain the corresponding coefficient of the polynomial p ∈ W. There are O(n k ) non-zero coefficients, so we consume O(kn k ) random bits to sample an element uniformly at random from W.
We use several facts to compute the G m -twirl. We prove these facts now. If G = AB is a direct product of groups A and B, then where the sum over T ∈ G/A is over a set of coset representatives, i.e. each T is any representative in the coset T A.
It is well-known that twirling any map over the Pauli group produces a Pauli channel. This can be seen from We have expressed each Kraus operator in the Pauli operator basis A k = Q γ (k) Q Q. We wrote P QP = (−1) ω(P,Q) Q where ω(P, Q) = 0 if [P, Q] = 0 and ω(P, Q) = 1 otherwise, and used the fact that 1 |P| P ∈P (−1) ω(P,Q)+ω(P,R) = δ(P, Q). Trace preservation ensures that Q k |γ (k) Q | 2 = 1. This twirl reduces a general map to one described by 4 n parameters. For convenience let η Q := k |γ (k) Eq. 13 can be derived from the definition of the twirl as follows For completeness, we conclude now with details about each step of the G 2 k -twirl: the cX -twirl, the cZ-twirl, the explicit orbits for each, theW 2 k /Z-twirl, and the finalW 2 k /cZ-twirl.
The cZ-twirl has an exponential number of orbits. Let u * ∈ [n] denote the first non-zero coordinate of a non-zero u ∈ F n 2 . The orbits correspond to representative elements C cZ := X ∪ (Z \ {I}) ∪ {X(u)Z(u * ) : u ∈ F n 2 , u = 0}. There are 2 n − 1 orbits of each non-identity type and |O X(u) | = |O X(u)Z(u * ) | = 2 n−1 for all u. Comparing the cZ orbits to the cX orbits, we see that Furthermore O X(u) contains one element from O X(1) and 2 n−1 − 1 elements from O Y (1,2) . Therefore parameters for O X(1) and O Y (1,2) are averaged into a new parameter β X−Y Y = (β X + (2 n−1 − 1)β Y Y )/2 n−1 that is simply the average over all parameters η Q with Q ∈ O X(1) ∪ O Y (1,2) upon substituting β X and β Y Y . After the cZ-twirl, the Pauli channel is described by 4 parameters.
These are the explicit orbits of the action of cX and cZ on P. The cX orbits are where J ⊆ [n] ranges over all subsets and v ∈ F |J| 2 ranges over all strings. It is straightforward to verify that these are the orbits since Λ(X) preserves the parity of the number of Y -type tensor factors and does not mix X-type and Z-type factors. The cZ orbits are where u ∈ F n 2 , u = 0, v ∈ F |J| 2 and J = supp(u). This is also straightforward to show from the action of Λ(Z) on P. TheW 2 k /Z-twirl is over coset representatives of Z 2 k (j) /Z. Let P = X(u)Z(v) ∈ P and P (w) = X(J)Z(w[J]+v) where J = supp(u). Twirling all tensor factors gives P ρP → 1 2 |J| w∈F |J| 2 P (w)ρP (w) † . (S13) Tensor factors that are X or Y become equal sums of both, whereas Z commutes. Therefore the cZ-twirl classes O X(u) and O X(u)Z(u * ) for fixed u (and fixed J = supp(u)) are averaged together equally. Since these orbits all have the same size and partition O X(1) ∪ O Y (1,2) and O Y (1) respectively, the parameters β X−Y Y and β Y are averaged into a new parameter β R := (β X−Y Y + β Y )/2. Again this is simply the average of all parameters η Q with Q ∈ O X(1) ∪ O Y (1,2) ∪ O Y (1) . After theW 2 k /Z-twirl, the Pauli channel is described by 3 parameters. The detailed calculations for theW 2 k /Z-twirl are as follows. Conjugating a single qubit Pauli operator X u Z v by a coset representative (Z 2 k ) a , such as occurs in a term of Eq. 12, gives cos(2πa/2 k )I + i sin(2πa/2 k )Z u X u Z v . (S14) In the (u, v) = (1, 0) case, when we twirl only the first tensor factor, P ρP is taken to 1 2 k−1 2 k−1 −1 a=0 P 2:n [ cos 2 (2πa/2 k )XρX − cos(2πa/2 k ) sin(2πa/2 k )Y ρX (S15) − cos(2πa/2 k ) sin(2πa/2 k )XρY + sin 2 (2πa/2 k )Y ρY P 2:n (S16) where P 2:n denotes remaining tensor factors of P on qubits {2, . . . , n}. Applying the identity cos 2 u+cos 2 (π/2+u) = 1, we can write 2 k−1 −1 a=0 1 2 k−1 cos 2 (2πa/2 k ) = 2 k−2 −1 a=0 1 2 k−1 cos 2 (2πa/2 k ) + cos 2 (π/2 + 2πa/2 k ) = 1/2. (S17) Similar arguments allow us to sum the remaining terms. The odd functions vanish and the even functions average to 1/2. The final steps of the twirl over G m are straightforward but we provide more details. The diagonal PTM representation of E(ρ) is R (E) QQ = β I + β Z P ∈Z\{I} (−1) ω(P,Q) + β R P ∈P\Z (−1) ω(P,Q) . (S18) Any non-identity Pauli operator commutes with exactly half of the Pauli group and anti-commutes with the other half. If Q ∈ Z \ {I} then R E QQ = β I + (2 n − 1)β Z − 2 n β R . Likewise, R (E) II = β I + (2 n − 1)β Z + (4 n − 2 n )β R = 1. When Q ∈ P \ Z, R (E) QQ = β I − β Z . Therefore the map is diagonal and has three blocks proportional to the identity. Consider a unitary operator U . If U commutes with Q or Q , then by the cyclic property of the trace R (U ) QQ = δ(Q, Q ). If U is a diagonal gate then R (U ) is block diagonal with identity matrices in the Z blocks, all off diagonal blocks equal to zero, and an arbitrary P \ Z block. Since R (U † ) = (R (U ) ) −1 and every element ofW 2 k /cZ is diagonal, the final twirled map is 1 |G| (S19)