Subspace benchmarking high-fidelity entangling operations with trapped ions

We present a new and simplified two-qubit randomized benchmarking procedure that operates only in the symmetric subspace of a pair of qubits and is well suited for benchmarking trapped-ion systems. By performing benchmarking only in the symmetric subspace, we drastically reduce the experimental complexity, number of gates required, and run time. The protocol is demonstrated on trapped ions using collective single-qubit rotations and the Molmer-Sorenson (MS) interaction to estimate an entangling gate error of $2(1) \times 10^{-3}$. We analyze the expected errors in the MS gate and find that population remains mostly in the symmetric subspace. The errors that mix symmetric and anti-symmetric subspaces appear as leakage and we characterize them by combining our protocol with recently proposed leakage benchmarking. Generalizations and limitations of the protocol are also discussed.


I. INTRODUCTION
Like any complex machine, the construction of a largescale quantum computer will require a rapid design-andtest cycle of individual components and their assemblage. Early recognition of this need by quantum information scientists led to widespread, (but not universal [1]), adoption of the process fidelity as the gold-standard for the quality of quantum operations [2]. While the quantity had solid theoretical motivations, experimentalists lacked a simple method for its measurement for two or more qubits. In practice, experimentalists were stuck choosing between comprehensive but resource-intensive process tomography [3] and convenient but less descriptive Bell-state tomography [4][5][6][7].
The advent of randomized benchmarking (RB) marked a powerful advancement in experimental validation of quantum operations [8,9]. The new tool provided algorithmic advantages, including robustness to statepreparation and measurement (SPAM) errors [10,11], multiplicative precision in estimates [12], and, under certain assumptions, an accurate direct estimate of the fidelity [13]. In the case of multi-qubit gate benchmarking, one disadvantage of the protocol is the requirement for individual addressing. While this demand is also needed for universal quantum computation, it can impose a significant engineering requirement on the design-and-test cycle, especially in trapped-ion systems. In this work, we alleviate this requirement and provide a benchmarking protocol that uses only global operations.
Other methods for RB with limited gate sets have also been considered in the context of logical RB [14,15], and as a way to simplify standard RB [16,17], but these methods suffer from added sampling overhead and require a detailed study of the gates' group structure to identify different decay rates. Our method circumvents these drawbacks.
The basic idea behind our procedure is to isolate a subspace of the total system to perform RB, hence we call * charles.baldwin@honeywell.com the method subspace randomized benchmarking (SRB). We then use leakage benchmarking [18][19][20] to identify errors that may drive population outside of this subspace. We show that these different methods can be combined into a single procedure and each effect can be individually measured. This analysis can identify many common errors, and we consider several examples numerically.
In addition to developing the theory behind SRB, we demonstrate the procedure in a trapped ion system. We benchmark the Mølmer-Sørenson (MS) interaction and collective single-qubit gates (identical single-qubit gates applied to all ions) using two ions in a four-ion crystal. The other ions are a different species and are used for sympathetic cooling to prevent the temperature from changing throughout the gate sequences. With this system, we are able to run gate sequences with as many as 400 MS gates. The SRB procedure returns an estimated entangling gate error of 2(1) × 10 −3 .
The remainder of this paper is organized into the following sections: II. Standard RB -A review of the standard randomized benchmarking protocol and how to interpret the experimental results.
III. Subspace RB -Our construction of SRB and how subspace leakage is considered.
IV. Extracting information from SRB -Analysis of how different errors are manifested in SRB.
V. Experimental platform -A description of the SRB demonstration.
VI. Results and discussion -Estimates of the subspace and leakage errors.
VII. Generalized subspace RB -Extensions to other systems.

II. STANDARD RANDOMIZED BENCHMARKING
We begin by outlining the standard RB procedure [11] with further details given in App. B. RB provides an estimate of the fidelity of a finite set of gates C = {C 1 , . . . , C N }. The standard RB construction requires that the gates form a representation of a finite group, which is usually the Clifford group. A Clifford gate is the unitary implementation of a Clifford group element and defined by their action on Pauli gates. A Pauli gate is the unitary implementation of a Pauli group element, which we label 1, X, Y , and Z and tensor products thereof. A Clifford gate maps Pauli gates to Pauli gates with a sign, e.g. CXC † = −Y . The Clifford gates are selected in RB for two reasons: (1) it is efficient to simulate their evolution by the Gottesman-Knill theorem [21], and (2) they form a unitary two-design which simplifies the analysis of an RB experiment [11].

with different
We typically pick exponentially distributed between 1 and max ≈ 1 where is the expected error [22]. We typically pick K ≈ 100 sequences, but rigorous estimates also exist based on expected fidelity [23,24].
Although the net action is chosen to be the identity operation, noise and errors in the state-preparation, gates, and measurement cause p( ) ≤ 1. The power of RB is to relate the average survival p( ) to the process fidelity of the Clifford gates. If each Clifford gate has the same error channel then randomization over the Clifford gates depolarizes each error channel due to the unitary two-design property [10]. The resulting average survival is where A and B relate to SPAM errors, and r is referred to as the standard decay rate. By evaluating p( ) for different values of , we can extract A, B, and r by leastsquares fitting. The process fidelity of the Clifford gate's error channel is then defined as In practice, Pauli and Clifford gates are implemented by compiling physically available gates. In trapped-ion Circuits used for standard RB and SRB. Blue gates are collective single-qubit gates, gray gates are entangling phase gates, and green gates are individually addressed singlequbit gates. (a) Standard RB for a transport-based trappedion architecture. Individually addressed gates require a transport operation and additional cooling, which involves significant engineering overhead and may increase the error rate per Clifford. Alternative individual addressing methods exist but also have engineering overheads and similarly affect the error rates. (b) SRB only has collective and phase gates that do not require transport and cooling can be implemented during the collective gates. This reduces total run time and other errors.
systems (and many other quantum systems) the physical single-qubit gates are implemented as rotations in the X −Y plane of the Bloch sphere with arbitrary amplitude . The entangling operations for many trapped-ion experiments, including ours, is the MS gate, As discussed later, we use the phase gate U ZZ = exp[−i π 4 ZZ], which is created by applying singlequbit "wrapper pulses" before and after the MS gate [25].

