Quantum circuits for exact unitary $t$-designs and applications to higher-order randomized benchmarking

A unitary $t$-design is a powerful tool in quantum information science and fundamental physics. Despite its usefulness, only approximate implementations were known for general $t$. In this paper, we provide for the first time quantum circuits that generate exact unitary $t$-designs for any $t$ on an arbitrary number of qubits. Our construction is inductive and is of practical use in small systems. We then introduce a $t$-th order generalization of randomized benchmarking ($t$-RB) as an application of exact $2t$-designs. We particularly study the $2$-RB in detail and show that it reveals self-adjointness of quantum noise, a new metric related to the feasibility of quantum error correction (QEC). We numerically demonstrate that the $2$-RB in one- and two-qubit systems is feasible, and experimentally characterize background noise of a superconducting qubit by the $2$-RB. It is shown from the experiment that interactions with adjacent qubits induce the noise that may result in an obstacle toward the realization of QEC.


I. INTRODUCTION
Randomness in quantum systems has been driving recent progress of quantum information science  as well as fundamental physics [35][36][37][38][39][40][41][42][43][44][45][46]. Theoretically, quantum randomness is often formulated by a unitary drawn uniformly at random, also known as a Haar random unitary. However, the Haar randomness is physically unfeasible in large quantum systems. From the viewpoint of applications, the unitaries that have similar properties of a Haar random unitary are of great importance since they can be used instead of the Haar one. When a random unitary has the same t-th order statistics as a Haar random unitary on average, it is called a unitary t-design. For instance, when a protocol exploits the tth power of the measurement probability after applying a Haar random unitary on any state, the protocol also works even if the Haar random unitary is replaced with a unitary t-design.
An important question about unitary t-designs is how to implement them by quantum circuits. Many implementations of unitary 2-designs, both approximate and exact ones, were proposed [53][54][55][56][57][58][59][60][61][62]. In contrast, only approximate implementations of unitary t-designs for general t were known [63][64][65][66][67]. Explicit constructions of exact unitary designs were left open except special cases [49,52,68]. Approximate ones typically suffice in applications, but exact designs are more preferable in certain protocols especially when they are used multiple times in a single run of the protocol. If this is the case, the error from each approximate implementation accumulates and eventually spoils the protocol.

arXiv:2102.12617v3 [quant-ph] 21 Sep 2021
One of such protocols is a randomized benchmarking (RB) protocol [22][23][24], a standard method for experimentally estimating quantum noise, where unitary 2-designs are used multiple times. Although the RB is experimentally-friendly and is widely used in various experimental systems, it reveals only the average gate fidelity. To obtain more information about the noise, a number of variants were proposed and experimentally implemented (see, e.g., Ref. [69] and the references therein), which are all based on 2-designs. It is highly expected that, by using higher-designs, much more information about the noise in quantum systems can be extracted. To this end, explicit constructions of exact unitary t-designs are important.
Constructing exact designs is, however, by far nontrivial. The difficulty is illustrated by a spherical tdesign, a random real unit vector analogous to a unitary t-design. The existence of exact spherical t-designs was proven in a non-constructive manner more than three decades ago [70]. Since then, more concise proofs and explicit constructions have been under intense investigation in combinatorial mathematics (see e.g., Refs. [71][72][73][74][75] and the references therein). In particular, it was only recently that constructions in general cases [76] and explicit constructions, in the sense that all the algorithms are given in a computable form [77], were proposed. Finding explicit constructions of exact unitary designs, since they are more complicated than spherical designs, is a rather non-trivial problem.
In this paper, we provide for the first time an explicit quantum circuit that generates an exact unitary t-design for any t on the arbitrary number N of qubits. More specifically, we show that an exact unitary t-design on d-dimensional Hilbert space, i.e., qudit, can be generated from those on smaller spaces, which is obtained based on the recent mathematical results by some of the authors [77]. Using this result, we provide an inductive construction of quantum circuits for exact unitary t-designs on N qubits: we first construct a unitary t-design on a single qubit and then extend it to N qubits. Unfortunately, the circuit fails to be efficient, but is still of practical use when the size of the system is small.
As an application of exact unitary designs, we introduce the t-th order RB, or the t-RB for short, that harnesses the power of exact unitary 2t-designs. The standard RB corresponds to the 1-RB. The t-RB enables us to experimentally characterize the higher order properties of quantum noises in the manner free from statepreparation and measurement (SPAM) errors. We especially investigate the 2-RB in detail and show that it reveals self-adjointness of the noise in the system. The self-adjointness is a new metric of the noise related to the feasibility of quantum error correction (QEC): small self-adjointness implies that the noise cannot be approximated by any stochastic Pauli noise. The noise on the system being stochastic Pauli is desirable both in theory and in practice. Stochastic Pauli noises are the commonly-used noise models in theoretical studies of QEC, since they are easy to numerically handle, and the properties of QEC, such as error thresholds and logical error rates, for stochastic Pauli noises are well-understood. Also, there is a practical advantage if the noise on the system is stochastic Pauli since they can be corrected simply by applying Pauli operators, making the error correcting scheme easier in general.
After numerically demonstrating the feasibility of the 2-RB, we perform the 2-RB in a superconducting system and estimate the self-adjointness of background noise, showing that the 2-RB experiments are feasible. From the experiment, we find that the interactions with adjacent qubits especially decrease the self-adjointness, which may lead to degradation of the performance of QEC with standard decoders. Hence, either improving the system or extending the noise model in theoretical studies of QEC, or both, is important for further experimental developments of quantum information technology.
This paper is organized as follows. In Sec. II, we provide a general introduction of unitary t-designs. Our main results are summarized in Sec. III for the quantumcircuit construction of exact unitary t-designs, and in Sec. IV for the t-RB protocols. A summary of the experiment of the 2-RB is provided in Sec. V. After we explain the structure of the remaining paper in Sec. VI, we provide a proof of the explicit construction in Sec. VII and the theory of the t-RB in Sec. VIII. The details of the experiment are provided in Sec. IX. We conclude our paper with summary and discussions in Sec. X. Technical statements are provided in Appendices.

II. UNITARY t-DESIGNS
Let U(d) be the unitary group of degree d < ∞. The Haar measure H on U(d) is the unique unitarily invariant measure on the unitary group, i.e., it satisfies Haar random unitary on average even when t copies of the unitary are given. To clasify this, let us define a quantum operation, i.e., a completely-positive and tracepreserving (CPTP) map, G µ t by for any quantum state on t qudits, where µ is either the Haar measure H(d) on a qudit or a uniform measure over a unitary t-design U t (d). Then, we can show that Definition 1 is equivalent to that (see e.g. [78]) This implies that, in any experiments that use t copies of a random unitary, no difference will be observed on average when a t-design is used instead of the Haar one. For instance, let us consider the probability distribution {p i (U ) := Tr[P i U ρU † ]} when a one-qudit state is measured by a given POVM {P i } i after the application of a unitary U . By setting the t-qudit state in Eq. (3) to ρ ⊗t and using Eq. (4), it follows that, for any s = 1, . . . , t, Thus, the distribution of the measurement outcomes for a Haar random unitary and that for an unitary t-design exactly coincide up to the t-th order on average. Note that this is merely an example, and Eq. (4) implies much more: a Haar random unitary cannot be differentiated from a unitary t-design even by more complicated experiments over t qudits.
The existence of an exact unitary t-design on U(d) for any t and d follows from the Carathéodoty's theorem and the fact that the dimension of the space, on which U ⊗t ⊗ U †⊗t is defined, is finite. Note however that the proof indicates only the existence of an exact unitary t-design. How to explicitly construct an exact unitary t-design has been a highly non-trivial problem.
In our construction, it is convenient to introduce a strong unitary t-design.
Clearly, a strong unitary t-design is a unitary t-design. Unlike standard unitary designs, strong unitary designs do not have a clear operational interpretation in quantum information processing, but we use it in the intermediate step of our construction.

III. MAIN RESULT 1 -QUANTUM CIRCUITS FOR EXACT UNITARY DESIGNS -
In this section, we provide explicit constructions of strong unitary t-designs for any t. In particular, a quantum circuit for a strong unitary t-design on N qubits is provided. We start with preliminaries in Subsec. III A, and provide the construction in Subsec. III B. We then comment on the circuit complexity of the construction in Subsec. III C.

A. Preliminaries
Unitary designs have been studied in terms of representation theory [49,51] since the operator U ⊗t ⊗ U †⊗t in the definition can be regarded as a representation of the unitary group. Our construction is based on representation theory, where irreducible decomposition of the operator plays an important role. A brief introduction of irreducible representations (irreps) will be provided in Section VII A. Here, we mention a couple of well-known facts that are necessary to state our main result.

