Leakage Benchmarking for Universal Gate Sets

Errors are common issues in quantum computing platforms, among which leakage is one of the most-challenging to address. This is because leakage, i.e., the loss of information stored in the computational subspace to undesired subspaces in a larger Hilbert space, is more difficult to detect and correct than errors that preserve the computational subspace. As a result, leakage presents a significant obstacle to the development of fault-tolerant quantum computation. In this paper, we propose an efficient and accurate benchmarking framework called leakage randomized benchmarking (LRB), for measuring leakage rates on multi-qubit quantum systems. Our approach is more insensitive to state preparation and measurement (SPAM) noise than existing leakage benchmarking protocols, requires fewer assumptions about the gate set itself, and can be used to benchmark multi-qubit leakages, which has not been achieved previously. We also extended the LRB protocol to an interleaved variant called interleaved LRB (iLRB), which can benchmark the average leakage rate of generic n-site quantum gates with reasonable noise assumptions. We demonstrate the iLRB protocol on benchmarking generic two-qubit gates realized using flux tuning and analyzed the behavior of iLRB under corresponding leakage models. Our numerical experiments showed good agreement with the theoretical estimations, indicating the feasibility of both the LRB and iLRB protocols.


Introduction
Quantum computation maps information processing into the manipulation of (typically microscopic) physical systems governed by quantum mechanics.Although quantum computation holds the promise to solve problems that are believed to be classically intractable, practical quantum computation suffers from various noise sources, ranging from fabrication defects and control inaccuracies to fluctuations in external physical environments.Such noise greatly hinders the practicability of quantum computation on unprotected, bare physical qubits beyond proof-of-concept demonstrations.
While any kind of error is unwanted and would possibly affect the quality of the computation processes, there is a significant difference between the harmfulness of different types of errors.The most "benign" error happens locally and independently on single qubits; such errors can, in principle, be compressed arbitrarily with quantum error correction under reasonable assumptions on the error rates [1,2].More malicious errors might introduce time correlations (e.g., non-Markovian errors) or space correlations (e.g., crosstalk) and are more challenging to mitigate.Of particular interest is the leakage error, where a piece of quantum information escapes from a confined, finite-dimensional Hilbert space used for computation, called computational subspace, to a leaked subspace of a larger Hilbert space.Such escaped information might undergo arbitrary and uncontrolled processes and is harder to detect, let alone correct.More seriously, typical frameworks of quantum error correction only deal with errors happening within the computational subspace and are either unable to apply or scale poorly with the leakage error.It is thus of great importance to be able to detect, correct, or even suppress leakage errors in order to conduct large-scale quantum computation.
This paper focuses on estimating the leakage error rate associated with a given quantum processor, preferably efficiently and accurately.This task is part of a process usually referred to as benchmarking, provides an estimate of certain characteristics of a piece of the quantum device before proceeding with subsequent actions.In the context of leakage benchmarking, the information can be used as a criterion to accept or abort a newly-fabricated quantum processor or as feedback information on leakage-suppressing gate schemes.
The first theoretical framework for RB-based leakage benchmarking was given by Wallman et al. [22].Without any prior assumption on the SPAM noise, this protocol was able to provide an estimate for the sum of the leakage rate and the seepage rate, i.e., the rate information in the leaked subspace comes back to the computational subspace.
Refs. [8,23] later gave a detailed analysis of the protocol and illustrated this framework with several examples relevant to superconducting devices.The authors were also able to differentiate the leakage from the seepage with reasonable assumptions on the SPAM noise.Based on these protocols, several experimental characterizations of single qubit leakage noise have been proposed in superconducting quantum devices [24,25], quantum dots [26], and trapped ions [27].
There are two major limitations to the existing protocols [8,11,22,23].First, all protocols require that the quantum gates act nontrivially on the leakage subspaces, in order to eliminate non-Markovian behavior originating from residual information stored in the leakage subspace.As most practical gate schemes only focus on their actions on computational subspaces rather than the leakage subspaces, leakage benchmarking schemes built upon them typically do not work in general, multi-qubit quantum systems.Second, most existing protocols can only estimate the sum of the leakage rate and the seepage rate without prior knowledge of SPAM noise, and the SPAM information is required if we need to get the leakage and seepage rates separately.As there is typically only one set of state preparation and measurement within one run of benchmarking, the SPAM errors do not get amplified and cannot be measured accurately [28].Such inaccuracy would further affect the accuracy of gate leakage rate estimation.A natural question arises: How can one characterize the leakage rate of a multi-qubit system without operating the leakage subspace while maintaining robustness to SPAM noise?
In this paper, we propose a leakage benchmarking scheme based on RB, dedicated to benchmark leakage rates on multi-qubit systems.Compared to existing protocols requiring the leakage subspace to be fully twirled, our scheme only requires having access to the Pauli group with gate-independent, time-independent and Markovian noise.Assuming each qubit has only one-dimensional leakage space, such a gate set does not twirl the leakage subspace as a whole, but instead twirls each invariant subspace of the Pauli group individually.This allows us to formulate the LRB process as a classical Markovian process between different invariant subspaces, which can be described by a Markovian Q-matrix [29].The leakage and seepage rates of the system can then be estimated by leveraging the spectral property of the Markovian process, which can, in turn, be estimated similarly to RB protocols on the computational subspace.
The Q-matrix has a dimension exponential with respect to the number of qubits in general and thus the spectral property is hard to be measured using LRB experiments.To further simplify the problem, we study the spectrum of the Q-matrix in two physically-motivated scenarios: The first model, named as leakage damping noise, assumes that leakage happens at most one qubit, and leakage does not "hop" from one qubit to another, which is the generalization of amplitude damping noise [30] in the computational subspace; the second model assumes that each qubit undergoes an independent leakage process.In both cases, the spectral property of the Q-matrix can be significantly simplified, and easier for data analysis.We also show how to calculate the corresponding average leakage rates on the above two noise scenarios of the proposed LRB protocol.As an illustration of the leakage damping noise model, we found the noise model of commonly used two-qubit gates such as iSWAP, SQiSW, and CZ gates all belong to this form.
Building upon the foundation of leakage randomized benchmarking (LRB) protocols, we delve deeper into the study of leakage benchmarking for specific multi-qubit gates, which is a crucial aspect of quantum hardware development.To this end, we propose an interleaved variant of the LRB (iLRB) protocol that allows for the benchmarking of individual gates, rather than a set of gates.We show that leakage rate can be extracted in general for arbitrary target gates with access to noiseless Pauli gates, and perform more careful analysis when Pauli gates are implemented noisy.In addition, we show that the leakage rate of the target gate can still be extracted under certain physically-motivated assumptions that inherently apply to flux-tuning gates in superconducting quantum computation.To demonstrate the applicability of the iLRB protocol, we apply it to the case of flux tunable superconducting quantum devices [31], construct its noise model, and benchmark the leakage rate of the iSWAP gate.
This paper targets both theorists and experimentalists, as it seeks to establish an experimental-friendly leakage benchmarking scheme.We offer a thorough theoretical analysis for multi-qubit scenarios, as well as numerical verification of the average leakage rate for the iSWAP gate.This is achieved by extracting the noise model of the iSWAP gate from its Hamiltonian evolution.
In Section 2, we introduce the fundamental concepts and notations.Section 3 presents our LRB protocol and analyzes the calculation of the average leakage rate using this method.In Section 4, we provide a detailed examination of the average leakage rate under two leakage models: single-site leakage and no cross-talk.Section 5 proposes the iLRB protocol for any target gate that commutes with the noise channel, focusing on a special leakage damping noise.In Section 6, we numerically validate the LRB and iLRB protocols.Additionally, we introduce the leakage damping noise model for iSWAP/SQiSW/CZ gates in flux-tunable superconducting quantum devices, based on their Hamiltonian evolution.We also test the iLRB protocols numerically using the noise model of the iSWAP gate.Finally, Section 7 concludes the paper with a discussion of our work and suggestions for future research directions.