III. SUBSPACE RANDOMIZED BENCHMARKING
Our procedure for RB, which we call subspace randomized benchmarking (SRB), follows the same general prescription outlined above except we consider a restricted set of physical gates. This restriction could be due to the physical limitations, i.e. the system cannot implement a universal set, or because some gates have much lower fidelity and must be avoided. In either case, the remaining set of gates is not universal, rendering standard RB impossible. We now construct a protocol with a restricted set of gates that retains most of RB's favorable attributes and still yields detailed information about the gate in question.
We introduce the procedure in the context of a trapped-ion system where the two qubits are encoded in two ions in the same potential well. Without tightly focused laser beams or complex transport operations, the available one-and two-qubit gates will be symmetric with respect to qubit exchange. This limited set of physical gates cannot generate the two-qubit Clifford gates needed for standard RB. Nevertheless, the goal of the setup is to determine a quantitative measure of the entangling gate's performance. This is a common setup for many trappedion experiments that are trying to demonstrate new entangling gates or high-fidelity performance [4,5,7].
The first step in SRB is to identify a subspace of the total Hilbert space where we can perform RB, which we call the RB subspace. For the trapped-ion system, the symmetric subspace of the two-qubits is a natural choice. States in this subspace are symmetric under the SWAP operation, with basis { 00⟩ , 1 √ 2 ( 01⟩ + 10⟩), 11⟩}. The remaining space is spanned by the anti-symmetric state 1 √ 2 ( 01⟩ − 10⟩), and we refer to it as the leakage subspace. It is easy to prepare the 00⟩ state via optical pumping with high fidelity [26]. The collective singlequbit gates and phase gate are block-diagonal in this basis, and therefore maintain population in each subspace. Measurements in the system return the number of qubits in the 1⟩ state (the number of ions fluorescing). If the population is isolated to the symmetric subspace then this measurement corresponds to a projection measurement in the basis for RB subspace defined above. The next step is to identify a set of gates in the RB subspace that form a unitary two-design, which was the property that allowed us to derive Eq. (1). In standard RB, we use Clifford gates, which map Pauli gates to Pauli gates. For prime dimensions, we can generalize the Pauli gates to Weyl gates and the corresponding Clifford gates map Weyl gates to Weyl gates [27]. These generalized Cliffords also form a unitary-two design, which will act to depolarize any error channel [28]. For d 1 = 3 (or qutrit) there are 9 qutrit-Weyl gates and 216 qutrit-Clifford gates.
To generate each qutrit-Clifford gate we use techniques from optimal control theory. In the RB subspace the collective gates act as spin-1 J x and J y angular momen-tum operators, and the phase gate acts as a J 2 z operator. Then, control of the symmetric subspace is analogous to control of a spin-1 particle. We make use of the Cartan decomposition [29,30] to decide the order of collective and phase gates and then perform a numerical search [31] for each of the 216 qutrit-Clifford gates to find amplitudes and phases. We find that only three phase gates are required per qutrit-Clifford arranged as shown in Fig. 1. This does require full simulation of the system, but for two qubits, we find that it takes a reasonable amount of time. Further details are given in App. D.
The protocol then proceeds the same as standard RB but with the qutrit-Clifford gates. After taking data with various sequence lengths, we estimate the RB decay parameter r, and calculate the subspace fidelity, which we denote with lower-case f to differentiate from the fullspace fidelity F , where d 1 is the dimension of the RB subspace. The subspace fidelity only provides an estimate for errors that are isolated to the RB subspace. We call the procedure given so far SRB-lite. An obvious problem with SRB-lite is that there may exist errors that are not isolated to the RB subspace. In a quantum algorithm, the input state to each gate is unlikely to be symmetric and only measuring errors on the symmetric subspace may miss errors that degrade the computation. Moreover, some errors might couple the RB and leakage subspaces, violating the assumptions we made to derive the RB decay curve. A similar problem arises in single-qubit RB if there are errors that cause population to leave the computational subspace, called leakage errors. Previously, several protocols to measure leakage rates were proposed, which we refer to as leakage benchmarking (LB) [18][19][20]. We now show LB can be performed along with SRB-lite to create a more complete picture of the errors in the MS gate. We call this combined procedure SRB.
In Ref. [19], leakage errors are quantified by a leakage rate L (rate that which population moves from RB to leakage subspace) and seepage rate S (rate that which population moves from the leakage to the RB subspace), where 1 1 and 1 2 are projectors onto the RB and leakage subspaces respectively. Our definitions differ slightly from Ref. [19] by factors d 1 (dimension of RB subspace) and d 2 (dimension of the leakage subspace) in order to simplify the derivations in App. C. Leakage and seepage errors cause additional terms to enter into Eq. (1) [18,19], where A, B and C relate to SPAM errors, r is the standard RB decay, and t relates to leakage and seepage rates.
We use a different technique that more effectively extracts the decay rates r and t than previous LB proposals. It requires including an additional random Weyl gate Q k , which can be compiled into the final gate C inv = Q k C † . . . C † 1 . This method was also proposed in Refs. [9,22] to catch pathological errors in standard RB. A similar method was explored in Ref. [12] to reduce the fitting parameters in order to derive multiplicative precision in standard RB. In addition to both of these advantages, we find that it also allows us to separate the analysis of the standard decay r from the effects of leakage errors and reduce the assumptions required in the derivation. Further details are given in App. C.
We again randomize over all sequences, which includes randomization over the final Weyl gate. Then, we process the data in two ways: • Standard analysis: For the standard RB decay, we only use the measurement outcomes m k such that ⟨m k Q k 0⟩ 2 = 1 in the absence of errors. Averaging over Weyl gates, with the related outcome m k , fixes B = 1 d 1 and C = 0, leaving a single exponential decay to fit, This differs from Ref. [19] since B is fixed but does not require assumptions on SPAM errors like in Ref. [20].
• Leakage analysis: To measure the leakage decay, we pick m = 0 outcomes, corresponding to the survival with the additional Weyl gate. Averaging over all sequences gives the average leakage decay In Ref. [19] L and S are used to quantify leakage errors, but we propose that t is more reliable and almost as valuable. This is because the estimates of L and S rely on the fits of C and B, which are limited by finite sampling and SPAM errors. The estimation method for t is almost identical to the estimation of r, and so follows the standard RB paradigm. The downside of only studying t is that we cannot differentiate between leakage and seepage errors, but our analysis finds that these rates are usually similar for standard errors in the trapped-ion system.
The results of both analyses then contribute to an estimate of the subspace fidelity of the qutrit-Clifford gates on the RB subspace. The leakage rate contributes to this estimate since the leakage errors reduce population in the RB subspace. We propose a variant of subspace fidelity that includes the seepage rate and is only dependent on r, t, and the dimensions, which we call the extended-sub fidelity, where the second equality is specific for the trapped-ion system. This quantity has the similar advantages of estimating t over L and S as well as capturing additional information about seepage errors in the fidelity.  2. Connections between measured quantities in SRB (standard decay rate r and leakage decay rate t), the extended-sub fidelity, and the full fidelity. The top row shows the connections between SRB and qutrit-Clifford gates. The bottom row shows the connection between interleaved SRB and the phase gate. If we assume the phase gate errors are depolarizing and dominate each Clifford gate then we can connect the top and bottom rows by the approximation To extract the extended-sub fidelity of a single phase gate we scale f ′ Cliff by the number of phase gates per Clifford, which is three, A similar procedure can be applied with f Cliff from SRBlite. This procedure is only exact when phase-gate errors dominate and each phase-gate has a depolarizing error channel. It is very likely that the phase gate errors dominate but unlikely that the errors are purely depolarizing in any real experiment. The quantity still serves as a rough estimate of the contribution of each phase gate to the error in the qutrit-Clifford gates. We report the estimated phase-gate error, which is 1 − f ′ ZZ . Another option to estimate phase gate fidelity is to perform interleaved randomized benchmarking (IRB) [32] to bound the phase-gate error. This is shown as the bottom row of Fig. 2. For SRB, such an interleaved procedure is more complicated since the phase gate is not a qutrit-Clifford gate. Instead of randomizing over the qutrit-Clifford gates we randomize over the larger group of arbitrary SU(3) unitary operators. The procedure is then the same as standard IRB; perform SRB and then perform SRB with interleaved phase gates. The standard IRB procedure also assumes that the errors are depolarizing and returns weak bounds on the fidelity if this is not true [32]. While this could be improved by including other experiments (like unitarity benchmarking [33]), we opt for the simpler procedure given above for estimating the phase-gate error and stress that the estimates are subject to the depolarizing assumption.