B. Inductive constructions
Our main technical result is to construct a strong unitary t-design on U(d) from those on U(d 1 ) and on U(d − d 1 ).
where Z (d1) λ is the zonal spherical function for λ ∈ Λ sph (d 1 , d, t). Let R λ be a unitary defined by where C(θ λ ) = diag(cos θ is a strong unitary t-design on U(d).
Theorem 3 follows from a more general result [77] shown by some of the authors, which works not only for the unitary group but also for a broader class of compact groups. For the sake of completeness, we provide a direct proof of Theorem 3 in Sec. VII.
We then claim that where ω = exp[ 2π t+1 ] is the (t + 1)-th root of unity, is a strong unitary t-design on U(1) for any t. This is easily checked by direct calculations: where δ rs is the Kronecker delta. Hence, we have implying that W 1 is a strong unitary t-design. From Theorem 3 and W 1 , a strong unitary t-design on a qudit can be inductively constructed.
In this construction, it is important to obtain zeros for the zonal spherical functions Z (1) λ . This is computationally feasible since they are polynomials of one variable and are explicitly given (see Appendix A of Ref. [77]). Furthermore, Λ sph (1, d, t) contains only t elements. Hence, we need to solve t polynomials with one variable, which is tractable as far as t is not too large.
We now consider a strong unitary t-design on N qubits. Again using Theorem 3, we obtain the quantum circuit on (N + 1) qubits based on that on N qubits. See also Fig. 1.
Corollary 5. Let Q N be a strong unitary t-design on N qubits, and Ctrl-Q N be a set of controlled-unitaries on N + 1 qubits, defined by Representing {0, . . . , D − 1} in binary form such as Then, is a strong unitary t-design on N + 1 qubits.
Corollary 5 implies that a quantum circuit for an exact unitary t-design can be inductively constructed from a strong unitary t-design on one qubit, i.e., that on U(2). Furthermore, a strong unitary t-design on U(2) can be constructed using Corollary 4. Thus, combining Corollaries 4 and 5, we obtain a quantum circuit for an exact unitary t-design for any t and on an arbitrary number of qubits.
Note that the circuit, constructed in this way, can be explicitly decomposed into two-qubit gates. The controlled unitary Ctrl-Q N part contains up to three-qubit gates, if the circuit Q N on N qubit is already decomposed into two-qubit gates. The three-qubit gates can be easily rewritten as a series of two-qubit gates. Also, the X-rotation controlled by N qubits, R X (θ λ ), can be decomposed into a sequence of two-qubit gates of polynomial length using sufficiently many number of ancillary qubits, which is based on a classical oracle that computes the angle θ (j) λ from j (see Appendix A). In special cases, we can find a much more concise construction based on a similar technique.
FIG. 1. The quantum circuit that generates an exact unitary t-design on N + 1 qubits from those on N qubits. The unitary QN is a quantum circuit for an exact unitary t-design on N qubits. The gate X(θ λ ) is the single-qubit X-rotation controlled by the other N qubits, which corresponds to RX (θ λ ) in the main text. Note that this gate can be decomposed into a sequence of two-qubit gates using a classical oracle that provides θ (j) λ from j. The rotation angles θ λ are obtained by solving Z is the zonal spherical function for the spherical representation λ ∈ Λ sph (D, 2D, t) with D = 2 N . The number of the controlled-X rotations is |Λ sph (D, 2D, t)|. By the inductive use of this quantum circuit in terms of N , we can decompose the circuit to that consisting only of two-qubit gates. Proposition 6. Let C(4) be the Clifford group on 2 qubits. There exists a fixed two-qubit unitary U c , such that C(4)U c C(4) is an exact unitary 4-design on 2 qubits.
Analytically, we can prove that there exist unitaries V 1 and V 2 such that C(4)V 1 C(4)V 2 C(4) is an exact unitary 4-design on 2 qubits [77]. Also, an algorithm for computing the unitaries V 1 and V 2 is given. It, however, turns out from numerics that it is not necessary to apply two extra unitaries if we choose a proper unitary U c , which leads to Proposition 6. An explicit form of the unitary U c is numerically obtained and is provided in Appendix B. Note that the existence of U c is confirmed numerically, so the statement holds up to the numerical precision.
This construction is only for a 4-design on 2 qubits, but the number of unitaries in the 4-design is much smaller than that of Corollary 5. It is an open problem whether a similar construction works for higher-designs on a larger number of qubits.

C. Efficiency and comparison with a Haar unitary
To quantitatively evaluate the complexity of the quantum circuit for an exact unitary t-design obtained in Corollary 5, we provide an order estimate of the number G(N ) of two-qubit gates in the circuit. Assuming 2 N t and using the fact that |Λ sph (D, 2D, t)| = O(e π √ 2t/3 ) due to the Hardy and Ramanujan formula for the asymptotics of the number of partitions, we obtain to the leading order of N . Hence, it is necessary to use exponentially many two-qubit gates as the number of qubits increases. This inefficiency of the quantum circuit may be intrinsic since the construction is inductive. There is another source of inefficiency. In Corollary 5, it is necessary to find zeros of zonal spherical functions (see Eq. (19)) for all λ ∈ Λ sph (D, 2D, t). The zonal spherical function is given in terms of the summation of the (normalized) Schur polynomials (see Appendix A in Ref. [77]). It is unlikely that the Schur polynomials have polynomial size algebraic formulas in general [82]. Moreover, the number of variables for each Z Hence, finding zeros of zonal spherical functions is computationally intractable.
In total, the construction for an exact unitary t-design on a large number of qubits is inefficient from both quantum-circuit and classical-computation viewpoints. We, however, think that our construction and our proof technique will form a solid basis of searching more efficient constructions of exact, as well as approximate, unitary designs. We also emphasize that, despite its inefficiency, our construction is of practical use on a few-qubit system, as we seek in the following sections.
Before we conclude this section, we comment on advantages of our construction of an exact unitary t-design over a direct implementation of a Haar random unitary. A naive way of implementing a Haar random unitary by a quantum circuit consists of three steps. First, we sample a Haar random unitary as a matrix by a classical computer. A classical algorithm for this is known [83], but it is trivially inefficient since the size of the matrix is exponentially large. We then classically compute a decomposition of the unitary matrix into a sequence of two-qubit unitaries, providing a classical description of a quantum circuit for the unitary. This step is also inefficient, and the resulting quantum circuit is almost surely composed of the exponentially many number of two-qubit gates. Finally, we implement the circuit in practice.
This quantum circuit for a Haar random unitary is inefficient in terms of the number of qubits, and so, cannot be of practical use in a large system. Even in a small system, this naive implementation has a crucial difficulty that, every time a unitary is sampled, the above protocol outputs a quantum circuit with a rather different sequence of various two-qubit gates. This implies that, in each sampling, one needs to significantly modify the quantum circuit. This is in a sharp contrast to our quantum circuit for a unitary t-design based on Corollary 5 since it has a fixed structure. In each sampling, only what one needs to do is to randomly choose single-qubit gates, or more precisely elements of U(1) from W 1 (see Eq. (12)), and to plug them into the quantum circuit with a fixed structure. This will help practical implementations of the circuit in small systems.
It should be also noted that the single-qubit gates in our construction can be sampled from a discrete set, though sampling from a continuous gate set is necessary in the direct implementation of a Haar random unitary. This is another advantage of our construction.

IV. MAIN RESULT 2 -HIGHER-ORDER RANDOMIZED BENCHMARKING -
We here introduce a higher-order generalization of the standard RB that uses exact unitary 2t-designs. We call it the t-th order RB, or simply t-RB. The standard RB corresponds to 1-RB. From the higher-order RB, more information about the noise can be extracted. In particular, we show that a new characterization of the noise, which we call self-adjointness, can be estimated from the 2-RB.
Before we proceed, we emphasize that exact unitary designs, not approximate ones, are of crucial importance in the RB-type protocols. This is because the protocol uses unitary designs multiple times. Hence, if each unitary has an error due to the approximation, it accumulates in the whole process and results in a large error at the end. Since the goal of the RB-type protocol is typically very high, such as benchmarking the fidelity > 95%, the error originated from the approximate designs would spoil the protocol. Hence, the use of exact unitary designs is of key importance. This point is more elaborated on in Subsec. IV D.
In Subsec. IV A, we overview a couple of metrics of the noise, i.e., the average fidelity and unitarity, and introduce the self-adjointness. The importance of the selfadjointness in QEC is argued in Subsec. IV B. We then introduce the t-RB in Subsec. IV C. We argue the importance of exact designs in more detail in Subsec. IV D. We focus on the 2-RB in Subsec. IV E and show that the self-adjointness and the unitarity of the noise can be estimated from the 2-RB at the same time. We briefly comment on the scalability of the t-RB in Subsec. IV F.

A. Characterizing noises
A noise E acting on a q-qubit system is formulated by a completely-positive and trace-preserving (CPTP) map. Let d be defined as d := 2 q . The average fidelity and the unitarity are defined by respectively, where E (ρ) := E(ρ − I/d), and ||A|| 2 = (Tr[A † A]) 1/2 is the Schatten 2-norm. The average fidelity satisfies 1/(d + 1) ≤ F (E) ≤ 1, and F (E) = 1 if and only if the system is noiseless, i.e., E is the identity channel, while the unitarity satisfies 0 < u(E) ≤ 1, and u(E) = 1 if and only if the noise is coherent, i.e., E is a unitary channel. The unitarity is an important metric in the context of QEC since coherent noise is known to be hard to correct in general [84][85][86].
In the RB-type protocols, it is more natural to use a fidelity parameter f (E) rather than the average fidelity itself. It is defined by and satisfies −1/(d 2 − 1) ≤ f (E) ≤ 1. We next introduce a self-adjointness of the noise. For any linear map E, an adjoint map E † is defined by which is equivalent to that all the Kraus operators of E are self-adjoint.
The self-adjointness H(E) of the noise E is defined by The normalization constant is chosen such that 0 ≤ H(E) ≤ 1. Obviously, H(E) = 1 if and only if E is selfadjoint, i.e. E = E † . Note that the self-adjointness has two contributions from the noisy map E, one is from the unital part and the other from the non-unital part. The non-unital part of the noise makes the self-adjointness less than one since, if E is not unital, then E † is not trace-preserving, which implies that E = E † .
To clearly separate the two contributions, we introduce a self-adjointness parameter h(E).
The self-adjointness parameter h(E) is related to the selfadjointness H(E) and the unitarity u(E) by where |α E | is a measure of the non-unital part of the noise (see Subsec. VIII A for the definition). We can clearly observe that H(E) consists of two factors, the unital part h(E) and the non-unital part |α E |. The three metrics of noises, namely, fidelity, unitarity, and self-adjointness, all capture different properties of the noises. The fidelity reveals the first order property of the noises, while the unitarity and the self-adjointness, which are independent to each other, reveal the secondorder. In order to improve noisy quantum devices, it is of crucial importance to obtain the information of noise as much as possible. Hence, it is certainly of practical use to introduce the self-adjointness as a new metric of noise. In addition, we argue in the next subsection that the self-adjointness has important implications for QEC.

B. Importance of self-adjointness in QEC
The most important family of self-adjoint noises is stochastic Pauli noises, whose Kraus operators are all proportional to Pauli matrices. In QEC, Pauli noises are the standard yet most important class of noises both in theory and in practice. From a theoretical perspective, Pauli noises are easy to numerically handle. Hence, most numerical calculations have been carried out by assuming Pauli noises, and it has been confirmed that QEC has preferable features, such as exponential decreases and threshold behaviors of logical error rates, if the noise is Pauli.
The noise being Pauli is also practically preferable in experimental realizations of QEC since it typically simplifies the decoding tasks. This is especially the case for stabilizer codes, such as surface and color codes, whose standard decoders are to estimate what types of Pauli operators should be applied on which physical qubits during recovery operations. For stochastic Pauli noises, if the estimation goes well, the state is fully retrieved with high probability by applying Pauli operators to the suitable physical qubits. In contrast, it is not possible to fully correct non-Pauli noises by applying Pauli operators since they generate undesired coherence between different code spaces. Thus, QEC of non-Pauli noises generally suffers from degradation of logical error rates when the standard decoders are used [86,87] or requires more complicated algorithms for retrieving the performance of QEC. Neither of them is preferable in practice since it induces additional experimental difficulties.
For these reasons, it is desirable to check that the noise on an experimental system is stochastic Pauli. To this end, the self-adjointness provides useful information since, if H(E) 1, then the noise is far from self-adjoint and cannot be approximated by Pauli noises. This implies that the practical situation differs from the standard assumption in theoretical studies of QEC and incurs additional difficulties on decoding procedure. Thus, the self-adjointness provides practical information about the feasibility of QEC using Pauli-based decoders.
Note that the difficulty of QEC for non-Pauli noises, captured by the self-adjointness, highly depends on the assumptions in quantum error correction schemes. When any decoding procedure is available, it would not be so important whether the noise is Pauli or non-Pauli. When this is the case, the unitarity will be a more suitable metric of noise relevant to the feasibility of QEC [84][85][86]. Note also that non-Pauli noises can be always transformed to a Pauli noise by Pauli-twirling. However, Pauli-twirling induces additional noise onto the system and, as a result, the performance of QEC will degrade. Thus, it is practically desirable to manufacture the system so that the noise is stochastic Pauli.
We also provide a pedagogical example of noise, where performance of QEC can be directly captured by the selfadjointness but not by fidelity nor unitarity. Consider a θ-rotation error around the X-axis on one qubit, i.e., exp[iθX/2], where X is the Pauli-X operator. The average fidelity F θ and the self-adjointness H θ can be obtained as The unitarity is 1 for any θ.
One may expect that the π/2-rotation error is easier to correct than the π-rotation since the former has higher fidelity than the latter. However, this is not the case since π-rotation is simply a perfect bit-flip that can be trivially corrected, while the π/2-rotation error is known to be particularly hard to correct [88]. Thus, neither the average fidelity nor the unitarity, which is 1 for both errors, is a good metric of the error correctability. In contrast, the self-adjointness clearly captures whether the error can be corrected, at least in this case, since H π/2 and H π are the minimum and the maximum values of the self-adjointness, respectively.

C. General description of the t-th order RB
We now introduce the t-RB using an exact unitary 2tdesign U 2t := {U i } i . As is the case for the standard RB, we assume that the noise is gate-and time-independent, so that the noisy implementation of U 2t is given by where E is the CPTP map that represents the noise, and we used the notation that U(ρ) := U ρU † .
Let O ini and O meas be the initial and measurement operators, respectively, which we assume to be Hermitian. We first apply a sequence of unitaries 2t . We then apply its inverse U im+1 := U −1 i , and measure O meas .
If the system is noiseless, E = id, this protocol results in a trivial expectation value that due to the inverse unitary U im+1 . However, when the system is noisy, the expectation value becomes The basic idea of the RB-type protocol is to extract some information about the noise E from the difference.
In the t-RB, we especially focus on the average of the t-th power of the expectation value over all choices of the unitary sequence. That is, (33) Using the representation-theoretic technique, it can be shown that V (t) (m, E|O ini , O meas ) is generally given in the following form: where λ labels the irreps of the unitary group,Â (t) λ and C (t) λ (E) are m λ × m λ matrices with m λ being the mul-tiplicity of the irrep λ. This is well-known in the literature of RB-type protocols, but we provide a proof in Sec. VIII C for completeness.
Despite its abstract expression, Eq. (34) has an important implication that the matrixĈ (t) λ (E) m depends only on E and m, but not on O ini and O meas . Hence, from the experimental data of V (t) (m, E|O ini , O meas ) for various m, it is in principle possible to estimate the matrix C (t) λ (E), which contains certain information of the noise E, in the way independent of O ini and O meas .
In practice, the most important situation is when the representation is multiplicity-free, i.e., m λ = 1 for any λ. In this case, V (t) reduces to a much simpler form: where is a bounded function. Hence, in this case, V (t) becomes a sum of some exponentially decreasing functions with respect to m.
To be more concrete, let us recall the standard RB, corresponding to the 1-RB. As shown in Ref. [22], V (1) is given by where A

D. Importance of exact designs in RB
In the RB protocol, it is important to use exact unitary designs because designs are used many times, sometimes a few hundreds to a thousand, in a single run of the protocol. To illustrate this, let us consider the 1-RB when the unitary 2-design in the protocol is -approximate.
Let m be the length of the unitary sequence as above. It is straightforward to show that, where E i 's are some constants that depend on O ini , E(O meas ), f , and how the design differs from the exact one. See Subsec. VIII D for the derivation. Compared to the 1-RB with exact ones, i.e., Eq. (36), fitting this function with respect to m is much harder since it is not a simple exponential decay.
The fitting may go well if f m /m. This requires a very high precision of the design since m can be a few hundreds in actual experiments. For instance, when f = 0.95, the degree of approximation of the unitary design should be order 10 −5 or so. Although it is possible to achieve this degree of approximation by a sufficiently long quantum circuit [64,65,67], the RB becomes unpractical if we use such a long circuit at every use of a unitary design in the protocol and repeat it a few hundreds times.
There might be a possibility to improve Eq. (37) by using different constructions of approximate unitary designs at every step, by which the differences from the exact design may become random so that they cancel out in total. This will be an interesting question, but at this point, it is not clear if such a technique works. Also, even if it works, we need to assume additional structures of approximate constructions.
The higher-order RB with approximate designs will incur more difficulty in practice. Since it uses higher moment of the outcomes, the fitting function becomes more complicated than Eq. (37) when one uses approximate designs. Similarly to the 1-RB with approximate 2-designs, much better degree of approximation, that is, longer quantum circuits, will be needed, which is not practical. Thus, we conclude that exact unitary designs are of crucial importance in a practical implementation of the t-RB.

E. Second-order RB
We next focus on the 2-RB using exact unitary 4-designs, and show that the 2-RB reveals the selfadjointness of the noise. To this end, we set the initial operator O ini to a traceless one, i.e., Tr[O ini ] = 0. This setting, together with the fact that the noise is trace-preserving, makes the representation multiplicityfree (see Appendix C). Hence, the expectation value V (2) (m, E|O ini , O meas ) for the 2-RB is given by a sum of exponentially decaying functions as shown in Eq. (35).
Note that the expectation value for a traceless initial operator can be obtained by performing the same experiment for two different quantum states ρ and ρ , and by taking the difference of the expectation values before they are squared. That is, Our second main result in this paper is about V (2) (m, E|∆, O meas ) as summarized in Theorem 7.
Theorem 7. In the above setting, V (2) (m, E|∆, O meas ) is given as follows. For single-qubit systems, where f (E), u(E), and h(E) are the fidelity parameter, the unitarity, and the self-adjointness parameter of the noise E, respectively. For multi-qubit systems, where See Subsec. VIII E for the proof.
In the single-qubit case, V (2) is a sum of two exponentially decaying functions with respect to m. Hence, from the double-exponential fitting of the experimental data of V (2) , we can simultaneously estimate u(E) and Since it can be shown that the former is not less than the latter, we can estimate which of the two decaying rates corresponds to which quantity without any ambiguity. It is also possible to estimate the fidelity parameter f (E) from the same data set by computing V (1) (m, E|∆, O meas ) because a unitary 2-design is also a unitary 1-design. Thus, from the experiment of the 2-RB on a single qubit, all of f (E), u(E), and h(E) can be estimated simultaneously.
In multi-qubit systems, V (2) has a little more complicated form and consists of four exponentially decaying functions. Also, the decaying rates do not directly correspond to neither the unitarity nor the self-adjointness parameter. We observe from Eq. (41) that h(E) can be obtained from a linear combination of the decaying rates C λ (E), the value of u(E), and f (E).
One may think that, in the case of multiple qubits, it is practically intractable to accurately fit four exponentially decaying functions from experimental data because each data point has an error. This difficulty can be circumvented by choosing appropriate initial and measurement operators. By doing so, we can set some of A λ zero in the ideal situation (see Tab. I). This allows us to estimate the decaying rates one by one. Note that the initial and measurement operators in Tab. I are all diagonal in the computational basis. Hence, it suffices to perform the experiments for the four initial operators |00 , |01 , |10 , and |11 , with the measurement in the computational basis. From the data of these experiments, it is possible to reproduce all cases listed in Tab. I by post-processing.   (40), for the 2-qubit case. The first row provides a pair of the initial and measurement operators (∆, Omeas).
We have assumed that the average fidelity of the noise E is close to 1, so that the inverse unitary Ui m+1 can be applied nearly noiseless (see Subsec. VIII E for the detail). The operator ZZ is Z ⊗ Z, and ρ− := |00 00| − |11 11|. By choosing proper operators, we can set some coefficients zero, so that experimental estimations of C λ (E) become easy.
In the multi-qubit case, the ambiguity remains to decide which of the decaying rates corresponds to which quantity. This is the case even when we use the above step-by-step estimation of the rates since, for instance, it is not clear if the unitarity u(E) is larger or smaller than C I . In this case, we need to additionally perform the unitarity benchmarking [69,89] to separately estimate u(E). If we have an estimated value of u(E), the step-by-step estimation allows us to decide all decaying rates without any ambiguity.
See Sec. V A and Sec. IX A for the performance of 2-RB in concrete cases.

F. Scalability
The t-RB for t ≥ 2 inherits most of the desired properties of the RB-type protocols. For instance, it is experimentally-friendly since, apart from using higherdesigns, the difference of the t-RB from the standard RB (1-RB) is only taking the t-th power of the expectation value before the average. It is also true that the t-RB is free from SPAM errors (see Eqs. (34) and (35)).
The property that the standard RB does have and the t-RB does not in general is the scarability. This is for two reasons. First, no efficient construction of exact unitary 2t-designs is known for t ≥ 2 so far. Second, in the t-RB protocol, it is necessary to apply the inverse unitary at the end of the unitary sequence. Hence, we need to beforehand compute the inverse of each sequence. When the system is large, the task is intractable in general. This difficulty is avoided in the standard RB by using the Clifford group, which is an exact unitary 2-design. Since the inverse is contained in the group, we can find the inverse relatively easily. One may expect that the difficulty of finding the inverse could be also avoided in the t-RB by using the 2t-design that is also a group, which is called a unitary 2t-group [52]. However, it is known that unitary 2t-groups do not exist for t ≥ 2 if the number of qubits ≥ 3. Thus, in the t-RB for t ≥ 2, the hardness of finding the inverse in a large system is inevitable.
Nonetheless, we emphasize that, in the current experimental situations, the RB-type protocols for more than three qubits are practically intractable due to the limitation of the coherent time. Thus, the experimental use of the RB-type protocols is currently aiming to characterize the noise on one-or two-qubit systems in a concise manner. Considering this fact, even if the t-RB is not scalable, it is practically useful and beneficial: it is as concise as the standard RB and provides more information about the noise, such as self-adjointness.

V. MAIN RESULT 3 -2-RB IN A SUPERCONDUCTING SYSTEM -
We finally implement the 2-RB in a superconducting system and estimate the self-adjointness of background noise. Unlike the analytical studies, the expectation values and the average over a unitary 4-design cannot be taken with arbitrary precision in experiments since the number of repetitions of experiment is practically limited. To check that this limitation does not cause any problem in the evaluation of the self-adjointness, we start with numerically investigating the feasibility of the 2-RB in Subsec. V A. We then provide a summary of experimental results in Subsec. V B.
In recent years, a number of experiments have been performed to characterize various noises on superconducting quantum systems in detail [90][91][92][93]. From our experiments, we show that the interactions with the adjacent qubits particularly decrease the self-adjoinenss and may cause problems toward realizations of QEC. In par-ticular, our result implies that there exists a gap between the superconducting system and the common noise model used in theoretical studies of QEC, and also that the standard decoders of stabilizer codes may suffer from degradation of logical errors. Hence, toward the realization of QEC, it is desired to further improve the system or to develop the theory of QEC.

A. Numerical evaluation
When the 2-RB is practically implemented, there are two additional concerns. One is originated from the fact that the expectation value O meas Oini,i is obtained from a limited number of measurements in the basis of O meas , resulting in an error due to a finite number of measurements. The other originates from the evaluation of the average E U i ∼U ×m

One-qubit cases
In the case of single-qubit systems, we consider a specific noisy map given by (45) which is characterized by three parameters p, q, θ. The first term of the right-hand side represents a unitary part and the second term represents a stochastic part of the noise. A parameter q determines a ratio between them. Hence, we can consider q as a coherent parameter of noise, e.g., noise is unitary when q = 1 and is a probabilistic Pauli noise when q = 0. The parameters θ and p represent the rotation angle of the unitary part and the error probability of the stochastic part, respectively. For simplicity, we choose θ such that the fidelity parameters of unitary and stochastic parts are equal, that is, p = sin 2 θ. Then, the fidelity parameter f (E 1 ) becomes independent of the coherent parameter q.
To perform the 2-RB for this noise, we may use the exact 4-design constructed in Corollary 5. However, it is known that the icosahedral group, which we denote by I, forms an exact 4-design on one qubit [49]. Since the icosahedral group has less cardinality than our inductive construction, we use it in the following analysis.
The numerical results for the 2-RB on a single qubit are shown in Fig. 2. For each sequence length m, we have taken 5000 random unitary sequences from I ×m and have had 5000 measurements to obtain a single data point of V (1) . A detailed fitting procedure is provided in Subsec. IX A.
To check the accuracy of the 2-RB, we consider the relative errors |y −ỹ|/(1 − y), where y andỹ are the theoretical value and the fitting value, respectively. Note that 1 − y ∼ 0 for all the fitting values when a fidelity close to unity is achieved. For almost all data points of F (E), u(E), and H(E), we find that the relative errors are less than 5.0%, except the case when p is large, or equivalently, when the fidelity is small. The relative error becomes moderately large, such as 35%, when p = 0.4 and q ≥ 0.1, corresponding to F (E 1 ) = 0.7. This is because the decaying rate of the second term in Eq. (39) is rather small, making the fitting difficult. However, such a case is not practically relevant since the fidelity is typically > 90%. Thus, we conclude that the 2-RB on 1-qubit systems works well in practice.
To analyze the dependence of the accuracy of the 2-RB on the number of measurements and samplings of random unitary sequences, we additionally perform the 2-RB on one qubit with the various numbers of measurements and samplings. The results are summarized in Tab. II, where we set the noise parameters to p = 0.02 and q = 0.02. From these results, it appears that setting the numbers of measurements and samplings of random sequences to a few hundreds is sufficient for a good estimate. These results further indicate that increasing the number of random sequences rather than the number of

Two-qubit cases
For two-qubit systems, we consider the noise given by which is similar to the one-qubit case. We choose θ as p = sin 2 θ, so that f (E 2 ) is independent of the coherent parameter q. In this case, we use the construction of exact unitary 4-designs given in Proposition 6.
In the two-qubit case, it is needed to fit the experimental data by a sum of four exponentially decaying functions, which is in general not easy especially when each data point has errors caused by the finite number of measurements and samplings of unitary sequences. To avoid this difficulty, we use the method explained in Subsec. IV E, and determine the four decaying rates, i.e., u(E 2 ), C I (E 2 ), C II (E 2 ), and C III (E 2 ) in Eq. (40), one by one.
The results are shown in Fig. 3. We have taken 10 4 random unitary sequences for each sequence length m and the parameters p, q.   scenario. Thus, the 2-RB works in actual situations also in the case of two-qubit systems.

B. Experimental implementations of the 2-RB
We demonstrate the 2-RB in a superconducting-qubit system. We first explain the setup of our experiments, and then verify the feasibility of the 2-RB experiment by comparing the unitarity obtained from the 2-RB with that from the unitarity benchmarking (UB) [69,89]. We finally characterize background noise of the system. As the background noise is gate-and time-independent, it satisfies the assumptions of the 2-RB (see Subsec. IX B for the detail).

Experimental setup
We use two superconducting qubits (Q 1 and Q 2 ) coupled with each other via an electric dipole interaction, which are a part of our 16-qubit device [94]. In all the experiments below, we use the qubit Q 1 as a target qubit of the single-qubit 2-RB and, in some experiments, Q 2 as an environmental qubit that induces additional error onto Q 1 .
The simplified system Hamiltonian H is formulated as follows, where ω i /2π is the eigenfrequency of the i-th qubit and χ ge /2π = −0.760 MHz is an effective interaction strength between the qubits [95]. It can be interpreted that the eigenfrequency of Q 1 switches depending on the quantum state of Q 2 . When Q 2 is in the |0 (|1 ) state, Q 1 has the eigenfrequency (ω 1 + χ ge )/2π ((ω 1 − χ ge )/2π). In the Bloch sphere representation, the state vector of the qubit rotates around the Z-axis with its eigenfrequency as the angular velocity. We use a local oscillator synchronized with the eigenfrequency of the qubit for observation. The state vector is stationary in a rotating frame of the local oscillator since the Z-axis rotation speed of the Bloch vector matches with that of the measurement basis. The rotation frame picture also holds when the qubit Q 1 couples to the adjacent qubit Q 2 when the qubit Q 2 is in the |0 or |1 state. For instance, when the qubit Q 2 is always in the |0 state, the eigenfrequency of Q 1 is (ω 1 + χ ge )/2π. We can detune the frequency of the local oscillator from the qubit frequency ω 1 by χ ge to make the state vector of Q 1 stationary.
It is, however, impossible to keep track of the eigenfrequency of the qubit when the state of the adjacent qubit varies. This results in an inevitable Z-rotation occurring in the quantum state. In an actual experiment involving multiple qubits, the frequency of the local oscillator is usually set to ω 1 /2π to minimize the average Z-rotation angle. See Subsec. IX B for the detail.

Comparison with the UB
In the experiment aiming to compare the 2-RB and the UB on a single qubit, we use only Q 1 and add an artificial noise after applying each gate. The isolation of the qubit Q 1 from the qubit Q 2 can be done by keeping the qubit Q 2 in the state |0 and by setting the frequency of the local oscillator to (ω 1 + χ ge )/2π, which effectively cancel the coupling between Q 1 and Q 2 . About the noise, we especially choose a single-qubit Z-rotation by angle 0.2π, denoted by R Z (0.2π). Both in the case of the 2-RB and the UB, we use the icosahedral group I and the Clifford group on a single qubits, respectively. Note that the former is an exact 4-design on a single qubit, and the latter is an exact 2design.
We have taken 100 and 1000 random sequences for the 2-RB and the UB, respectively. This is because the UB with the Clifford group converges slower than the 2-RB with the icosahedral group, which is likely due to the fact that the former and the latter are based on unitary 2-and 4-designs, respectively. A higher-design typically leads to a quick convergence since it is more concentrating around the average [96]. A faster convergence of the UB with 4-design is expected, which highlights the potential use of a higher-design also for the UB. We have taken 10 4 measurements for each random sequence to obtain a data point of V (2) . The results are summarized in Tab. III.
From the results, we observe that the unitarity characterized by the 2-RB matches with that by the UB. This indicates that the 2-RB on our single-qubit system works to characterize the gate performance.
Note that the difference between the unitarity from the 2-RB and that from the UB is slightly beyond the standard deviation. This is likely because the noise property varies in the UB experiment. As mentioned, we have taken 1000 random sequences in the UB to ensure the convergence of the statistical average, which has taken more than 10 hours in total. Since the noise in the experimental system drifts in such a long timescale, the situation of the experiment deviates from the ideal situation, where time-independence of the noise is assumed. Indeed, unlike the theoretical prediction of the UB, the data is slightly different from a single-exponential decay. This deviation is expected to be the origin a less precise value of the unitarity estimated from the UB.

Characterizing background noise
We next perform the single-qubit 2-RB, aiming to characterize background noise of the qubit Q 1 in the experimental system. We intentionally insert a delay time t after each application of a gate to extract the information of the background noise.
In the following experiments, we have taken 100 random unitary sequences from I ×m , where I is the icosahedral group, for each sequence length m and have had 10 4 measurements for each random sequence to get a data point of V (2) .
In the first experiment, we set the frequency of the local oscillator to (ω 1 + χ ge )/2π and treat the qubit Q 1 as a target qubit isolated from the qubit Q 2 . The background noise of the isolated qubit is often phenomenologically modeled by the Lindblad Master equation given by where L 1 =â/ √ T 1 represents the energy dissipation with the relaxation time T 1 ,â = (X + iY )/2 is an annihilation operator of the qubit, and L 2 = Z/ 2T φ represents the phase dissipation with the relaxation time T φ = 1/(1/T 2 − 1/2T 1 ). By solving the Eq. (48), we can obtain phenomenological predictions about the background noise E t corresponding to the delay time t.
We sweep the delay time t from 100 to 500 ns. The value V (2) obtained from the experiments is shown in Fig. 4 (a). We estimate the unitarity and the selfadjointness from V (2) through a fitting based on a sum of two exponentially-decaying curves given in Eq. (39). However, we observe single-exponential decays from the results. This indicates two possibilities. One is that the two decaying rates are nearly the same. The other is that one of the two decaying rates is much smaller than the other, so that one exponentially-decaying curve becomes quickly negligible as m increases.
In our experiment, the former is the case because the average fidelity is high, which is confirmed from the 1-RB. We can analytically show that the two decaying rates typically coincide when the fidelity is sufficiently high. More specifically, we have (see Eq. (101) in Subsec.VIII B) to the first order of , where is the infidelity 1 − F (E). This implies that the two decaying rates in Theorem 7 are approximately greater than 1−4 and 1−6 , respectively. Thus, if 1, which is indeed the case in our system, the two decaying rates are hard to distinguish, making the curve of V (2) a single-exponential decay.
We, hence, estimate the single-exponential decay rate from the experimental data of the 2-RB and derive u(E t ) and h(E t ) from Here, u(E t ) is obtained from the estimated decaying rate, and f (E t ) from the 1-RB (See Eq. (39)). The obtained fidelity, unitarity, and self-adjointness are summarized in Fig. 4 (b).
In calculating selfadjointness, we solved Eq. (28), where we substituted α E of the phenomenological prediction. They reveal that the background noise of the isolated qubit has the unitarity u(E t ) that slowly decreases as the delay increases, while its self-adjoingness H(E t ) is nearly independent of the delay. As we have explained in Subsec. IV A, a problem may occur when the unitarity is high and the self-adjointness is low, which is not observed in this experiment. Hence, we conclude that the background noise in this case would not cause any problem toward the realization of QEC.
In the second experiment, we set the frequency of the local oscillator to ω 1 /2π and treated Q 1 as a target qubit exposed to the noise induced by the adjacent qubit Q 2 . In this experiment, no control pulses are applied to Q 2 , so that Q 2 is expected to remain in the |0 state. This leads to a continuous rotation of the state vector of Q 1 by the interaction Hamiltonian term of χ ge Z/2.
In this case, the background noise with the interaction Hamiltonian is modeled by the Lindblad Master equation written as follows, providing a phenomenological model. Similarly to the first experiment, we sweep the delay time t from 60 ns to 180 ns. The delay time is set to a shorter time than the first experiment because the fidelity deteriorates due to the Z-rotation error. Since the Z-rotation error does not affect the unitarity, we conclude that the decay rate, which is less sensitive to the delay time than the other, corresponds to the unitarity. The results of the experiment are shown in Fig. 5 (a). As seen from the results, the curves V (2) obey double-exponential decay. From the two decaying rates, we obtain the unitarity u(E t ) and the self-adjointness H(E t ) as a function of the delay time as depicted in Fig. 5 (b). Note that, although the unitarity may seem different from the former experiment, it is merely due to the different time scale of the horizontal axis. The unitarities in the two experiments indeed coincide within the standard deviation (see, e.g, the delay time 100 (ns)).
The experimental results qualitatively coincide with the phenomenological predictions obtained from Eq. (51) (see Fig. 5 (b)). However, the experimental values tend to be smaller. This indicates that there exist noise sources not included in the phenomenological model. The candidates of the additional noise sources are calibration errors in the R X (π/2) gates, the initial thermal excitation rate of Q 2 (7.2 %), and the interaction of Q 1 with adjacent qubits other than Q 2 . Note that the initial thermal exci- Compared to the first experiment (Fig. 4), we observe from Fig. 5 that the fidelity F (E t ) and the self-adjointness H(E t ) quickly decrease as the delay time t increases. The latter decreases especially quickly: H(E t ) ≈ 0.48 at t = 140 (ns). This implies that, even if the fidelity is moderately high (F (E t ) ≈ 0.91 at t = 140 (ns)), the extra Z-rotation induced by the interaction with another qubit radically changes the property of the noise and makes the noise far from self-adjoint. Consequently, as the delay increases, the noise quickly becomes the one that cannot be approximated by any stochastic Pauli noise.
This result has an important implication toward a realization of QEC. As mentioned in Subsec. IV A, theoretical studies of QEC commonly assume stochastic Pauli noises to numerically compute error thresholds and error rates. Our result implies that, when the interaction with another qubit is non-negligible, we cannot directly apply the theoretical predictions based on Pauli noises. This problem will be more prominent when the system size grows since, in a large system, a qubit interacts with more qubits in an uncontrolled manner, making the noise much less self-adjoint and much far from Pauli noises. To circumvent this, effective cancellation of the dipole interaction is of great importance in the further improvement since the dominant interaction between qubits should be originated from the electric dipole interaction.
This feature of the noise, i.e., interactions with other qubits induce small self-adjointness and the difficulty of approximating the noise by a Pauli noise, is expected to be common in any experimental systems. The 2-RB experiment and the self-adjointness offer a useful method and measure, respectively, to experimentally evaluate the noise in the system from this perspective.

VI. STRUCTURE OF THE REMAINING PAPER
The remaining of this paper is organized as follows. In Sec. VII, a proof of Theorem 3 is provided. A brief introduction of representations of the unitary group is also provided before the proof. We then explain the higherorder RB in Sec. VIII, including the proof of Theorem 7. The methods used in the numerical analysis, and the experimental demonstrations are provided in Sec. IX. After we summarize the paper in Sec. X, we prove technical statements in Appendices.

VII. CONSTRUCTIONS OF EXACT DESIGNS
In this section, we provide a proof of Theorem 3. We start with a brief introduction of representations of the unitary group in Subsec. VII A, and prove Theorem 3 in Subsec. VII B.

A. Unitary t-designs and representation theory
Unitary t-designs are closely related to representations of the unitary group since the operator U ⊗t ⊗ U †⊗t in the definition can be regarded as a representation ρ of U ∈ U(d) on H ⊗2t d with H d being the Hilbert space with dimension d, i.e., ρ(U ) = U ⊗t ⊗ U †⊗t . It is natural to consider irreps of the unitary group.
A well-known fact is that each irrep can be indexed by a non-increasing integer sequence λ := (λ 1 , λ 2 , . . . , λ d ), i.e. λ 1 ≥ λ 2 ≥ · · · ≥ λ d , of length d. In particular, each irrep in U ⊗t ⊗ U †⊗t can be indexed by an element of a set Λ(d, t) defined by where λ + and λ − are the absolute value of sum of positive and negative λ i 's, respectively. Using this notation, the representation space H ⊗2t where m λ is the multiplicity of the irrep λ. Accordingly, the map ρ is also decomposed into the irreducible ones ρ λ . Based on the irrep (ρ λ , V λ ) of the unitary group, a unitary t-design U t (d) can be characterized in a representation-theoretic manner: for any λ ∈ Λ(d, t), The strong unitary t-designs are similarly characterized in terms of irreps [49]. To this end, let Λ ≤ (d, t) be where λ + is not necessarily equal to λ − . Then, a strong unitary t-design U ≤t (d) satisfies for any λ ∈ Λ ≤ (d, t).
One of the merits in this characterization is that the right-hand-sides of Eqs. (54) and (56) are zero for all nontrivial irreps due to the Schur's orthogonality relation, which states that, for any unitarily inequivalent irreps λ and λ , for any i, j, i , j , where (ρ λ (U )) ij is the (i, j) element of the matrix. By setting the irrep ρ λ to a trivial irrep, i.e., ρ λ (U ) = 1 for any U ∈ U(d), we have for any non-trivial irrep λ. On the other hand, for any trivial irrep λ, it is trivial that From these facts, (strong) unitary t-designs can be defined in terms of representation as follows: Definition 8 (Unitary designs in representation theory). An ensemble U t (d) of unitaries is an exact unitary tdesign if it holds for any irrep ρ λ with λ ∈ Λ(d, t) that An ensemble U ≤t (d) is a strong unitary t-design if Eq. (60) holds for any irrep ρ λ with λ ∈ Λ ≤ (d, t).

B. Proof of Theorem 3
We now prove Theorem 3, which states that W d defined by is a strong unitary t-design on U(d). Here, where U ≤t (d) and U ≤t (d−d 1 ) are strong unitary t-designs on U(d) and U(d−d 1 ), respectively, and R λ is constructed by solving the zonal spherical function Z λ . It suffices to show for all non-trivial irreps indexed by λ ∈ Λ ≤ (d, t). Note that the average over U ∼ W d consists of the independent averages over all W d1⊕d−d1 , further consisting of those over the strong unitary t-designs U ≤t (d 1 ) and U ≤t (d−d 1 ). Let us first fix a non-trivial irrep λ ∈ Λ ≤ (d, t) and consider W λ defined by Since we consider only irreps λ ∈ Λ ≤ (d, t), this average can be replaced with the averages over the prod- To investigate W λ , we consider the irreps of K. Since K is a subgroup of U(d), each irreducible space V λ of U(d) is decomposed into a direct sum of those of irreps of K. For the same reason as in Definition 8, every nontrivial irrep of K becomes zero by taking the average over d 1 ). Hence, if the non-trivial irreducible representation space V λ of U(d) does not contain trivial irreps of K, W λ = 0. In contrast, if a non-trivial irrep λ of U(d) contains trivial irreps of K, then the matrix elements of W λ corresponding to the trivial irreps of K are one, and the others are zero.
Trivial irreps of K in a non-trivial irrep λ of U(d) were studied in a great detail since (K, U(d)) is an example of a Gelfand pair [97,98]. It is known that the irreps of U(d) indexed by λ ∈ Λ sph (d 1 , d, t) contains only one trivial irreps of K, and that other irreps of U(d) contain no trivial irrep of K [80]. Since trivial irreps are onedimensional, we denote by |w λ ∈ V λ a unit vector that spans the trivial irrep of K in the spherical representation λ ∈ Λ sph (d 1 , d, t) of U(d). Then, we have Importantly, for any λ ∈ Λ sph (d 1 , d, t), there exists at least one R λ ∈ U(d) such that w λ |M λ (R λ )|w λ = 0. This follows from the fact that where we have used the unitary invariance of H(d) and that the irrep ρ λ is non-trivial, so that E U ∼H(d) [ρ λ (U )] = 0. Due to the intermediate value theorem, there always exists at least one unitary R λ ∈ U(d) such that w λ |M λ (R λ )|w λ = 0. Using such R λ ∈ U(d) and Eq. (66), it is straightforward to observe that W λ M λ (R λ )W λ = 0. Furthermore, it follows that We, hence, obtain Thus, the finite set of unitaries W d1⊕d−d1 R λ W d1⊕d−d1 satisfies the condition for the design, i.e., Eq. (63) for any non-trivial irrep λ ∈ Λ ≤ (d, t), which leads to the statement that the set of unitaries defined by is a strong t-design on U(d).
Finally, let us clarify the relation between R λ and the zero of the zonal spherical function. To this end, we first observe that the matrix element w λ |M λ (U )|w λ of M λ is the zonal spherical function Z λ (U ). This can be checked by a simple calculation: for any W 1 , W 2 ∈ K and U ∈ U(d), we have where the averages are all taken over . Thus, Z λ (U ) is bi-K-invariant, and so, is the zonal spherical function. This implies that R λ is indeed a zero of the zonal spherical function. Based on this fact, we can provide a matrix form of R λ in the fixed basis in which a unitary in W d1⊕d−d1 is represented as U ⊕ V . To this end, it is important to notice that the bi-K-invariance of the zonal spherical function implies that it is characterized by the cosets of The cosets can be further identified with d 1 -dimensional subspaces corresponding to the support on which U(d 1 ) acts. For instance, the identity element in the coset of K corresponds to the subspace V 0 spanned by the first d 1 vectors of the fixed basis. The matrix form of R λ is obtained by specifying the relation between V 0 and the subspace corresponding to another representative of the coset.
To characterize the relation between two subspaces, we use the principal angles. For two subspaces X and Y , let us refer to θ = min argcos| x|y |, where the minimum is taken over all unit vectors |x ∈ X, |y ∈ Y , as the minimum angle between X and Y . The principal angles (θ 0 , . . . , θ m−1 ) between two m-dimensional subspaces X and Y are then defined as follows: θ 0 is the minimum angle between X and Y , and θ i+1 is the minimum angle between X ∩ span{|x 0 , . . . , |x i } and Y ∩ span{|y 0 , . . . , |y i }, where (|x j , |y j ) is a pair of the unit vectors that leads to θ j .