Notations
In order to characterize leakage, we assume that the quantum states lie in a Hilbert space H with finite dimension d that decomposes into a computational and a leakage subspace, denoted as H c and H l respectively.Let d c := dim (H c ) and d l := dim (H l ) = d − d c be the dimensions of H c and H l .Unless explicitly specified, we assume throughout the paper that a single qubit (site) lies in a three-dimensional Hilbert space with basis {|0⟩, |1⟩, |2⟩}, where the computational subspace is spanned by {|0⟩, |1⟩} and the leakage subspace by {|2⟩}.In other words, higher-level excitations of a qubit can be ignored.We call such a system a single qubit with leakage.
A composite system of n qubits with leakage lies in a Hilbert space , where H c k (H l k ) represents the computational (leakage) subspace of the qubit k.We define the computational subspace of H be where no qubits leaks, that is, H c = n k=1 H c k .Hence d = 3 n , and d c = 2 n .The projector on the computational subspace Π c = ⊗ n k=1 Π c k is a tensor product where Π c k is the projector onto the computational subspace on the k-th qubit.Note that the projector onto the leakage subspace on the k-th qubit is Π l k = |2⟩⟨2| and the projector onto the leakage subspace Π l := I − Π c ̸ = ⊗ n k=1 Π l k , where I is the identity operator on H.
, and H l = i̸ =c n H i .For each Hilbert space H i , denote Π i := Π i / dim(H i ) the trace-normalized projector associated to the projector Π i .
We assume the noise of interest to be Markovian and time-independent throughout this paper.Given an ideal unitary U ∈ U(d c ), we denote U(•) := (Π l ⊕ U ) • (Π l ⊕ U † ) as the corresponding ideal unitary channel acting on the whole space.Given a completely-positive trace-preserving (CPTP) channel Û characterizing the noisy implementation of U, we further denote Λ := U † • Û as the noise information of U accounting leakage.Note that Û = U • Λ as U is a unitary channel.The average leakage and seepage rates of a channel Λ are defined as [23] L ave (Λ) = Tr Π l Λ( Π c ) ; (1) We often write L ave and S ave when the noise channel Λ being referred to is unambiguous.Unless explicitly specified, we use the term "leakage noise" to represent both leakage and seepage errors.The Pauli group with phase P < U(2) is defined as ± {1, i} × {I, X, Y, Z}, where X, Y, Z are 2 × 2 Pauli-X/Y/Z matrices respectively.Let P n := P ×n < U(2) ×n .For an element P = i P i ∈ P n , its corresponding ideal unitary channel in the full space is defined as P := i P i .For sake of simplicity, we identify the element P with its corresponding ideal channel P, and use P as a shorthand for the corresponding noisy implementation P • Λ.
Inspired by the Pauli-transfer matrix (PTM) representation [32], here we define the condensed-operator representation |•)) of linear operators as the Liouville representation [33] with respect to the orthonormal operator basis I = {Π i / dim(H i )} i∈{c,l} n .The basis is not complete in the sense that it does not span L(H); for a linear operator ρ not lying in the span of I, |ρ)) is understood as the projection of ρ onto the span of I followed by the vectorization, that is, |ρ)) := | P(ρ))), where P : L(H) → span(I); P(ρ) := i Tr (Π i ρ) Π i is the twirling projector from L(H) to span(I).For sake of clarity, in the following, we represent the condensed operator representations under the basis {| Π i ))} i , and the adjoints under the basis {((Π i |} i .Note that Π i | Π j = δ ij .Under such basis choice, for a generic linear operator A ∈ L(H), we have For a superoperator Λ, the corresponding condensed operator representation is then Since I does not form a complete basis, compositions of condensed operator representations do not directly translate to compositions of the corresponding linear operators; rather they translate to compositions of the twirled versions of the corresponding linear operators through the twirling projector P.More specifically, we have We denote [n] := {1, . . ., n} throughout the paper.

Leakage randomized benchmarking protocol
Here we present a leakage randomized benchmarking protocol that does not require actions on the leakage subspace or assumptions about SPAM errors.Our protocol is based on the assumption that the noise, represented by the operator Λ, is Markovian, time-independent, and gate independent.We further assume we have access to a noisy measurement operator Πc close to the projector to the computational subspace Π c .
(1) Given a sequence length m, sample a sequence of m Paulis P 1 , . . ., P m from P n uniformly i.i.d., and perform them sequentially to a fixed (noisy) initial state ρ0 , obtaining Pm The average leakage rate L ave and seepage rate S ave are estimated with the fitted exponents λ i .The number of exponents for p Πc (m) depends on the specific noise model of Λ.In the following, we will show the explicit representation of E [p Πc (m)] and λ i .
The Pauli group P n can twirl any quantum state in computational subspace to the maximum mixed state [34,35], i.e., 1 |Pn| Pc∈Pn P c (ρ c ) = Tr (ρ c ) Π c , where ρ c ∈ L(H c ) is a quantum state in computational subspace.Here we expand the twirling of a Pauli group from computational subspace to the entire Hilbert space, as shown in Lemma 1.
Lemma 1.Let P be the twirling projector such that P(ρ) = i Tr (Π i ρ) Π i for any quantums state ρ.Then it can be equivalently represented as the expectation of all the Pauli channels, Lemma 1 can be obtained from the twirling properties of Pauli group P n in the computational subspace.We postpone the proof of Lemma 1 into Appendix A. With Lemma 1, we can construct the connections of L ave , S ave and the multi-exponential decay curve p Πc (m), as shown in the following theorem.
Theorem 1.Given Pauli group P n with gate-independent leakage error channel Λ, the average output probability in LRB protocol p Πc (m) = (( Πc |Q m−1 |ρ 0 )), where Q := Q Λ is the condensed-operator representation of Λ and ρ0 is some noisy state determined by the input state ρ0 .The average leakage rate equals L ave = 1 − Q c n ,c n and the average seepage rate equals S ave = 1 Proof.Let P 1 , . . ., P m be the ideal gate elements sampled from P n .Then the expectation of the probability for measuring computational basis equals where ρ0 := Λ (ρ 0 ) , ρ0 is the input state with state preparation noise.Eq. ( 10) holds by Lemma 1; Eqs. ( 13) and ( 14) follows from Eqs. ( 6) and ( 4) respectively.
By the definition of Q, we have Q i,j = Tr Π i Λ( Π j ) .Moreover, for every j it holds that i Q i,j = Tr Λ( Π j ) = Tr Π j = 1 since Λ preserves the trace.This indicates that Q is a Markov chain transition matrix.By the definitions of L ave and S ave in Eq. ( 2), we have and Theorem 1 demonstrates that Pauli-twirled quantum channels with leakage can be represented as Markov chains operating on distinct leakage subspaces, including the computational subspace itself.The leakage properties can be inferred from the spectral characteristics of the transition matrix, akin to analyses of RB protocols in the computational space [36].However, this framework does not directly provide an easily applicable LRB scheme, as the transition matrix Q typically has a dimension of 2 n , resulting in complex matrix exponential decay behavior as the number of qubits increases.
Nonetheless, estimating the leakage rate can be significantly simplified in scenarios where the number of qubits is small enough to allow manageable matrix exponential decay or when additional assumptions can be made about the leakage behavior.In the subsequent sections, we propose several physically relevant leakage noise models with straightforward theoretical exponential decay curves suitable for experimental implementation.

Average leakage rate for specific noise
In this section, we present two specific leakage noise models -single-site leakage damping noise and cross-talk-free leakage noise.We also provide the respective average leakage rates for each model.
In the following, we investigate the average leakage rate for specific leakage noise where leakage only happens on a single site (qubit).For any 1 ≤ i ≤ n, we define such that {|k⟩} k∈Bi forms a basis of the specific leakage subspace H c i−1 lc n−i where only the qubit i is leaking.Let H l,(1) := i H c i−1 lc n−i be the leakage subspace that exactly one qubit is leaking, with the corresponding basis set B := i B i .We propose a single-site leakage damping noise model as a generalization to the amplitude damping noise [30]: Define the Kraus operators where probabilities The single-site leakage damping noise model is defined as a CPTP map Λ such that for any input state ρ.Denote the average leak and seep probabilities associated with the i-th site as respectively.
In the above definition, the parameters p kk ′ can be understood as the probability of the state |k⟩ flipped to |k ′ ⟩ after the leakage damping noise, and Π H\Hc∪H l,(1) in E 0 denotes that the noise model has no effect on the Hilbert space with leakage happens on more than one site.It is easy to check that i E † i E i = I, hence Λ is a CPTP map [30] in Hilbert space H. Additionally, we introduce Eq. ( 21) to simplify the representation, and we will find that the average leakage and seepage rates are only related to p i and q i for all of i ∈ [n].The prefactor 1/2 n is added to fit the definition of "average" leakage and seepage rates in Eq. (2).

Single-site leakage noise
For the particular noise model described in Definition 1, we can simplify the average leakage rate from Theorem 1 as stated in the following theorem.
Theorem 2. Let Λ be a single-site leakage damping channel as described in Definition 1.Let p i and q i be as defined in Eq. ( 21), and assume that p i > 0 for all i and q 1 ≥ • • • ≥ q n .Then after performing n-site LRB protocol, the expectation of the probability for measuring computational basis p Πc (m) = The average leakage and seepage rates of Λ are L ave = i p i and S ave = 2 n 3 n −2 n i q i respectively.
Single-site leakage model described as a Markov chain.Self-loops are omitted.Here 0 denotes computational subspace H c , i where 1 ≤ i ≤ n denotes the subspace where only the ith qubit is leaked, i.e., H c i−1 lc n−i .
Proof.If the noise model is described by Definition 1, the corresponding condensed-operator representation only acts non-trivially on the n + 1-dimensional subspace spanned by where p i , q i are defined in Eq. ( 21).This transition matrix can be illustrated in Fig. 1.Eq. ( 22) holds since , and similarly we can get other elements of Q.
Although the spectrum of the transition matrix Q cannot be explicitly solved in the general case, it is possible to derive bounds on all its eigenvalues by examining its characteristic polynomial.For simplicity, we prove the theorem under a generic scenario where n ≥ 2 and q 1 > q 2 > • • • > q n .In this case, it can be demonstrated that all eigenvalues of Q are distinct, making Q inherently diagonalizable.A detailed analysis of situations where algebraic multiplicities arise can be found in Appendix B. Denote where then the roots of function f (x) are meanwhile the eigenvalues of Q.Note that As It can be seen that f (x i ) and f (x i+1 ) always have different signs, indicating a zero in (x i , x i+1 ) for all i ∈ [n − 1].As deg(f ) = n, there is only one zero left to be determined, which is guaranteed to be real since all the other zeros are real.Let When x < x 1 , h(x) and f (x) have the same sign, and .
Therefore we have To summarize, we have a complete characterization of all eigenvalues By Theorem 1, the average leakage and seepage rates for the Pauli group with this specific noise equal i p i and 2 n−1 3 n −2 n i 2q i respectively.We assume in Theorem 2 that p i > 0 for all i.When p i = 0 for some i, the matrix Q might not be fully diagonalizable, requiring more complex data processing schemes.From a physical perspective, such complications can be mitigated by preparing the initial state such that the initial leakage on qubit i is negligible.Theorem 2 shows that when the seepage probability of all qubits are close to each other and close to leakage probability, i.e., p i ≈ p j ≈ pP and q i ≈ q j ≈ qP for all of i, j ∈ [n], then the multi-exponential decay will approximately collapse to two-exponential decay with λ 1 ≈ 1−2q P , λ 0 = 1−2q n − i p i ≈ 1−2q P −np P .With the properties of the eigenstates for eigendecomposition of the transition matrix Q, we can further simplify the exponential curve to a single decay since the coefficient of λ 1 equals zero when the state preparation noise is negligible.The leakage and seepage of the n-qubit system can be consequently derived according to L ave = i p i ≈ np P , and S ave = 2 n 3 n −2 n j q j ≈ n2 n 3 n −2 n pP , as shown in the following corollary.
Corollary 1.Let the leakage noise Λ be as described in Definition 1 such that p i = p j = pP > 0 and q i = q j = qP for different i, j ∈ [n], and assume state preparation is noiseless, then after performing n-site LRB protocol, the expectation of the probability for measuring computational basis p Πc (m m , where A, B are some real constants.The average leakage and seepage rates of Λ are L ave = i p i ≈ np P , and The decay rate 1 − 2q P − np P obtained from the LRB experiment does not provide sufficient information to fully determine L ave and S ave .Rather, additional prior knowledge is required, such as the ratio of the leakage and the seepage rates.We postpone the proof of this corollary into Appendix C.

Cross-talk-free leakage noise
Previous studies have indicated that cross-talk in real devices can be significantly minimized [19].In the subsequent subsection, we demonstrate that the exponential decay can be simplified under the condition that leakage noise occurs independently and locally across different qubits.We make the assumption that the local noise adheres to Definition 1 for each individual gate.It is important to note that in this context, the noise is inherently single-site, as each qubit possesses only one leakage site.
Corollary 2. By performing the LRB circuit in n-site cross-talk free system for the Pauli group, the expectation of the output probability for the computational subspace of the k-th qubit is equal to p Πc k (m) = A + Bλ m k , and the average leakage and seepage rates where λ k = 1 − p k − 2q k , and A, B are some real numbers, p k , q k are leakage rates associated with Eq. ( 21) in the k-th qubit.
We postpone the proof of the corollary into Appendix D. This corollary can be obtained by restricting the noise in Theorem 2 to be the tensor product form of each local noise on a single qubit.Then if p k ≈ q k or we know the relationship between p k and q k with the analysis of the system, we can estimate L, S by fitting p k from p Πc k for all of k ∈ [n] independently.We note that the cross-talk-free noise is different from the noise defined in Definition 1, since only a single qubit can leak in Definition 1.By Corollary 2, the fitted curve associated with p Πc will not follow a single exponential decay, since We can check that when n equals to 2, there will be 3 exponents,

Interleaved LRB protocol for specific target gates
In this section, we focus on benchmarking specific target gates.Benchmarking the leakage rate of an arbitrary target gate T differs from benchmarking the leakage rate of the Pauli group, as the target gate does not readily form the Pauli group.We propose an interleaved variant of the leakage for the previous interleaved randomized benchmarking protocol [10], named iLRB (interleaved leakage randomized benchmarking).We note that the target gate channel T can be any gate scheme, provided that the associated leakage noise model conforms to the form discussed in this section.The iLRB protocol is outlined as follows: (1) Sample a sequence of m Paulis P 1 , . . ., P m from P n and perform them sequentially to the noisy initial state ρ0 interleaved by target gate T to get Pm . Measure the output states and estimate (2) Repeat Step (1) multiple times to estimate p Πc (m), the expectation of p Πc (P 1 , . . ., P m ) under random choices of P 1 , . . ., P m .
When the leakage noise of the Pauli gates is negligible compared to that of the target gate T , we can benchmark any target gate T = T • Λ T where Λ T has the same leakage noise as in Definition 1, by only performing the first two steps of the above iLRB protocol.In this case, we can directly leverage Theorem 2 to get the average leakage rate of Λ T .
When the leakage noise of the Pauli gates is not negligible, however, steps (3) and ( 4) are needed to separate the target gate leakage from the Pauli gate leakage, and more assumptions on the target gate leakage noise are needed.We assume a specific case of the noise model in Definition 1, where the target gate T has the noisy implementation T = T • Λ T = Λ T • T and the noise Λ T is defined in Definition 2 with the same value for all of p i , q j in Eq. (21).
Similarly, we assume that P = P • Λ P with noise channel Λ P as defined in Definition 2 also having same value for all of p i , q j in Eq. ( 21).
Definition 2. We define the simplified single-site leakage damping noise model as a CPTP map Λ such that for any input sate ρ, where where Definition 2 is to be regarded as a particular case of Definition 1, with at most a single leak happening between each H c i−1 lc n−i ⊆ H l,(1) and H c .Such a simplified noise model has important applications such as measuring leakage for two-qubit gates on superconducting quantum chips.See more details in the next section.The requirement of the noise model for iLRB protocol can be further relaxed to more than a single leak between each H c i−1 lc n−i and H c with the same leak probability.
With the assumption of the above noise model, the average leakage rate of target gate T can be estimated with exponential decay curves of p Πc (m), p Πc,P (m) obtained from iLRB protocol, as shown in the following theorem.Usually, the state preparation noise is negligible compared with gate and measurement noise, we also show that assuming state preparation is noiseless, we can further simplify the iLRB protocol to single-exponent decay curves in the following theorem.Theorem 3.For any n-site target gate T where its noisy implementation T = T • Λ T = Λ T • T , Λ T and the noise of Pauli group P n both have the formations as in Definition 2 with noise parameter p be ϵ T , pP respectively, after performing the iLRB protocol, the expectation of the output probabilities where and the average leakage and seepage rates for target gate T equal L T = nϵ T , S T = 2 n nϵ T 3 n −2 n respectively.Assuming state preparation is noiseless, we can further simplify Eqs.(34)- (35) to Proof.By Theorem 2, we have λ P1 = 1 − 2p P , and λ P2 = 1 − (n + 2)p P .Since T • Λ T = Λ T • T , and Λ P , Λ T both have formations as in Definition 2, then where for i ∈ {0, 1, . . ., n}, then we have We also provide the details for the representation of Q in Appendix E. Let λ be the eigenvalue of Q, with the representation of Q we have which implies we have eigenvalues and and λ 3 = 1.Specifically, the multiplicity of λ 1 equals n − 1.We postpone the proof of Eq. ( 46) into Appendix F. The average leakage rate for T gate can then be determined as L T = Tr (Π l Λ T (Π c /2 n )) = nϵ T , and seepage rate The single exponential decay result of this theorem for noiseless preparation noise can be obtained from Theorem 3 and the properties of the eigenstates for the eigendecomposition of the transition matrix Q.We postpone the proof into Appendix G.
By leveraging of Theorem 3, we can estimate L T , S T using the fitted λ i estimated from iLRB protocol.

Numerical results
In this section, we carry out the numerical experiments for the average leakage rate of the multi-qubit Pauli group with the LRB protocol proposed in Section 3. Our iLRB protocol proposed in Section 5 can be applied to few-qubit cases, which is experimentally important to test the leakage and seepage of quantum gates.To support this, we show average leakage rates for iSWAP/SQiSW and CZ gates with prior noise according to the Hamiltonian of superconducting quantum devices in Appendix I.

Average leakage rate for multi-qubit Pauli group
In this subsection, we numerically implement the LRB protocol introduced in Section 3 to estimate the average leakage rate of the Pauli group.We list two examples to show the robustness of our protocol.Example 1 presents a simple noise model where the amplitude damping only happens in a pair of qubits between the set B i and {0, 1} n for any i ∈ [n] in noise model 1.To show the robustness of our protocol, in Example 2 we give a more complex noise model that contains all of the amplitude dampings of the qubit pairs between the set B i and {0, 1} n for any i ∈ [n], and we additionally add the amplitude damping for qubit pairs both in the same B i or {0, 1} n .Example 1.For an n-qubit circuit, we select a specific form of the noise Λ from the noise model in Definition 1.Let f i , g i be n-trit string denoting basis from computational subspace and leakage subspace respectively, and f i := 0 . . .01 i 1 i−1 0 . . .0, g i := 0 . . .02 i 0 i−1 . . .0 when 2 ≤ i ≤ n, and f 1 := 0...011 = f 2 , g 1 := 0...02.We define the noise model as where where where p i , q i are uniformly randomly picked from [2.5 × 10 −5 , 3.75 × 10 −5 ] for i, j ∈ [n].We take the number of qubits n = 4 in the numerical experiment.
To demonstrate the SPAM robustness of the LRB protocol, we choose a specific form of noise in the state preparation and measurement processes.Assume the state preparation process has the depolarizing noise in H c and H l .The resulting initial state can be denoted as where ρ 0 = |0⟩ ⟨0|, and p c , p l are the depolarize probabilities with p c + p l ≤ 1.The measurement noise is modeled as a perfect computational basis measurement followed by independent classical probabilistic transitions on each individual site.The probability transition matrix associated with site j is denoted as where η j0 , η j1 are 0-flip-to-1, and 1-flip-to-0 probabilities respectively, η lji , η sji are i-flip-to-2 and 2-flip-to-i probabilities respectively, where i ∈ {0, 1} for the j-th qubit.We set parameters p c = p l = 0.0001, and η j0 = 0.05, η j1 = 0.1, η lj0 = η sj0 = 0.0001, η lj1 = η sj1 = 0.0005 for any j ∈ [n].The number of qubits n = 4.We depict the probabilities of measuring outcomes in the computational subspace with the circuit size of Pauli gates as in Fig. 2. Note that here we regard a Pauli gate as a series of Pauli X/Y/Z gates without interaction.The theoretical leakage rate L ave = i p i = 3.4 × 10 −5 , and seepage rate S ave = i q i = 8 × 10 −6 .We fit the exponential decay curve with the LRB protocol proposed in Section 4. By Theorem 2 and Corollary 1 we see that if p i and q j are close to each other for all of i, j ∈ [n], then the probability p Πc will be approximately collapse to a two-exponential decay with λ 1 ≈ 1 − 2p P , and λ 0 ≈ 1 − (n + 2)p P .When the state preparation noise is small, we can approximate pP via a single exponential decay pP (m) = A + Bλ m 0 for some constants A, B. We fit the experimental data to a single exponential decay curve to obtain λ0 = 0.999957±1.2×10−5 .Then the average error p = (1 − λ0 )/(n + 2) = (7.14 ± 2.00) × 10 −6 , and thus the estimated average leakage rate Lave = np = (2.86 ± 0.80) × 10 −5 and seepage rate Ŝave = n2 n p 3 n −2 n = (7.03± 1.97) × 10 −6 .The estimated results are consistent with the theoretical ones within the errors of statistics, which verify the validity of the LRB protocol.
Example 2. The SPAM noise is set the same way as in Example 1.To show the robustness of the LRB protocol, we choose noise Λ which contains all of the flips (1) between subspace H l,(1) and computational subspace H c , and (2) inside each subspace H k for all k ∈ {c, l} n .We choose the number of qubits n = 3.The noise strength p ij is picked uniformly and randomly from interval 10 −3 [1, 1 + 10 −5 ] for (i, j) in (H c , H l,(1) ) ∪ (H l,(1) , H c ) and p ij is picked uniformly and randomly from interval 10 −6 [1, 1 + 10 −5 ] for i and j both in H k and i ̸ = j, k ∈ {c, l} n .By Theorem 2, the theoretical average leakage and seepage rates are L ave = n i=1 p i = 1.51 × 10 −5 and S ave = 2 n 3 n −2 n n i=1 q i = 6.41 × 10 −6 respectively.By Corollary 1, we fit the experimental data using a single exponential decay curve and obtain λ = 0.999974 ± 1.046 × 10 −5 .Then the average error p = (1 − λ)/(n + 2) = (5.19± 2.09) × 10 −6 , and the estimated average leakage rate Lave = (1.56 ± 0.63) × 10 −5 and average seepage rate Ŝave = (6.56 ± 2.64) × 10 −6 .The numerical results validate the LRB protocol and demonstrate that the noise in the computational subspace does not affect the average leakage rate.We depict the probabilities of measuring outcomes in the computational subspace with the circuit size of Pauli gates as in Fig. 3.

Average leakage rate for specific gates
One important application of the iLRB protocol in Section 5 is measuring leakage of experimentally realized twoqubit quantum gates.Noise in real quantum gates can be very hard to characterize due to the complexity of gate schemes.For example, in the flux-tunable superconducting quantum devices, to implement a two-qubit iSWAP gate, the two qubits are brought to resonance adiabatically, left alone to evolve for some time duration, and finally detuned adiabatically back to their normal working frequencies [31].Both the adiabatic evolution and the resonant evolution might lead to leakage and seepage.If one carries out the iLRB protocol for some specific gates proposed in Section 5, one would theoretically get one decay curve that consists of multiple exponents.A general multi-exponential decay curve is hard to fit due to statistical errors in real quantum experiments.To simplify the problem, we focus on the leakage damping noise models given in Definitions 2(The explicit form of the two-qubit case is given below).It can make the data fitting and processing more manageable.These simplified noise models are supported by the Hamiltonian evolution of the target two-qubit gates.

Average leakage rate analysis
The leakage damping noise model for iSWAP/SQiSW gate is shown below.This noise model is supported by qubits' Hamiltonian evolution.See more details in Appendix I.
where S = {(02, 11), (11,02), (20,11), (11,20)}, and where ϵ ∈ [0, 1/2].This noise model contains one parameter ϵ that remained to be fitted by the iLRB experiment.Since this noise model belongs to the noise model in Def. 2, the average leakage rate of these gates can be formalized with Theorem 3. Another commonly realized two-qubit gate in flux-tunable superconducting quantum devices is the CZ gate, of which leakage damping noise model reads (See Appendix I for more details) where S = {(02, 11), (11,02), (20,11), (11,20)} and Similar to iSWAP/SQiSW gates, the noise model of the CZ gate learned from Hamiltonian evolution can be represented as noise model (56).Since usually, the noise for single-qubit gates is much lower than that of the two-qubit gates, we make the assumption that Pauli gates are noiseless.Comparing Eq. ( 55) with Eq. ( 57), one finds the leakage damping noise model for iSWAP/SQiSW gate can be treated as a special case with ϵ 1 = ϵ 2 = ϵ.Thus for the more general leakage damping noise model in Eq. ( 57), we provide the following corollary for the data analysis after carrying out the iLRB protocol.
Corollary 3.For two-qubit target gate T with noisy implementation T = T • Λ G , where Λ G has the form defined in Eq. (56), and we assume the Pauli gates are noiseless.Then by performing the iLRB protocol, the expectation of the output probability is , and the average leakage rates L = ϵ1+ϵ2 4 , and S = ϵ1+ϵ2 5 .We postpone the proof of this corollary in Appendix J.Here we only need the assumption that the noise of T gate is right hand side of T , since . By Corollary 3, we can get the average leakage rates by fitting λ 1 , λ 2 from the exponential curve p Πc .

Discussion
In this paper, we proposed a framework of leakage randomized benchmarking that addresses the limitations of previous proposals and is more versatile in its applicability to a wider range of gates.The LRB protocol is particularly suitable for multi-qubit scenarios in the presence of SPAM noise.We presented an interleaved variant of the LRB protocol (iLRB) and conducted a thorough analysis of the leakage and seepage rates under various noise models, with a focus on the leakage-damping noise model and two-qubit gates in superconducting quantum devices.We carried out numerical experiments and see a good agreement between the theoretical leakage/seepage rates and the numerical ones for multiple gates.As the iLRB protocol is sensitive only to leakage, rather than the specific logic gate in computational subspace, it can be easily extended to other two-qubit gates realized in experiments.We leave the experimental demonstration of the iLRB protocol for future work.
One major difference between LRB and RB protocols is that single-exponential decays are guaranteed under general assumptions for RB protocols if the computational space is sufficiently twirled.However, leakage subspaces are hardly affected by any gate schemes designed on purpose for the computational subspaces, causing LRB to exhibit much more complicated decay behavior.Alternatively, gates that can twirl the leakage subspace might lead to cleaner decay behavior, but would pose somewhat unrealistic assumptions on the gate implementation that might not be experiment-friendly.In our work, we choose not to pose assumptions about the gates themselves, but instead require prior knowledge of the leakage noise models.Such prior knowledge facilitates data processing and interpretation, but their validity needs to be established either experimentally or through first-principle error analysis.Although we have proposed two simplifications under which the LRB behaviors are better understood, a more case-by-case study might be needed for other physically oriented noise models.As a complement, in Appendix J we also analyze the leakage information we can gain in a more complex noise model.
We have posed several intriguing open questions for future exploration: (1) Could we apply the LRB protocol discussed in Section 3 to compute the average leakage rate for the Pauli group, considering other Markovian and gate-independent, time-independent noise types aside from leakage damping noise?
Looking ahead, is it feasible to extend our protocol to address non-Markovian noise within the entire space ?
(2) Is the requirement for the noise channel and gate operation to commute essential when benchmarking any gate?
(3) The iLRB protocol aims to evaluate the leakage rate of the iSWAP/SQiSW gate.Experimental verification is anticipated as the next step in future work.
The middle two terms vanish since U U = U U † = 0; the first term follows the twirling property of the Pauli group in the computational subspace and the last one from that H l is one-dimensional.For general n, we then have for any n-qubit quantum state ρ.We here denote by C k a quantum operation C acting on the kth qubit.

B Complete Proof of Theorem 2
Proof.We prove Theorem 2 under more general cases where n = 1 or q i = q i+1 for some i.While the eigenvalues of such matrices can be derived from the continuity of roots of polynomials with respect to the coefficients, we prove here that Q is always diagonalizable, even if algebraic multiplicities occur.The matrix Q is defined by the two vectors ⃗ p and ⃗ q.In the following, we use the notation Q(⃗ p, ⃗ q) in case ⃗ p and ⃗ q are to be explicitly specified.
• n = 1: In this case we have where x i = 1 − 2q i ; the two eigenvalues are 1 and x 1 − p 1 .
Note that x 1 − p 1 = 1 iff Q = I, and therefore Q is always diagonalizable.
• There exists i such that q i = q i+1 .We prove that Q is similar to , where p ′ i = p i + p i+1 , p ′ i+1 = 0, q ′ i = q i , q ′ i+1 = 0, and p ′ j = p j , q ′ j = q j for all j ̸ ∈ {i, i + 1}.Such similarity is given by the following transformation Q ′ = AQA −1 , where .
By repeatedly applying such similar transformations and rearranging the rows and columns, we can reduce Q to a canonical form Q * 0 0 Σ , where Σ is diagonal, and Q * = Q(⃗ p * , ⃗ q * ) where no pairs of entries in ⃗ q * collide.Q * is diagonalizable and hence so is Q.

C Proof of Corollary 1
By Theorem 1, after performing n-site LRB protocol, the expectation of the probability for measuring computational basis equals By eigen-decomposing matrix Q Λ , where Σ is the diagonal matrix contains all of the eigenvalues of Q Λ , and V is the matrix contains all of the associated eigenstates.By Theorem 2 we see that there are three different eigenvalues, (1) λ 0 = 1 − 2q − np with multiplicity one.Let the associated eigenstate be ⃗ v = (v 0 , v 1 , . . ., v n ); (2) λ 1 = 1 − 2q with multiplicity n − 1.Let the associated eigenstates be u (s) = (u n ) T , where s ∈ [n − 1]; (3) λ 2 = 1 with multiplicity one.Let the associated engenstates be ⃗ w = (w 0 , w 1 , . . ., w n ).

H Gate representations
Here we provide the matrices representation of two-qubit gates iSWAP, SQiSW, and CZ operating on the entire Hilbert space H as follows.
When implementing the iSWAP gate, the two qubits are working at the same frequency, which means the energy spectra of the two qubits are the same.In that case, by assuming the ground state energy ω where η is conventionally called anharmonicity which quantifies the difference between energy gap ω and energy gap ω . Then, the double excitation Hamiltonian can be rewritten as With this effective Hamiltonian, we can see explicitly the symmetry between states |20⟩ and |02⟩.With these intuitions, we can write the leakage damping noise model of iSWAP gate as where S = {(02, 11), (11,02), (20,11), (11,20)}, and To see how the leakage damping noise model can be obtained from the Hamiltonian Eq. ( 127), within the double excitation subspace, we assume the initial density matrix can be parameterized as Here we assume the coefficients for |02⟩ and |20⟩ are the same, due to the symmetric structure in the Hamiltonian (127).Additionally, we assume that the evolution of the initial state follows the Schrödinger equation.Here we use ( * ) to denote some irrelevant off-diagonal entries in the density matrix.If we focus on the diagonal entries of the resulting density matrix, compared with Λ iSWAP (ρ 0 ), it can be realized that we can identify ε = ϵ kk ′ , ∀(k, k ′ ) ∈ S. Considering the definition of the average leakage and seepage rate in Eq. ( 2), if all the above approximations and assumptions hold, the leakage damping noise model and the real quantum gate would have the same average leakage and seepage rate, which verifies the reasonability of the leakage damping noise model Eq.(54) for iSWAP gate in the main text.The SQiSW gate is similar to the iSWAP gate, with an only difference at the evolution time t in Eq. (132).Thus the leakage damping noise model of SQiSW gate can also be described in Eq. ( 54), where the free parameter ϵ is different from the one for iSWAP gate.
To quantify the magnitude of ε, notice that for iSWAP gate, the evolution time t = 2π/g [31], so ε in Eq. ( 132) is determined by the anharmonicity η and coupling strength g of the two-qubit system.Taking experimental data from, e.g., Ref. [18], where η = −2π × 1.87GHz and g = 2π × 11.2MHz, we have εiSWAP ∼ 2.8 × 10 −4 .We will use this magnitude of ϵ in our iLRB numerical experiments for the iSWAP gate.
The leakage damping noise model of the CZ gate can be obtained by generalizing that of the iSWAP gate.When implementing the CZ gate in superconducting quantum devices, the effective Hamiltonian is still in Eq. (125) under RWA, thus the corresponding leakage damping noise model involves states |11⟩ , |20⟩ , |02⟩.Different from the iSWAP gate, here, we do not take the two qubits to resonance.To realize a CZ gate, by tuning the eigenenergy of one of the qubits, one would bring the state |11⟩ to resonate with, e.g.|20⟩, to accumulate phase on |11⟩ at the CZ resonance point (See figure 5), which means the energies of |11⟩ and |20⟩ are very close during the tuning and at the resonance point.Thus state |11⟩ has larger probability of leaking to |20⟩, compared with leaking to state |02⟩ [37].Thus, different from iSWAP gate, if we still write the leakage damping noise model of CZ gate as the form in Eq. ( 128), we have ϵ 20,11 ̸ = ϵ 02,11 and ϵ 11,20 ̸ = ϵ 11,02 .Further, the explicit Hamiltonian evolution for the CZ gate can be described by the Landau-Zener transition [37,38].It tells us that we can identify ϵ 11,02 = ϵ 02,11 and ϵ 11,20 = ϵ 20,11 .Thus the operators in the Eq. ( 128) can be parameterised as

Figure 2 :
Figure 2: The probabilities of measuring outcomes in computational subspace with circuit size m.Here the probability is estimated over 200 randomly selected circuits.The vertical axis denotes the estimation for p Πc , and the horizontal axis denotes the size of Pauli gates sampled from n-qubit Pauli group.(b) is the zoom-in figure of the red box curves of (a).

Figure 3 :
Figure 3: The probabilities of measuring outcomes in computational subspace with circuit size m.The estimation is over 200 randomly selected circuits.The vertical axis denotes the estimation for p Πc , and the horizontal axis denotes the Pauli gates sampled from n-qubit Pauli group.(b) is the zoom-in figure of the red box curves of (a).

Figure 4 :
Figure 4: The probabilities of measuring computational subspace with the number of circuit size of iLRB protocol.(b) is the zoom-in view of (a) with an error band.Here the probability is estimated over 500 randomly selected circuits.The vertical axes for (a) and (b) denote the estimation for p Πc , p Πc,P respectively, and horizontal axes denote the size of Pauli gates sampled from n-qubit Pauli group.(b) is the zoom-in figure of the red box curves of (a).

Corollary 4 .
For any two-qubit target gate T with noise model Λ in Eq. (134), with the assumption that Pauli group is noiseless, after performing the iLRB protocol, the expectation of the output probability p Πc (m) = A+B 1 λ m 1 +B 2 λ m 2 , where 5(3λ P − 2) .