IV. EXTRACTING INFORMATION FROM SRB
To understand how errors affect SRB, we preformed an extensive study of typical errors in trapped-ion experiments. We found that each error falls into one of two categories: (1) Errors that preserve the subspaces and (2) errors that mix population between RB and leakage subspaces. The results are summarized in Table I and further details of each error channel are given in App. G.
Table I also shows how errors have different characteristics that are measurable by SRB. For example, intensity errors fall into type-1 since they preserve the subspaces, and therefore we expect an SRB experiment to measure r < 1 and t ≈ 1. Spontaneous emission causes mixing between the subspaces so that we expect r < 1 and t < 1. An interesting result from the study is that several errors can be tuned between type-1 and -2 by adjusting ion spacing and the laser/motional phase. This could also be used to further isolate dominant errors in the phase gate.
While SRB can detect most of the common errors in the phase gate, it will miss one type of error. We refer to these as 'worst-case' errors, and they correspond to errors that change the relative phase between the RB and leakage subspaces. For example, if the ions reordered this would look like a logical SWAP error, which generates a −1 phase between the two subspaces. This type of error will go unnoticed in an SRB experiment but could cause serious problems in a full computation. In general though, these errors are unlikely as they do not arise from any of the known MS error mechanisms [34] and our idling SWAP rate, (from background gas collisions), is observed to be on the timescale of tens of minutes.

V. EXPERIMENTAL PLATFORM
We implement our SRB procedure in a linear RF Paul trap with a chain of two 171 Yb + data qubits and two 138 Ba + ions for sympathetic cooling. The motion of the ions is described by shared motional modes with transverse (X,Y) and axial (Z) harmonic confinement at frequencies {ω x , ω y , ω z } = 2π × {2.72, 3.24, 1.15} MHz for a single 171 Yb + . The axial modes of motion and the three lowest frequency radial modes of the four ion chain are cooled ton < 1 via resolved sideband cooling of both Yb + and Ba + before the start of an SRB sequence [35]. The axial in-phase mode is cooled during an SRB sequence using sideband cooling on the Ba ions to mitigate effects of motional heating that takes place during the SRB sequence.
The data qubits are encoded in the 2 S 1 2 hyperfine levels of the 171 Yb + and initialized and read out using standard techniques [26]. The qubit levels are coupled via Raman transitions driven by a pair of 375 nm phaselocked lasers with a 12.6 GHz frequency offset to match the hyperfine frequency ω 0 of 171 Yb + in a magnetic field of 5.2 Gauss.
Collective single-qubit operations are preformed with co-propagating, σ + -polarized Raman beams. Two qubit operations are realized via theσ φ gate in the phasesensitive geometry [25,34]. The ion chain is illuminated with lin⊥lin-polarized Raman beams focused to a spot size ≈ 35 µm and propagating at a relative angle of 90 degrees, such that the wavevector difference ∆ ⃗ k ∥ẑ. The beatnote frequency of the two beams satisfies δω = ω 0 ± ω z,4 ± δ, where ω z,4 is the highest frequency axial mode and 2π δ ≈ t gate = 30µs.
The detection is performed by state-dependentresonance florescence, which returns photon counts proportional to the number of ions in the 1⟩ (bright) state. For a detection duration of 400 µs, we detect an average of 9 photons from each ion in 1⟩ and 0.1 photons from each ion in 0⟩. Empirically, we find that the threshold measurement technique [26] results in large measurement errors for 1 and 2 excited ion outcomes, which degrades the RB fitting. Instead, we perform detector tomography and process the RB data to reduce the measurement error [36]. This allows us to reconstruct relative frequencies of the desired projection operators.
Uncertainties correspond to one-sigma confidence intervals estimated by a semi-parametric bootstrap method [22]. In this experiment, the qutrit-Cliffords were found with a slightly different method that resulted in an average of 3.056 phase-gates per qutrit-Clifford gate. The SRB-lite method does not require detector tomography to reconstruct the relative frequencies needed. This makes the experiment simpler but does not identify the leakage error or reduce the number of fit parameters by fixing the B parameter.
Next, we implemented SRB with sequences of qutrit-Clifford gates [3,17,40] each with 50 random sequences (results plotted in Fig. 3b, solid lines with open markers). Each sequence was repeated 50 times to build statistics. For the leakage analysis, we used the even-parity outcomes. The results were fit to the decay curves derived in Sec. III with the least-squares method and plotted in Fig. 3a. We found a standard decay of r = 0.9934(8) and a leakage decay of t = 0.985 (23). These estimates translate to an extended-sub fidelity for the qutrit-Cliffords of f ′ Cliff = 0.993(3) and a phase-gate error To validate that phase-gate errors dominate in SRB and SRB-lite, we preformed standard single-qubit RB and measured a Clifford error of We also note that the qutrit-Clifford gates have a very similar form to methods for generating arbitrary twoqubit unitaries, as demonstrated in Ref. [37], which makes the extended-sub fidelity for the qutrit-Cliffords a relevant metric as well.
We also experimentally explored how SRB behaves when we inject errors into the system. We chose to inject over-rotation errors in the MS gate which translate to over-rotations in the phase gate plotted in Fig. 3c. For these tests, we used the same sequence lengths, reps, and trials per sequence as in Fig 3b. The standard decay r decreases quadratically with the rotation magnitude while t remains constant. The measured values match closely with the dashed lines, which are from a theoretical model described further in App. F. The zero over rotation corresponds to the solid lines in Fig. 3b and the maximum over rotation value corresponds to the dashed line in Fig. 3b. Due to the close ion spacing we were not able to inject errors that cause leakage to the anti-symmetric subspace but we do study them numerically in App. F.
In certain cases, the leakage analysis may not be a reliable quantitative measure of leakage errors. This can happen if there are other leakage subspaces or there is not adequate control of the leakage subspace (see App. C). These effects violate the assumptions made in the derivation of Eq. (8) and cause other exponential decay terms to contribute to p L ( ). We also observe from the semiparametric bootstrap of the numerical and experimental tests that the t decay rate has much higher variance. This is due to a much shallower decay, which can make fitting difficult, and the large variance in the values of p L ( ). All of these issues will degrade the estimate of t.
However, the standard analysis is robust to all of these effects by including Weyl gates with the final inversion gate. This means the decay rate r is still a valuable quantity and the leakage analysis could be used for qualitative information about leakage. The same analysis could be applied to any LB experiment to still extract information about gate fidelity in the presence of unknown or uncontrollable leakage.