VIII. HIGHER-ORDER RB
In this section, we investigate the higher-order RB in detail. We begin with a preliminary in Subsec. VIII A and explain several basic properties of the self-adjointness in Subsec. VIII B. We consider the t-RB for general t and the 2-RB in Subsecs. VIII C and VIII E, respectively.

A. Liouville representation
2 be normalized Pauli operators on one qubit, where normalization is in terms of the Hilbert-Schmidt inner product. For q qubits, we introduce a vector n = (n 1 , n 2 , . . . , n q ) (n i ∈ {0, 1, 2, 3}) and use the notation that We also denote 2 q by d in this section. The Liouville representation is a matrix representation of quantum channels, also known as the Pauli transfer matrix. See, e.g., Refs. [84,89,99]. Let |· be a linear map from a set of all linear operators on a ddimensional Hilbert space to a d 2 -dimensional vector space that specifically maps σ n to a canonical orthonormal basis vector e n . Since the map is linear, we have for any linear operator A. Note that A|B = Tr[A † B]. Based on this vector representation of linear operators, a linear supermap E can be represented by a matrix. The Liouville representation of a linear supermap E is defined by which is a regular matrix of size d 2 . The matrix element in the canonical basis of {|σ n } n is given by The vector and Liouville representations satisfy the following properties: Properties of a linear supermap E can be also expressed in terms of the Liouville representation. For instance, the linear map E is TP if and only if (L E ) 0 0 = 1 and (L E ) 0 n = 0 for any n = 0. Since we are interested in the CPTP map E that represents a noise, its Liouville representation is always in the form of where 0 is a row vector of length d 2 − 1 with all elements being zero, α E is a column vector of length d 2 −1, called a non-unital part of the noise, andL E is a (d 2 −1)×(d 2 −1) matrix. The non-unital part α E of the noise is the zero vector if and only if the map E is unital, i.e., E(I) = I with I being the identity operator.
In the Liouville representation, the fidelity parameter f (E) and the unitarity u(E) of a noisy CPTP map E are given by respectively.