VII. GENERALIZED SRB
SRB can be generalized to many other systems that are restricted to a non-universal set of gates. In this section we give a brief description of the general SRB procedure and outline the requirements.
In general, the set of available gates in a system is generated by a small set of interactions described by Hamiltonians {H i }. We assume that each Hamiltonian can be implemented with arbitrary amplitude i.e. U i = exp[−iθ i H i ]. These Hamiltonians form a representation of a Lie Algebra under commutation. If that algebra is equal to su(d) for a d-dimenstional Hilbert space then the system is said to be controllable [29] and any unitary in SU(d) can be created with proper choices of {θ i }.
For SRB, we consider a set of Hamiltonians do not generate su(d). Instead they generate a representation of another Lie algebra g, which is a sub algebra of su(d).
To perform SRB, we identify a subspace of the Hilbert space with dimension d 1 where g forms su(d 1 ). (If there are multiple such subspaces it is best to pick the largest one to capture the most information about the system.) The system is then said to be subspace controllable, and we can generate any unitary in that subspace.
Identifying such a subspace may be a non-trivial task. For the trapped-ion example, we were aided by standard results on angular momentum addition in quantum mechanics, i.e. two spin 1 2 systems add to form spin-1 (triplet) or spin-0 (singlet) states. Such relations could be applied to similar systems with arbitrary dimension subsystems. However, there may be extra leakage subspaces, which complicate the leakage analysis. In these cases, the leakage analysis will only provide identification of leakage rates and not quantitative measurement of their magnitude. Luckily, we can still extract some information on the quality of the gates from the standard decay analysis with the inversion procedure discussed in Sec. III. We consider this case in App. III.
Once the required subspace is identified, we then generate a set of unitary gates that forms a two-design. For the procedure discussed in Sec. III, d 1 must be a prime or a power of a prime. This restriction is due to the construction of the Pauli or Weyl operators that are used to generate the Clifford group. If d 1 is not prime or a power of a prime, we could consider other unitary two-designs or randomize over Haar-random unitaries.
The selected gates are then generated by optimal control techniques [31]. Again, optimal control requires a numerical search with a full simulation of the unitary. This limits the d 1 to be small enough to simulate the corresponding unitary operations classically. While this is a limitation of SRB, standard RB is also limited to small system sizes (the largest published result is on three qubits [38]). This is for two reasons: (1) the size of the Clifford group grows exponentially with the number of qubits and the optimal decomposition for each Clifford is found by exhaustive search, and (2) the difference in depth between Cliffords blows up with qubit dimension, which violates the standard RB assumption that each Clifford has the same error.
The remainder of the general procedure then follows standard RB and the trapped-ion example considered above. Random gates are applied for different sequence lengths and the total decay rate is measured. We can again determine the standard and leakage decays to derive the extended-sub fidelity.

VIII. CONCLUSIONS
We proposed a new method for benchmarking entangling qubit gates in trapped-ion systems, which is experimentally easier to implement and provides more information than standard methods used in such systems. While the method does not capture as much information as standard RB, it can be combined with LB to identify most errors that occur. We also proposed a modified LB procedure that allows for separating the leakage decay from other error rates that should be useful with any LB procedure.
We implemented this technique on a trapped-ion system to measure the quality of entangling gates. We verified a low error phase gate with 1−f ′ ZZ = 2(1)×10 −3 . For our system, SRB is faster than both standard RB with two qubits (since it does not require transport) and Bellstate fidelity measurements (since it does not have standard quantum limit scaling with precision). This makes it an excellent method to use in the design-and-test cycle of tuning up high-fidelity entangling gates.

ACKNOWLEDGMENTS
The authors would like to thank Patty Lee for her helpful comments on an earlier version of this manuscript.

Appendix A: Python implementation
A Python implementation of SRB-lite, SRB, and detector tomography is also included as ancillary file to this submission. These methods are demonstrated in two Jupyter Notebooks that include simulation methods. They can also be edited to work with experimental implementation. Further details are given in the Jupyter Notebooks.