B. Properties of the self-adjointness
For a CPTP map E, the self-adjointness H(E) and the self-adjointness parameter h(E) are defined by where E (ρ) = E(ρ − I/d), and the last line is shown in Appendix D.
We first show the relation between H(E) and h(E), i.e., Eq. (28) in Subsec. IV A: From the definition of H(E), we have By rewriting E with E , the first term in the right-hand side is expressed in terms of the unitarity u(E), such as which leads to Similarly, we obtain from the facts that L E † = L † E and that |α E † | = 0 for any TP map E.
From the definition of E , it is straightforward to show that the self-adjointness parameter h(E) is given by Combining these altogether, we arrive at The self-adjointness parameter also satisfies the following properties. They are all shown in Appendix D. 4. the average gate fidelity F (E) is bounded from above by u(E) and h(E):

C. A general expression for the t-RB
We here show that the expectation value V (t) (m, E|O ini , O meas ) in the t-RB has a general form of whereÂ λ is a regular matrix depending on O ini and E(O meas ), andĈ λ (E) is a regular matrix depending only on E. As we will see below, λ labels the irreps of a t-copy representation of the unitary group, and the size of the matrices is equal to the multiplicity of each irrep.
The expectation value is defined by where E U i is the average over all unitary sequences In terms of the Liouville representation, we have Noticing the t-th power and the fact that each unitary is independently chosen from a unitary 2t-design U 2t , we obtain where L av is defined by The last line follows since U 2t is an exact unitary 2tdesign.
To write down L av explicitly, let us consider the tensort Liouville representation given by where GL(K) is the general linear group acting on the d q -dimensional vector space K defined by |σ ns : n s ∈ {0, 1, 2, 3} q , s ∈ [1, t]}.
(110) We denote the irreducible decomposition by where λ labels the irreps, and m λ is the multiplicity of the irrep labeled by λ.
The key observation is that which simply follows from the unitary invariance of the Haar measure. This implies that L av ∈ End U (K), where End U (K) is a set of all endomorphisms of K that commute with the tensor-t Liouville action of U(d). It is wellknown that End U (K) is isomorphic to the direct sum of matrix algebras: where M (m λ , C) is a set of all m λ × m λ matrices over C. Thus, the operator L av ∈ End U (K) can be represented by a direct sum of matrices.
To obtain the explicit form of L av , let K whereĈ λ (E) ∈ M (m λ , C). Each elementĈ λ (E) is given by where we have used the irreducibility in the fourth line. Consequently, it follows that Substituting this into Eq. (105), we obtain where the m λ × m λ matricesÂ λ are given by This completes the proof.
D. The first-order RB Let us briefly overview the 1-RB using an exact unitary 2-design, namely, the standard RB. We also explain how the result changes when the 2-design is an approximate one rather than the exact one.
In the 1-RB, the representation space is given by We need to find a irreducible decomposition of K under the action of a unitary group U(d) as U → L U . The Liouville representation L U is defined by L U |ρ = |U ρU † . Hence, K is irreducibly decomposed to where Denoting by Π 0 and Π 1 projectors onto K 0 and K 1 , respectively, we have where f (E) is the fidelity parameter. Note that U 2 is an exact unitary 2-design. We thus obtain that where When the 2-design is an approximate one U ( ) 2 , Eq. (128) holds only approximately. The degree of approximation depends on how we measure it, but we here assume that the design is -approximate when Eq. (128) holds up to -approximation. That is, we assume that where ∆ is some operator of O(1). Note that standard definitions of approximate designs require harder criteria (see, e.g., Ref. [78]). In this case, instead of Eq. (129), we have where Comparing Eqs. (129) and (132), we observe that using approximate unitary 2-designs result in more complicated form or the fitting function.