Appendix B: Standard randomized benchmarking derivation
In this appendix, we describe standard RB with a final Pauli/Weyl gate included. As discussed in the main text, including this gate (1) catches total amplitude damping errors [22], (2) reduces the fitting parameters [12,22], (3) allows one to show multiplicative precision [12], and (4) reduces the assumptions needed for LB.
We work in the Liouville (or transfer matrix) representation of quantum processes. In this representation, density matrices and measurement operators are translated to supervectors A → A⟫, and processes are translated to superoperators U [⋅]U † → U. The inner product of two supervectors is equal to the trace overlap between the corresponding matrices ⟪A B⟫ = Tr(A † B). A supervector may be unnormalized, ⟪A A⟫ = a; for example, a mixed state does not have a normalized supervector. However, at times we consider an orthonormal basis of supervectors, that is a basis of orthonormal operators, ⟪P i P j ⟫ = δ i,j for a d-dimensional Hilbert space. The superoperators are related to the Kraus decomposition by If A is a trace preserving (TP) process, then ⟪1 A = ⟪1 .
Each applied Clifford gate can be separated into its ideal action C i and an error channel Λ i . We also include the effects of a state preparation error channel Λ P and a measurement error channel Λ M . For zeroth-order RB we make the unrealistic assumption that Λ i = Λ for all i. Recent work has shown that RB is robust to this assumption, although the final estimate will no longer have the same relation to the process fidelity [13].
The extra Pauli/Weyl gate Q k is compiled with the final inversion gate at the end of each sequence. We then study outcomes m k such that ⟨m k Q k = ⟨0 , and define the survival probability as where Q k = Q * k ⊗ Q k , ψ⟩ ⟨ψ → ψ⟫ for any pure state ψ⟩. The superopertors {C i } form a representation of the Clifford group under matrix multiplication. Therefore, in each of the i sequences we expand C i,j in terms of relabeled elements {D i,j } defined as C i,1 = D † i,1 , C i,j = D † i,j D i,j−1 , and inversion gate C i,inv = D i, . Next, we average over all possible sequences i, which translates to an average over all Clifford gates at each position, yielding where N C is the number of Clifford gates, and the posi- which is a super-super operator, (an operator acting on a superoperator). The twirl T maps any process to the depolarizing processΛ d [10], which mixes an input state ρ with the identityΛ d [ρ] = rρ + 1−r d 1. The parameter r is the called the depolarizing parameter and relates directly to the process fidelity. In the Liouville representatioñ Λ d = rP + 1⟫⟪1 , where 1⟫ is the trace-normalized identity and P is the projection onto the orthogonal subspace 1⟫⟪1 . Substituting this superoperator into the previous equation gives where A k = ⟪m k Λ M ΛQ k PΛ p 0⟫ and B k = ⟪m k Λ M ΛQ k 1⟫⟪1 Λ p 0⟫. Next, we average over all d 2 values of Q k and define A = 1 d 2 ∑ k A k . The asymptote reduces where we used the relations by line number: by the definition of a POVM, and (2 → 3) ⟪1 Λ = ⟪1 for any trace preserving map Λ and ⟪1 0⟫ = d −1 2 . We then get the fixed-asymptote decay function Subgroups of Pauli/Weyl gates, e.g. {1, X}, will also fix the asymptote [12] but will require additional assumptions for the leakage analysis. We then estimate p( ) for different values of and fit the results to Eq. (B5) to estimate A and r. The process fidelity of the depolarized process is equal to the process fidelity of the unknown error channel [39].
In practice, we do not measure sequences with every possible final Pauli gate, but only sample from an even distribution over all final operators. This is similar to the standard RB practice of sampling from all sequences and, following similar arguments as in standard RB [11], we expect to approach the exact distribution with a reasonable number of samples.

Appendix C: Derivation of subspace randomized benchmarking decay
In this appendix, we repeat the previous derivation with leakage/seepage in the different erorr channels. We derive the decay curves in Eqs. (7) and (8). We first follow similar steps as the derivation for leakage benchmarking given in Ref. [19]. We then consider expansions of this derivation.

Basic decay derivation
Let d 1 be the dimension of the RB subspace, d 2 be the dimension of the leakage subspace, and d = d 1 + d 2 . The projector on the RB subspace is defined as 1 1 , with tracenormalized supervector 1 √ d1 1 1 → 1 1 ⟫, and the projector onto the leakage subspace is defined as 1 2 , with tracenormalized supervector 1 √ d2 1 2 → 1 2 ⟫. Consider an ideal Clifford unitary C acting on the RB subspace. Even in an ideal experiment, implementing this Clifford will also implement a unitary process V on the leakage subspace. Let U = C ⊕V be the total unitary, or in the Liouville representation U = (C ⊕V ) * ⊗(C ⊕V ).
We expand this superoperator, where C = (C ⊕ 0) * ⊗ (C ⊕ 0), V = (0 ⊕ V ) * ⊗ (0 ⊕ V ), and A = (C ⊕ 0) * ⊗ (0 ⊕ V ) + (0 ⊕ V ) * ⊗ (C ⊕ 0) (from now on we drop the direct sum with zero notation e.g. C ⊕0 → C for brevity). The total superoperator U is block diagonal since CV = CA = AV = 0. However, the error channels Λ, Λ P and Λ M may not be block diagonal in this basis. Some of these off-diagonal effects correspond to the leakage/seepage errors we wish to quantify. Each SRB sequence i consists of + 1 gates U i,j , which can be expanded as shown above. The survival probability is There are a few differences between this survival probability and Eq. B2. For one, we do not fix the outcome m with respect to the final Pauli/Weyl gate. This will later allow us to do the leakage analysis discussed in the main text. Also, the measurement projectors {⟪m } are slightly more complicated since we now include the effects of the leakage subspace. The definition of POVMs implies that some POVM elements will include projection onto the leakage subspace, ∑ m ⟪m = √ d⟪1 where 1 is the identity on the entire system. For example, in the trapped-ion system the odd-parity measurement projects onto both the symmetric and anti-symmetric states.
To simplify the derivation, Ref. [19] makes the following assumptions about the experiment: • Assumption 2 : A i,j = 0, which can be enforced by averaging over the relative phase between subspaces.
• Assumption 3 : The ability to independently average over C i,j and V i,j .
• Assumption 4 : The action on the leakage subspace {V i,j } forms a unitary one-design.
These assumptions require some level of control over the leakage subspace, which is possible in the trapped-ion system, but may be difficult in other systems. For now, we make the same assumptions but consider their relaxation later. The result is that each gate is expanded U i,j = C i,j +V i,j where C i,j and V i,j can be varied independently. As in standard RB, we expand each RB subspace superoporator C i,j = D † i,j D i,j−1 , and additionally each leakage sub- The final gate becomes more complicated when considering the additional leakage subspace because it contains both an inversion and the final Pauli/Weyl gate. We associate any action on the leakage subspace with the inversion part since the action on the RB and leakage subspace can be varied independently. The inversion gate is selected to relate the sequence to a series of twirls C i,inv = D i, , but in the leakage subspace there is an additional term when we expand We combine this with the extra term on the leakage subspace Q ′ k,i = Q * k ⊗Q k +W † i , where i, + 1 → i since it is a random element from the set of leakage subspace gates. Now, we are ready to average over all possible sequences where N V is the number different gates applied to the leakage subspace and For this analysis, we define a new twirl super-super operator and apply the one-and two-design definitions where P 1 is the projection onto the part of the RB subspace orthogonal to 1 1 ⟫⟪1 1 . The coefficients are defined r = Tr(P 1 Λ), L = ⟪1 2 Λ 1 1 ⟫, S = ⟪1 2 Λ 1 1 ⟫, Our definitions of the leakage and seepage rates L and S differ slight from Ref. [19] by factors of d 1 d 2 .
The next step is to solve for T [Λ] by diagonalizing T [Λ]. The first term in Eq. (C4) is orthogonal to all other terms P 1 1 i ⟫⟪1 j = 0 for i, j ∈ {1, 2}, but the remaining terms span a two-dimensional subspace { 1 1 ⟫, 1 2 ⟫} decomposed as with eigenvalues λ ± and eigenvectors Π ± , Then T [Λ] = r P 1 + Π + + λ − Π − . In the main text, we defined t = λ − . We are now ready to combine our reductions to derive the survival probability with leakage. Substituting in the expressions for T [Λ] , We refer to the first term as the standard decay, the second as the asymptote, and the final term as the leakage decay.
To measure the standard decay (standard analysis in the main text) we study the outcomes m = m k defined by the final Pauli/Weyl gate ⟨m k Q k 0⟩ = 1, and average over all d 2 1 Pauli/Weyl elements. This is the same approach as in App. B. The coefficient for the standard decay cannot be simplified further so define A = 1 d 2 1 ∑ k ⟪m k Λ M ΛQ k PΛ P 0⟫. We can reduce the asymptote in the following steps: where we used the relations by line number: (1→2) Q k Π + = Π + ∀ k and then the only k dependence remains is in the outcomes ∑ k ⟪m k = ⟪1 d 1 √ d, (2→3) ⟪1 Λ M ΛΠ + Λ P = ⟪1 when Λ M , Λ, and Λ P are TP on the full space and ⟪1 Π + = ⟪1 by inspection of Eq. (C5). A similar analysis can be applied to the leakage decay by noting Q k Π − = Π − and ⟪1 Π − = 0 to show that This means that the method derived in App. B is robust to leakage/seepage errors.
To quantify leakage/seepage (leakage analysis in the main text), we return to the distribution in Eq. (C7) and study the outcome m = 0, which is chosen to have no projection onto leakage subspace. For example, in the trapped-ion system, the even-parity outcomes are possible choices, whereas odd-parity outcomes have support in the anti-symmetric state. We then average over the Pauli/Weyl gates, which is equivalent to character benchmarking [17]. The Pauli/Weyl operators form a unitary one-design and project onto the trivial irreducible representation Π Q = 1 However, in this case the standard decay term is equal to zero since Π Q P 1 = 0. The asymptote is easily calculated in the subspace spanned by { 1 1 ⟫, 1 2 ⟫}. The ideal projection of the prepared state in this subspace is 1 √ d1 1 1 ⟫. State preparation and measurement each may have leakage/seepage errors which cause them to overlap with the leakage subspace, where L P and S P are leakage and seepage rates for state preparation, respectively, and L M and S M are leakage and and seepage rates for measurement, respectively. The final step is to note that Π Q is the identity in this reduced subspace so Π Q Π ± = Π ± . Then the asymptote is reduced to where B = O(L M , S M ). The leakage decay coefficient is found similarly, but we omit the full definition for brevity and instead give the simplified form where C = O(L M , S M , L P ) and is related to SPAM errors. The net decay function is then, In practice, we only fit the measured decay with the first two terms since the final two terms are likely much smaller with high-fidelity state preparation and measurement.

Trapped-ion specific leakage detection
In the trapped ion system d 1 = 3 and d 2 = 1, and Assumptions 3 and 4 are automatically enforced since W i,j trivially forms a unitary one-design W i,j = 1 2 ⟫⟪1 2 ∀ i, j. Assumption 2 can be enforced by creating a phasereversed set of qutrit-Cliffords gates by including one more phase gate (further details given in App. D). Averaging over both the standard qutrit-Cliffords gates and the phase-reversed qutrit Cliffords then eliminates the A i,j terms. Alternatively, in App. C 3 we consider the derivation above but with Assumption 2 relaxed.
Additionally, trapped-ion systems have more leakage sublevels that are not treated in the previous deriviation. For 171 Y b + there are two other hyperfine levels per qubit ion that population may leak to due to sponetaneous emission. We consider such effects on the SRB derivation in App. C 4

Averaging without the phase-reversed set
In the experiment, we chose not to average over the phase-reversed qutrit-Cliffords gates, discussed in the previous subsection, in order to fix the number of phase gates per Clifford. This choice is experimentally convenient since it means each Clifford has a very similar implementation. It also makes the experiment closer to the standard RB assumption that each Clifford has the same error channel. However, it means that Assumption 2 is not enforced. Ref. [19] suggested that the randomization over the relative phase was not necessary and A i,j ≈ 0. However, we find in this case that A i,j is not small and in fact has some eigenvalues equal to one. This is an extreme case where the relative phase is the same for each Clifford gate since it is only changed by phase gates and each qutrit-Clifford gate is compiled with exactly three phase gates.
We repeat the analysis for the trapped-ion system but without Assumption 2 to verify that the our procedure will still work. Assumptions 3 and 4 are still enforced since W i,j = 1 2 ⟫⟪1 2 ∀ i, j.
We return to Eq. (C3) but include A i,j . We expand the extra term Including these extra terms gives, where the second line contains the terms we considered earlier and the third line contains the new terms.
We numerically tested the effect of each term in Eq. C14 by forming a matrix representation of the supersuper operator, e.g.
where ⋅) is the vectorized superoperator. Next, we take the eigendecomposition and analyze the magnitude of the eigenval-ues. We see that all new terms except for the final term have small eigenvalues (≤ 0.1). The small eigenvalues are exponentially suppressed in and we do not expect them to contribute to the decay curve. These are from an odd number of tensor products of D i unitaries, which average out to be small. The final term has an even number of D i unitaries, We see the final two terms have small eigenvalues (although slightly larger than before, ≤ 0.2). The first two terms, however, each have a single eigenvalue equal to one. The corresponding eigenvectors are can be derived in a swapped basis, where a⟩ = 1 √ 2 ( 00⟩ + 11⟩). The resulting operator projects the error channel ∑ j,k ⟪k, a Λ k, a⟫ j, a⟫⟪j, a = Tr(ΛΠ upper )Π upper where Π upper = ∑ k k, a⟫⟪k, a . A similar derivation can be applied to the second term giving Tr(ΛΠ lower )Π lower where Π lower = ∑ k a, k⟫⟪a, k . These two projectors correspond to parts of Λ that map coherences between the RB and leakage subspace since k, a⟫ → k⟩⟨a .
The Pauli/Weyl gate Q ′ k has the same product with each of the previous projectors since it still acts as the identity over the { 1 1 ⟫, 1 2 ⟫} subspace. The only change is that it has the additional B ′ † k term, which overlaps with Π upper and Π lower . Further analysis of the extra term may show ways to eliminate these contributions.
The rest of the analysis is the same as above but with these two extra terms. The new terms do not contribute to the standard analysis since ⟪1 Π upper/lower = 0. Unfortunately, in the leakage analysis the new terms persist. This means the leakage decay derived in Eq. (C9) has two extra decay rates from Π upper/lower . These terms do not relate to leakage/seepage errors and we cannot differentiate the effects. However, we expect that the state prepared and the measurements used have very little coherence between the two subspaces, which in practice likely limits the contribution of these two terms. This requires additional assumptions, but we believe they are reasonable at least in the trapped-ion case.