E. The second-order RB
We now focus on the 2-RB. Although the representation space in this case is where we have used the notation that σ n1⊗ n2 = σ n1 ⊗σ n2 , it is not necessary to consider the whole space because we assume that the initial operator ∆ is traceless. This, together with the fact that the noise map is tracepreserving, implies that the operator remains traceless during the whole process. We also observe that the whole process is symmetric under the exchange of the first and the second spaces, each labeled by n 1 and n 2 in Eq. (137). Hence, in the analysis of the 2-RB, the relevant space is only the traceless symmetric subspace defined by K T S := span{|σ n1⊗ n2 + σ n2⊗ n1 : n 1 , n 2 ∈ {0, 1, 2, 3} q , ( n 1 , n 2 ) = ( 0, 0)}, (138) where 0 = (0, . . . , 0). The irreducible decomposition of K T S can be obtained by an extensive use of the result in Ref. [100] (see Appendix C), based on which we explicitly compute V (2) (m, E|∆, O meas ). It turns out that the situation differs depending on whether q = 1 or q ≥ 2. We, hence, deal with the two cases separately.

2-RB in a single-qubit system
When q = 1, the irreducible decomposition of K T S is given by which is multiplicity-free. Here, K 0 and K 1 are respectively, with S n,m := (σ n ⊗ σ m + σ m ⊗ σ n )/ √ 2. It is obvious that d 0 := dim K 0 = 1 and d 1 := dim K 1 = 5. This decomposition implies that the expectation V (2) (m, E|∆, O meas ) is in the form of where both A λ and C λ are given by with Π λ being the projections onto K λ . Since the projections can be explicitly constructed from Eqs. (140) and (141), we can compute C λ (E). First, we have = u(E).
For C 1 (E), we start from the relation that where Π sym is the projection onto the symmetric subspace of K ⊗2 . The projection Π sym is also expressed by (I+F)/2. Here, I is the identity operator on K ⊗2 and F is the swap operator on K ⊗2 defined by 3 n,m=0 |σ n σ m |⊗ |σ m σ n |. Using the swap trick, we have Moreover, from the direct calculations, we obtain Tr[|S 0,n S 0,n |L ⊗2 E ] = L 00 L nn + L 0n L n0 , where L nm = σ n |L E |σ m . Further using the relations we obtain from Eq. (148) that Altogether, we obtain 2. 2-RB in a multi-qubit system For a multi-qubit system (q ≥ 2), the traceless symmetric subspace is decomposed into four irreducible subspaces: where K 0 = span{ n = 0 |σ n⊗ n } and the others are given in Appendix C. Each irrep is multiplicity-free. We denote by D λ the dimension of each subspace, which are Since the decomposition is multiplicity-free, V (2) (m, E|∆, O meas ) is a sum of four exponentially decaying functions. Furthermore, from the fact that K 0 = span{ n = 0 |σ n⊗ n }, we obtain that C 0 (E) = u(E). Hence, we have with Π λ being the projections onto the irrep K λ . It is not clear whether each C λ (E) (λ = I, II, III) has a clear physical meaning, such as C 0 (E) = u(E) being the unitarity. However, a linear combination of them does. To see this, we use the relation that where S 0 n := (σ 0⊗ n + σ n⊗ 0 )/ √ 2. From this relation, we can show, by a calculation similar to the one-qubit case, that We finally note that Tab. I in Subsec. IV E is obtained by constructing the orthonormal basis in each subspace K I , K II , and K III (see Appendix C). We also assume that the noise E is weak, so that E(O meas ) ≈ O meas . Based on this assumption, we have enabling us to compute A λ for given initial and measurement operators.

IX. EXPERIMENTAL REALIZATION OF 2-RB
Based on the former sections, we explain in detail how we have experimentally implemented the 2-RB and estimated the self-adjointness of the noise in the system. In Subsec. IX A, we provide the details of the numerical evaluation of the one-and two-qubit 2-RB discussed in Sec. V A. The details of experiment is given in Subsec. IX B.
A. Numerical analysis

Single-qubit systems
We explain the fitting procedure of the 2-RB on onequbit systems in detail. The noise we consider is given by the following CPTP map: which is characterized by three parameters p, q, and θ. We particularly choose θ as p = sin 2 θ for the fidelity parameter f (E 1 ) to be independent of the coherence parameter q. Using the Liouville representation, it is straightforward to compute the fidelity parameter, the unitarity, and the self-adjointness parameter of this noise. They are, respectively, given by As shown in Theorem 7, V (2) for one qubit is To obtain Fig. 2 in Subsec. V A, we first estimate the fidelity parameter f (E 1 ) from the 1-RB, i.e., by fitting V (1) (m, E 1 ||0 0|, |0 0|) and then obtain u(E 1 ) and h(E 1 ) from the fitting results of V (2) (m, E 1 |Z, |0 0|) and the estimated value of f . The fitting of V (2) is first performed based on Eq. (173) by regarding the coefficients A 0 , A 1 and the two exponential decaying rates as free parameters. Then, we subtract the first exponential curve of Eq. (173) from the data and carry out the fitting of the second exponential curve again where we consider A 1 and the base of the second exponential curve as free parameters. This procedure is redundant, but turns out to improve the accuracy of the fitting because, in most cases, the first exponential decaying rate is larger than the second one, and thus the second exponential curve is clearly visible in the region of long sequence length.
We also estimate how many measurements and unitary sequences suffice to obtain a good estimate of the noise parameters. In Fig. 6, we provide V (2) for various numbers of measurement and samplings of unitary sequences. We observe that 1000 for both suffice to obtain a good estimate when p = q = 0.02, which corresponds to the gate fidelity 97.3%.
When the number of measurements is small, experimental values are positively biased from theoretical val-ues with the infinite number of samplings (see the figures in the top line of Fig. 6). This difference can be understood as follows. Let O meas Oini,i,n be the expectation value for a random sequence described by i with the finite number n of measurements. We describe the expectation value and variance of this random variable averaged over all choices of the unitary sequence as µ and σ 2 /n, respectively. Note that this mean value is independent of n, and this variance is inverse proportional to n since O meas Oini,i,n is a linear combination of binomial distribution. V (2) is the expectation value of squared random variable O meas 2 Oini,i , and its expectation value over random sequences is derived as Therefore, an experimentally obtained value with a finite number of measurements is positively biased by σ 2 /n. In practice, we can remove the effect of this bias by increasing the number of sampling n and using the region satisfying µ 2 σ 2 /n for fitting. We finally check the robustness of the 2-RB on single- qubit systems against SPAM errors. Although the 2-RB is ideally SPAM-error free, the fitting may become harder with the existence of SPAM errors. Our analysis, however, reveals that this is unlikely the case. Here, we model the state-preparation error η prep as ρ = η prep |0 0| + (1 − η prep )|1 1| and ρ = (1 − η prep )|0 0| + η prep |1 1|, and measurement error η meas as readout bitflip error, i.e., the POVM is {Π x |x ∈ {0, 1}} where Π x = x∈{0,1} η meas (x |x)Π x and η meas (x |x) is conditional probability. In this numerical experiments, we assume that η = η prep = η meas (0|1) = η meas (1|0) and get 1000 samples and 1000 random sequences. We set the parameters p = q = 0.02. The results are provided in Fig. 7 and Tab. IV. The relative errors of estimates for F, u, and H are within 5% except for the estimate for u when η = 0.3, thus the 2-RB is likely to work well even in realistic situations with SPAM errors as expected from the analytical studies.