Considering more leakage subspaces
For systems that have n ≥ 2 leakage subspaces the analysis becomes more complicated, and we are not able to extract the leakage decay exactly. For this sub-appendix we enforce Assumptions 2-4, which now require additional control of the phase between each leakage subspace. Then T [Λ] can be reduced to a similar form as in Eq. (C4), with the addition of leakage L i,j and seepage S i,j rates between the n + 1 subspaces, With n leakage subspaces, we must consider a larger operator subspace spanned by { 1 i ⟫} and diagonalize the corresponding (n + 1) × (n + 1) matrix. We can identify one eigenprojector under the assumption that Λ is TP, which implies T [Λ] is TP. We label this eigenprojector Π 1 = V 1 ⟫⟪1 where V 1 is a dual vector to ⟪1 , and the corresponding eigenvalue is λ 1 = 1 since ⟪1 T [Λ] = ⟪1 . (In the previous subsection this was Π + .) The remaining eigenvalues λ i for i = 2, . . . , n+1 are complicated functions of the leakage/seepage rates between subspaces, and have corresponding eigenprojectors Π i . We then substitute into Eq. C3, The standard analysis remains the same since ⟪1 Π 1 = ⟪1 and ⟪1 Π i>1 = 0. The asymptote corresponds to the known eigenvalue λ 0 and is again B = 1 d 1 .
Things do not reduce as nicely for the leakage analysis. Eq. (C18) has n+1 terms but, as in the two-subspace case, the standard decay is zero since Π Q P 1 = 0. The asymptote is nonzero and a function of the leakage/seepage rates. There are additionally n leakage decay parameters, which would require a multi-exponential fitting method to extract. Even if we use such a method, the different eigenvalues that define the leakage/seepage rates between subspaces could not be discerned.
This contributes another two terms in the summation for full fidelity. The extended-sub fidelity is an average of these d 2 1 + 1 terms. To weakly bound the full fidelity, we bound the magnitude of the other d 2 −d 2 1 −1 terms, which come from basis elements that mix population between the RB and leakage subspaces. These terms are each in the range [−1, 1], since the operator subspace is spanned by a set of traceless Hermitian operators with eigenvalues ±1. This gives the bounds, The lower bound is proportional to the extended-sub infidelity, up to dimensional constants, which means it approaches zero for high-fidelity gates. However, the upper bound has a constant factor, which means there may be a process that has a high extended-sub fidelity but low full fidelity. The upper-bound is saturated when all other terms in the full fidelity are equal to -1. Consider the operator basis elements P i,1 = i 1 ⟩⟨i 2 and P i,2 = i 2 ⟩⟨i 1 , where i 1 ⟩ are basis elements for the RB subspace and i 2 ⟩ are basis elements for the leakage subspace. If ⟪P i,1 Λ P i,2 ⟫ = −1 for all i then Λ generates a -1 relative phase between the RB and leakage subspace. This is the worst-case error we considered in Sec. IV corresponding to SWAP errors.

Appendix F: Simulations with different error channels
In this appendix, we compare analytic models of different error channels with simulations of SRB. We focus on the estimated quantities r, t, and f ′ ZZ to validate the expected behavior with the selected error channels.

Analytic models
While the extended-sub fidelity only weakly bounds the full fidelity, we find that for many typical error channels the values are close. In this section we directly calculate the SRB decay parameters and difference between the extended-sub fidelity and full fidelity, for intensity errors, optical pumping, and inhomogeneous fields. We chose these error channels since they are parameterized by a single error magnitude and are easy to simulate, while still demonstrating different behaviors. The results of the analysis are shown in Table II. Below is a brief description of each: • Intensity errors: Errors that cause the implemented phase gate to have a different amplitude than the expected value of π 4. This could be due to laser intensity fluctuations, Debye-Waller errors, or uneven intensity on the MS gate beams. We model these errors by applying an extra gate U e ( ) = exp[−i ZZ] after each phase gate.
• Optical pumping: Errors that cause population to be pumped deterministically to one state (in our case modeled as 00⟩). This could be due to stray light interacting with both ions. We model the errors as an amplitude damping channel with rate .
• Inhomogeneous fields: Errors due to stray magnetic fields that fluctuate slowly with respect to the gate time. These errors could be due to trap electronics or other magnetic field sources in the lab. They are modeled as a phase-dampening channel with rate .