Two-qubit systems
For two-qubit systems, we investigate the noise given by where we set θ to p = sin 2 θ. From Theorem 7, V (2) in this case shall be in the form of and C λ (E 2 ) satisfy With this setting, theoretical values are derived as follows.
from which the theoretical value of the self-adjointness parameter h(E) can be computed from Eq. (178). A step-by-step fitting process of two-qubit systems is shown in Fig. 8 for p = 0.01. The numbers of measurement and samplings of unitary sequences are both set to 10 4 . The fidelity parameters f (E 2 ) can be obtained from V (1) (m, E 1 ||00 00|, |00 00|), which are plotted in the top line of the figure. In the figure, dashed lines are theoretical values, and the shaded area represents a standard deviation of each data. When q is near to unity, standard deviations are large even when the sequence length increased. This is because the final quantum state is nearly pure state when q ∼ 1, and probability distributions fluctuate randomly according to chosen random sequences. On the other hand, when q ∼ 0, the final quantum state quickly converges to the maximally mixed state, and thus probability distribution becomes independent of the chosen random sequences. We then fit four values u, C I , C II , C III step by step where initial and measurement operators ∆ and O meas are chosen according to Tab. I. The fitting results for several q are shown in Tab. V.
The obtained results are shown in the second line of Fig. 8. The sampled data points and fitting results are shown as blue points and dashed lines, respectively. These lines are linear combinations of two exponentially decaying function. Exponential decays with coefficient u and C I are shown as orange and green lines, respectively. While we can clearly see two exponential decays for coherent noise, i.e., in the case of q ∼ 1, an exponential decay of C I part where these initial and measurement operators are chosen from the third and fourth columns of Tab. I. The obtained results are shown in the third and fourth lines of Fig. 8. In each figure, numerical data is plotted as blue circles and fitting results are shown as dashed lines. In each fitting process, only a single exponentially decaying term is unknown in advance. We showed unfitted exponential decay as orange circles and fitting results as dashed lines. Although accuracy of orange data becomes not reliable when its value becomes much smaller than the others, we can fit C II and C III reliably.
We calculate the averaged fidelity F (E), unitarity u(E), and self-adjointness H(E) from the fitting results. The processed values are shown in Tab. VI. We evaluate relative errors for all the plots with the same method as the single-qubit 2-RB, and confirm that the relative errors are below 4% for all the cases of p = 0.01, 0.02, 0.1, 0.2 except the case when the theoretical value is exactly zero. While the relative errors of the self-adjointness become a few tens of percent in the case of p = 0.4, we can say that reliable values can be obtained when the fidelities of oper- ations are sufficiently high. Note that, although standard deviations of the fitting results become large when u is almost equal to C I , we confirmed that the fitted results are close to theoretical values even in such cases. Thus, we conclude that 2-RB works also for two-qubit systems.

B. Details of the experiments
In this section, we provide the details of the experiments in Sec. V B. A superconducting qubit can be regarded as a sort of the LC resonant circuit, where a Josephson junction is an effective inductance, and has the Hamiltonian equivalent to that of the one-dimensional free particle trapped in anharmonic potential. The parameter fields of the qubits are summarized in Tab. VII.
All unitary gates required in the 2-RB for a singlequbit case were implemented by two R X (π/2) gates and three R Z gates with an arbitrary rotation angle, which are implemented by the shaped microwave pulse (Half-DRAG) with the length 11.70 ns [101] and the Virtual-Z gates [102], respectively. The pre-measured averaged gate fidelity of the single-qubit Clifford gate is 0.991.
The single-shot qubit readout was performed via the impedance-matched Josephson parametric amplifier [103], and the assignment fidelity of the readout was 0.943.
As supplemental experiments, T1 decay and Ramsey oscillation were observed to clarify the background noise Ramsey oscillation experiment. The horizontal axis represents the detuning from the qubit eigenfrequency, and the vertical axis represents the power spectrum of the Ramsey oscillation. The points connected by the blue line represent the experimental data, and the orange line represents the fitting curve.
source of Q1. The experimental results of the T1 decay is shown in Fig. 9, where the horizontal and the vertical axes represent the delay time and the projected value of the IQ readout signal, respectively. The blue and orange dots represent the experimental results with and without a flip operation just before measurement, respectively. The lines provide the fitting curves. As seen from the result, the T1 decay of Q1 follows exponential behavior, which is consistent with the expected behavior in isolated qubits.
The experimental results of the Ramsey oscillation is also given in Fig. 10, where the horizontal and vertical axes represent the detuning from the qubit eigenfrequency and the power spectrum of the Ramsey oscillation, respectively. The points connected by a blue line represent the experimental data, and the orange line provides a fitting curve. In the fitting, we did not take the data in the small power spectrum region (< 1 V 2 ) into account. This is because the noise floor derived from the white noise is dominant there. As seen from the result, the power spectrum of the Ramsey oscillation has no peaks other than the qubit eigenfrequency. From the fitting curve, it was found that the detuning from the qubit eigenfrequency ∆ (MHz) and the power spectrum of the Ramsey oscillation P S(∆) (V 2 ) are related as P S(∆) ∝ ∆ −2.004 .
This is consistent with the expected behavior when the transmon qubits are isolated well.
From the results of these supplemental experiments, we conclude that Q1 is not in the strong coupling regime with any noise source, which implies that the background noise of Q1 is time-independent. Thus, the requirements for 2-RB are met in our experimental system.

X. SUMMARY AND DISCUSSIONS
In this paper, we have provided an explicit constructions of exact unitary t-designs for any t. In particular, quantum circuits for exact unitary t-designs on N qubits have been provided for the first time. Our construction is inductive with respect to the number of qubits. Hence, all constructions obtained in this paper are inefficient when the number of qubits is large, implying that it is of practical use only when the size of the system is small.
As an application of exact unitary 2t-designs on a small system, we have proposed the t-RB, which enables us to experimentally estimate higher-order properties of the noise on a quantum system. Since the unitary designs are used in multiple times in a single run of the protocol, it is important for the design to be exact. After providing a general scheme of the t-RB, we have studied the 2-RB in detail. It was shown that the 2-RB reveals the self-adjointness of the noise, a new characterization of the noise that we argue to play an important role in QEC especially when decoders are based on applications of Pauli operators. Our results have been demonstrated numerically, which shows that the 2-RB is experimentally tractable. We have then experimentally implemented the single-qubit 2-RB on the superconducting qubit system. From the experimental results, we found that the characteristics of the background noise of a qubit changes depending on the presence of the interaction with the adjacent qubits.
Our results open a number of future problems. Regarding the implementation of t-designs, it is important to improve the efficiency. Despite that the inefficiency in our construction is likely to be intrinsic due to an inductive nature of the construction, the representationtheoretic method provides a way to searching more efficient ones. More specifically, the key in the construction is the relation between the representation of the whole unitary group and that of a certain subgroup of the unitary group. This indicates that finding the construction of exact unitary designs may be reduced to the problem of searching for a subgroup whose representation has a good relation to that of the whole unitary group.
It is also important to further develop the theory of the t-RB protocol. In this paper, we have analyzed only the 2-RB in detail. It then turns out that self-adjointness of the noise can be revealed. It is of great interest to concretely investigate what characterization of the noise can be generally obtained from the t-RB. In the context of QEC, it is also important to comprehensively analyze quantitative relations between the self-adjointness and the feasibility of QEC. Another promising future problem is to use exact higher-designs in the other RB-type protocols. The t-RB is a straightforward generalization of the standard RB. However, there are numerous variant protocols [69], most of which, if not all, are based on the Clifford group that is an exact unitary 2-design. By extending such protocols to those with higher-designs, the noise on the system can be characterized in more detail.

h(E) = u(E) if and only ifL
where F := n σ ⊗2 n is the swap operator on K ⊗2 . Note that, in the last line, we have used the fact that E is TP, also implying that E † is unital. We now use another expression of the swap operator, which is