Simulations
We studied SRB's performance with the three error channels described above. The results are plotted in Fig. 4. We simulated SRB experiments with injected errors in the phase gate with 1 − F ∈ [10 −3 , 10 −1 ], singlequbit depolarizing errors of magnitude 10 −5 per gate, and depolarizing SPAM errors of magnitude 10 −3 . For most values of 1 − F the sequence lengths are exponentially distributed up to 1 (1 − F ). For larger values of 1 − F , we used the sequence lengths 1, 2, ..., 10, which perform better. For each sequence length, we generated 100 random sequences of gates and sampled each of those sequences 100 times to build statistics. We simulated the trappedion parity measurement, which detects the number of ions in the bright state by state-dependent-florescence, and performed detector tomography to better differentiate the outcomes. For the leakage analysis, we used the even-parity outcomes. Uncertainties are calculated by the semi-parametric bootstrap method [22] with onesigma confidence intervals plotted. An example of a single SRB experiment is plotted in Fig. 4b with fits to the standard decay model (solid line) and leakage decay model (dashed line). Fig. 4a shows the full infidelity vs. the estimated extended-sub fidelity of the phase gate from SRB. For the errors simulated, we see that the extended-sub infidelity estimated by SRB is very close to the full infidelity.
Figs. 4c-d, shows the extracted standard and leakage decay rate from each experiment with (c) intensity errors, (d) optical pumping, and (e) inhomogeneous fields. We calculated the expected values of r ZZ (solid lines) and t ZZ (dashed lines) based on the error models described above. We then scaled the measured valuesr andt from SRB assuming that the phase errors dominate and are depolarizingr ZZ =r 1 3 andt ZZ =t 1 3 . This assumption is clearly not true since we inject error channels that are not depolarizing, but it allows us to relate the models for phase gate errors directly to SRB measurements. We observe that the values ofr ZZ andt ZZ are similar to that predicted from the analytic model described above. The intensity error hasr ZZ decay but not ZZ decay while the optical pumping and inhomogeneous magnetic field errors have visible decay for both.
For larger error magnitudes the value oft ZZ differs more from the expected value. This is due to fitting troubles in the leakage analysis that are discussed more in Sec. VI but the analysis still gives an indication of leakage errors. In this appendix we categorize some common error channels in the phase gate according to their action on the RB (symmetric) and leakage (anti-symmetric) subspaces. To determine if the error channels are subspace preserving (type-1) or leakage inducing (type-2) we look at the leakage and seepage rates defined in the App. C. When L = S = 0 then the error is type-1 otherwise the errors are type-2.
The phase gate is created by applying wrapper pulses before and after an MS gate, and, therefore, most errors we study are MS gate errors. Let U M S be the unitary that describes the ideal MS gate. Within the Lamb-Dicke approximation and the resolved side-band limit, the propagator has a simple closed form that can be found by considering the Magnus expansion ]. The first term in this expression is usually referred to as the spin-dependent force and the second is referred to as the spin-dependent phase. Ideally, the spin-dependent force integrates to zero and one is left with U M S = exp[−iŜ 2 ( ηΩ 2δ ) 2 (δτ g − sin(δτ g ))] whereŜ is the collective spin operator, η is the Lamb-Dicke parameter, Ω is the Rabi frequency, δ is the detuning from resonance and τ g is the gate time.
Intensity/Lamb-Dicke/Debye-Waller errors: Fluctuations in the Rabi frequency Ω can usually be attributed to either laser intensity noise, the Debye-Waller effect, or the Lamb-Dicke nonlinearity [34]. These errors re-scale the rotation angle of the spin-dependent phase, but the operatorŜ 2 remains symmetric under qubit exchange, implying that these errors are subspace preserving.
Trap frequency error : Trap frequency fluctuations are another common error source [40], and can be thought of as fluctuations in the detuning parameter δ. While these errors manifest themselves as fluctuations in the spin-dependent phase, they also alter the structure of the spin-dependent force component. Ideally the spin-dependent force component is given by exp[−i ∫ dt ′ H(t ′ )] = exp[Ŝ(α * (t)a † − α(t)a)] and α(t gate ) → 0. However, trap frequency fluctuations will lead to a non-zero α andŜ is generally not symmetric with respect to qubit exchange [25]. In order forŜ to be symmetric under qubit exchange, the laser phase on the two different ions must be the same as their relative motional phase in the motional eigenmode used for the gate. For example, if one opts to use the center-of-mass mode for the gate operation,Ŝ is only symmetric if the ion-spacing x 0 satisfies kx 0 = mπ where m is even with k being the laser's wave-vector component that is parallel to the relevant principle axis of motion. If one chooses to use a motional mode where the two ions' motion is out-of-phase, then m is required to be an odd integer. Therefore trap frequency errors can be either subspace preserving or leakage inducing.
Spectator mode coupling: Spectator mode coupling is close in structure to trap frequency fluctuation errors. The inclusion of spectator modes into the propagator changes the amount of spin-dependent phase that is accumulated, but also introduces additional spin-dependent forces. In the case of a 4-ion crystal, no matter which motional mode is used for the gate, there will always be spectator modes of both types, those where the qubit ions' motion is in-phase and those where it is out-ofphase. This implies that no matter what the ion-spacing is, there will be residual spin-dependent force operators that are leakage inducing.
Carrier coupling: Another off-resonant transition to be considered is that of the so-called carrier transition, which changes the spin without changing the motional state. As discussed in [34], this error channel can be made to vanish if the gate time is chosen to be commensurate with the mode frequency, but will be non-zero in the presence of small timing errors. Unlike spectator-mode transitions, this error will introduce operators whose exchange symmetry only depends on the laser phase at the location of the ions and is subspace preserving when m is even. Therefore, carrier coupling can be either subspace preserving or leakage inducing.
Carrier frequency fluctuations: Fluctuations of the carrier laser frequency are equivalent to collective qubit frequency fluctuations. In general, qubit frequency fluctuations are leakage inducing if there are differential shifts on different qubits. However, in the case where the two data ions are only 2.5µm apart, common-mode field fluctuations should dominate, and therefore are subspace pre-serving to a good approximation.
Uneven illuminations: Counter-intuitively, uneven illumination of the two ions is also subspace preserving. Under uneven ion illumination, the collective spin-operator gets perturbed asŜ =ŝ 1 +ŝ 2 → (1+δΩ 2)ŝ 1 +(1−δΩ 2)ŝ 2 , which is not symmetric under qubit exchange. However, the spin-dependent force still integrates to zero in this case and the square of this operator is still symmetric under qubit exchange, meaning that this error is subspace preserving.
Heating: Heating of the motional mode can be thought of as instantaneous momentum kicks randomized in time. These momentum kicks do not prevent the phase space loops from closing, rather they only change the amount of area that is traced out. Therefore, this error channel is subspace preserving.
Spontaneous emission: Spontaneous emission drives population outside the computational subspace and is therefore leakage inducing.
In order to generate Table I which summarizes which phase gate errors are subspace preserving, we go through the same analysis for the wrapper pulses and then XOR those results with those from the MS gate. In our case, where the wrapper pulses are collectively applied, the only single-qubit gate errors that couple the subspaces come from uneven illumination of the two qubits, which makes this error leakage inducing.