Quantum circuit simulation of superchannels

Quantum simulation is one of the central discipline to demonstrate the power of quantum computing. In recent years, the theoretical framework of quantum superchannels has been developed and applied widely as the extension of quantum channels. In this work, we study the quantum circuit simulation task of superchannels. We develop a quantum superchannel simulation algorithm based on the convex decomposition into sum of extreme superchannels, which can reduce the circuit cost. We demonstrate the algorithm by numerical simulation of qubit superchannels with high accuracy, making it applicable to current experimental platforms.

Quantum simulation is one of the central discipline to demonstrate the power of quantum computing. In recent years, the theoretical framework of quantum superchannels has been developed and applied widely as the extension of quantum channels. In this work, we study the quantum circuit simulation task of superchannels. We develop a quantum superchannel simulation algorithm based on the convex decomposition into sum of extreme superchannels. We demonstrate the algorithm by numerical simulation of qubit superchannels with high accuracy, making it applicable to current experimental platforms. Our study stands as an expansion of the superchannel theory to the field of quantum simulation and algorithm, as well as an extension of quantum simulation from channels and open-system dynamics to superchannels and processes with manifest quantum memory effects.

I. INTRODUCTION
In modern quantum physics, the control and engineering of complex quantum dynamics is essential. In recent years, quantum computing has become a powerful paradigm to achieve this. With the standard quantum circuit model, the field of digital quantum simulation is expected to be powerful to realize and study quantum dynamics.
A quantum circuit contains a sequence of gates realized at discrete time slots. The original and still an important task is quantum Hamiltonian simulation, with each gate in a circuit usually described by a local Hamiltonian term [1]. At the same time, quantum circuit simulation relies on quantum gates directly and no explicit Hamiltonian. It is more flexible and of digital nature, and applies to broader settings such as quantum optical systems.
There are scenarios beyond the common setting of channels. A mathematical milestone is the development of quantum superchannel (and comb) theory [22][23][24] based on channel-state duality [3,25]. It has been used * wds@itp.ac.cn in quantum metrology [26], computing schemes beyond circuit model [27], and recently used in resource theory of channels [28], non-Markovian dynamics with quantum memory effects [29], and quantum algorithms [30]. Therefore, it is worthy to explore more applications of superchannels in quantum computing.
The quantum channel simulation approach developed in [7][8][9] is based on the convex-decomposition of channels into sum of extreme ones. The convexdecomposition of channels relies on the convex geometry of channels acting on a quantum system. The nontrivial feature is that extreme channels can have ranks higher than one [3,[31][32][33][34][35][36], making the convex-decomposition difficult. Actually, there is a notable open problem called Ruskai's conjecture [33] which concerns the upper bound for the number of (generalized) extreme channels needed. This has motivated mathematical investigation on the convex geometry [34,35]. The algorithmic approach developed in [7][8][9] has numerically verified Ruskai's conjecture for low-dimensional cases (up to four). Also the circuit cost is smaller than a direct Stinespring's dilation method. If the ancilla is recyclable and refreshable, the ancilla can be as small as a qubit, which adaptively controls a sequence of operations to simulate a channel [37][38][39][40].
In this work, we extend quantum channel simulation to quantum circuit simulation of superchannels, or quantum superchannel simulation. It is the expansion of superchannel theory to quantum simulation and algorithm on the one hand, and on the other hand, it is also the extension of quantum channel simulation to superchannels. In particular, we develop the simulation scheme based on convex-decomposition of extreme superchannels. This requires the description of extreme superchannels, which has not been systematically studied before. We make a parameter counting for superchannels and combs, and present an analog conjecture with Ruskai's for superchannels. We develop numerical optimization simulation algorithm for qubit superchannels, and present simulation results for a few cases. Our results can be used to study more general superchannels and employed for various experimental implementation.
This work contains the following parts. In section II, we review the representations of quantum channels, and the convex-decomposition algorithm of channels. We revisit the exact decomposition of channels [33,41], and argue that the decomposition is non-unique and most likely can only be done numerically. In section III, we develop the representations of superchannels and extreme ones, and also the convex-decomposition algorithm, with numerical simulations for the case of qubit superchannels. We also discuss the simulation of unital channels and superchannels, and the physics of extremality (in the appendix). We conclude in section IV.

II. QUANTUM CIRCUIT SIMULATION OF CHANNELS
A. Channels and extreme channels In this section, we describe the theoretical framework for quantum circuit simulation of superchannels by reviewing the case of channels which employs a convexdecomposition optimization algorithm.
We focus on finite-dimensional Hilbert spaces. We first recall representations of quantum channels briefly. For a Hilbert space H , quantum evolution is in general described by CPTP maps [2], or known as quantum channels for K i as Kraus operators [2] with i K † i K i = 1, and r as the rank of E, ρ ∈ D(H ). Recall that the Kraus representation above is not unique, and the rank is the minimal one among all equivalent representations. As an example, unitary evolution is rank 1 with U † U = U U † = 1 for U ∈ SU (d), with d = dim(H). The notable channel-state duality [3,25] maps a channel E ∈ L (D) into a quantum state i=0 |i, i known as a Bell state or ebit. The rank r is concisely characterized as the rank of the Choi state ω E . The Stinespring's dilation theorem [42] guarantees that a channel can be described as an isometry V with V = i |i K i , and further dilated to a unitary U with V = U |0 as the first block column of U , and |0 as the initial ancillary state. We will refer to this as the quantum circuit representation of channels.
An important feature of channels is the convexity of their set. A convex set has extreme points, which cannot be written as any convex combination of others. It is clear to see unitary operators are extreme, but there are more extreme channels with higher ranks. For the set of qudit channels S d , Choi [3] proved that a channel E is extreme in S d iff {K † i K j } is linearly independent for K i as its Kraus operators. A consequence is that the rank of a channel necessarily satisfies r ≤ d in order for being extreme, while for general channels the rank is upper bounded by d 2 . Define the set of rank up to d channels as S g d , and such channels are called generalized extreme channels [33], or "gen-extreme" channels following Ref. [43]. A gen-extreme channel that is not extreme is called a quasi-extreme channel, which clearly can be written as a convex sum of at least two extreme channels.
The mathematical conditions above for being extreme or gen-extreme are not enough for applications. We still need to find formulas of them. Here we review the representations of extreme channels, including the Choi state form [8] and quantum circuit form [8,9,40,43,44].
Choi state is a bipartite matrix, which allows a natural block-matrix form representation. There is a concise block-matrix form of gen-extreme channels [8] as for C i ≡ C ii ≥ 0, and U i , U j ∈ U(d). For general channels, the unitary operators above would be replaced by contractions [8]. Also a peculiar feature of gen-extreme channels is that the form above can be used directly to find the quantum circuit for the channel, without solving the Kraus operators [43]. This is generically not possible for general channels without knowing the Kraus operators. Namely, define an isometry V : for ρ = ij ρ ij |i j|. Now the channel is realized by a quantum circuit W with for W = swap · U * , and the swap gate, which is applied after U * , is on the system and ancilla which are of the same dimension. The unitary U is defined such that U |0 = V . The Kraus operators are derived as K i = i|W |0 .
Knowing the quantum circuit W is not enough yet since we still need to decompose it into product of elementary gates. In quantum computing, we often use the universal gate set {cnot, U(2)}, for the controlled-not gate cnot and the group U(2) on qubit. Matrix decomposition techniques have been employed to design quantum circuits for channels, including the cosine-sine decomposition (CSD) and QR decomposition [8,9,40,44]. A qubit gate can be further decomposed as product of Hadamard gate and T gate for fault-tolerance [4], but in this work we do not require this. To quantify the complexity or cost for executing a circuit, we refer to the circuit cost as the size of required ancilla, and the total number of cnot and qubit gates in a circuit, which scales linearly with the number of free parameters in the unitary operator for the circuit [44]. Also for the setting of our algorithm, we mostly use the number of free parameters to characterize the complexity of a circuit or circuit family.

B. Convex decomposition of channels
Given the convex body of channels, we expect a channel can be decomposed as a sum of extreme channels, just as a state decomposed as a mixture of pure states. However, convex decomposition of channels is difficult, and it is still an open problem in the field of quantum information since 2007 [33]. Namely, Ruskai and coauthors conjectured that any channel E ∈ S d can be decomposed as The uniform probability in the sum (7) can be relaxed to general mixing with probability p i ∈ [0, 1] and d i=1 p i = 1. The case for d = 2 has been proven [41], also revised in the next section, the method of which can be extended to the case for dimension-altering channels from a qudit to a qubit [33], a useful fact that will be used below for our simulation of qubit superchannels. This conjecture is numerically supported for dimensions up to four [8,9]. The classical optimization algorithm aims at for θ as the angle parameters in ω E , which is a sum of gen-extreme channels, and ω E as a random input Choi state for a channel E. The objective function D t (ω E , ω E ) is the trace distance between states, upper bounded by 1. Practical simulation finds that p i can be chosen to be uniform, which is suitable for experimental implementation [16]. It is unlikely that there exists an exact decomposition for general channels. The Carathéodory theorem on convex sets [45] provides the upper bound d 4 − d 2 as the number of extreme channels, which is too loose in this case. As the rank of extreme channels is not fixed, it is hard to characterize them rank by rank. Ruskai's conjecture is also meaningful from a heuristic parameter-counting. A general qudit channel contains parameters, hence it is reasonable to require a mixing of d gen-extreme channels. Compared with eigenvalue decomposition of matrices, the convex decomposition of channels can be viewed as a sort of generalized eigenvalue decomposition. A gen-extreme channel can be mapped to an isometry V as described above, which can be further treated as a pure state. The difficulty of this decomposition shall be comparable with the eigenvalue decomposition of Choi state directly to find the set of Kraus operators. The merit for quantum simulation is to reduce the circuit size, from O(d 4 ) gates to O(d 3 ) gates and from three qudits to two qudits, which would benefit experimental implementation in the near future.

C. Example: qubit channel decomposition
Two decades ago, Ruskai, Szarek, and Werner [41] first proved that a qubit channel can be decomposed into a sum of two gen-extreme qubit channels. We will refer to it as the RSW method. Much later, the method based on cosine-sine decomposition (CSD) is developed [9] for an optimization algorithm. Within this framework, we emphasize that the convex decomposition is not unique and generically can only be done numerically, i.e., cannot be done analytically for arbitrary channels.
The RSW method mainly contains two ingredients: an exact decomposition based on the Choi state and a canonical form of extreme qubit channels based on the affine representation of channels. Up to a pre and a post unitary gate, it shows that a gen-extreme qubit channel is specified by two Kraus operators with α, β as rotation angles. Based on the block-matrix form (see Sec. II A), any qubit channel E can be written as an average of two gen-extreme qubit channels since a two-by-two contraction matrix A can be decomposed by singular-value decomposition as an average of two unitary matrices The two gen-extreme qubit channels in the RSW decomposition have the same diagonal blocks C 0 and C 1 in their Choi states. From parameter-counting, a qubit channel needs up to 12 parameters, but only 8 for genextreme ones. From the Choi state, there are 4 parameters for C 0 and 8 parameters for a contraction A, but 4 parameters for a unitary U ∈ U(2). This shows how the parameters split under the decomposition. In addition, this method extends to channels that map qudit to qubit systems [33] due to the two-by-two block-matrix form of Choi states, but cannot be extended further to qudit channels [8].
With the map between Choi state and circuit for genextreme channels defined in Sec. II A, we find a more explicit form Despite this, the forms between F (9) and K (12) do not match. If one wants to use the F form, it still needs to find the pre U pr and post U po unitary rotations on the qubit, and also possibly a basis transformation U = [u ij ] on the ancilla. Namely, K i = j u ij U pr F j U po , which apparently cannot be analytically solved in general. However, due to the diversity of channel representations and the rich structure of the set of channels, it is still left open whether there could be a better analytical method for the convexdecomposition of generic or special types of channels. The theoretical analysis above is beneficial to understand the decomposition problem, but for the simulation algorithm we need to follow the algorithmic approach that we have developed [7][8][9]. In order to motivate the circuit simulation of superchannels, now we recall the CSD method which directly represents gen-extreme channels by quantum circuits [9,40,44]. For the qubit case, we need the unitary operators U ∈ U(4) which can be represented as . Then the Kraus operators are The above form (14) was used as an ansatz in an optimization algorithm for the convex-decomposition of qubit channels, which takes a desirable qubit channel as input and yields the parameters in the circuits directly [7][8][9]. The CSD method extends to qudit channels and we will use this for superchannels, as well.

III. QUANTUM CIRCUIT SIMULATION OF SUPERCHANNELS
A. Superchannels and extreme superchannels In this section, we review the theoretical framework of quantum n-combs [22][23][24], which include superchannels for n = 2. We will elucidate the circuit representation and present an analog of Ruskai's conjecture for superchannels. As channels can be represented as Choi states, the mappings on Choi states are further defined as superchannels. Similar frameworks have also been introduced [46][47][48][49], while here we employ the most well developed framework of quantum combs. By definition, channels are 1-combs, and states are 0-combs. However, in this work we will change the terminology slightly. We will define n-superchannels as (n + 1)-combs to simplify the notions. So channels are 0-superchannels, the usual superchannels are 1-superchannels. From Fig. 1, n-superchannels will have n circles for the input, while n + 1 unitary 'boxes' as the 'teeth' in the comb.
We focus on Choi state and quantum circuit representations, while other forms such as the affine representation can also be developed. We consider the qudit case. Label the input space for the channels as H 0 , H 2 , . . . , H 2n , and the output space for the channels as H 1 , H 3 , . . . , H 2n+1 . Let an orthonormal operator basis for each space L (H ) k be {E i,[k] } i∈ [1,d 2 ] with E 1,[k] = 1, and all others traceless. As we use ω E as Choi state for a channel E, we use ω (n) as the Choi state of an n-superchannel, which is shown [50] to be the form with arbitrary hermitian operators B i,[k] acting on ⊗ k j=0 H j . The form is recursive with (16) Before we present the quantum circuit form, we make a parameter counting for ω (n) , which has not been presented explicitly before. The free parameters are contained in the B matrices above. Let y n be the number of parameters for ω (n) , then For channels, y 0 = d 2 (d 2 − 1), so we find For superchannels y 1 = (d 4 + 1)d 2 (d 2 − 1). We find from this the dilated quantum circuit is to add one additional ancilla of dimension d 2 , see Fig. 1. Given a Kraus form of a channel E = {K i }, it is changed by a superchannelŜ aŝ and a S † a S a = 1. The new channel is represented by Kraus operators with ia F a † i F a i = 1. Note we put a hat on the symbol for superchannels, and the superscript T is transpose. The operators K m v and K ma w are Kraus operators, and their meaning can be seen clearly from the quantum circuit form of a superchannel, see Fig. 1. We have K m v = m|V |0 and K ma w = m|W |a for unitary operators V and W .
The circuit form in Fig. 1 takes the channels as input directly, plugged in the circuit 'board' [23] at suitable locations instead of starting from the initial time. This can also be equivalently realized by a circuit acting on their Choi states as input at the initial time [51], based on the channel-state duality, in particular, the map between the identity operator and the ebit state. It is apparent that superchannels are special channels, so they can be simulated by circuits for channels. However, the simulation cost might be larger. For a rank r superchannel, the size of pre V and post W unitary is of size at most rd. If we use a channel to simulate its action of Choi state, then it needs a unitary dilation U of size rd 2 . The particular form of Kraus operators S a (20) needs to be taken into account to obtain an optimal circuit. For the goal of quantum simulation in this work, it is suitable to employ the circuit form in Fig. 1.
There is no direct relation between U and the pair of V and W . However, if V is of special type, we can obtain U from the pair. From CSD, V can be written as a block form, which is a product of multiplexers [9]. Recall that a multiplexer is a type of uniformly-controlled gates, which is block-diagonal in a certain basis. Here V is bipartite acting on the system and an ancilla. A multiplexer either has the system or the ancilla as the control. So, if V only contains one multiplexer, then it can be shuffled to delete the ebits. The U is of the form U = (1 ⊗ W )(1 ⊗ V ), with V obtained from V , and the two 1s each acts on one part of the Choi state.
Furthermore, the adaptive-circuit approach developed in [37][38][39][40] reduces the ancillary cost to a single qubit for any channels, and no qubit ancilla needed for POVM, yet this approach does not apply to superchannels. The reason is that the pre and post unitary operators cannot be reduced to isometries due to the memory wire in between, which transfers quantum information from the pre to the post unitary operators. This highlights the peculiar effect of quantum memory realizing superchannels which proves to be useful in various tasks [26][27][28][29][30].
The extremality condition of n-superchannel has been proven by D'Ariano and coauthors [50], which is an algebriac condition that is more complicated than the case for channels. While we find the necessary condition is that the set {K † i K j } is linearly independent for K i as the Kraus operators for an n-superchannel. This further reflects on the ancillary dimensions in the quantum circuit form. For clarity, let us denote the ancilla dimensions by a vector, then we find the ancilla dimension for extreme n-superchannel is alternating between d 2 and d, namely, If n is even, we need to modify the final ancilla dimen-sion. For example, if n = 2, the ancilla dimension is (d 2 , d 2 ). From the parameter counting, we find a fullrank gen-extreme n-superchannel needs parameters in the order O(d 3(n+1) ). As a result, we conjecture that which extends Ruskai's conjecture for qudit channels in Eq. (7). The above conjecture is one of the central result in this work. Below we numerically support this conjecture for qubit superchannels with n = 1 and d = 2. Also it is worthy to mention that the RSW method does not directly extend to the case of qubit superchannels.

B. Convex decomposition of superchannels
We first propose a few quantum circuit forms for genextreme superchannels, and then apply them to the qubit case. In the above we know that qudit superchannels contain O(d 8 ) parameters, while gen-extreme ones contain O(d 6 ) parameters. This means it needs an ancilla of dimension d 2 for gen-extreme superchannels, as shown in Fig. 2. The post unitary W contains O(d 6 ) parameters, while the pre unitary V contains O(d 4 ) parameters. This is named as type-I circuit. We also introduce type-II and type-III circuits as shown in the figure as subclasses of type-I. The type-II circuit contains O(d 5 ) parameters, and type-III circuit contains O(d 3 ) parameters. Beyond decomposition task, these circuits can be used to generate different types of gen-extreme channels. The type-I is general, while type-II only contains a memory ancilla of dimension d, and type-III does not have a memory wire in between, namely, the two unitary V and W realize two independent channels. As such, they can be further reduced to gen-extreme channels. In terms of Choi states, it is a product state for type-III, while entangled for the other two. The blockmatrix form of Choi states is not apparent for superchannels compared with channels, so we will not present them here.
For qubit superchannels, the convex-decomposition needs four gen-extreme ones of type-I. We use CSD to represent the pre and post unitary operators V ∈ U(8) and W ∈ U(8). The quantum circuit can be easily designed using iterative CSD as in [8,9,40]. A circuit form can also be found in [52]. If we count the number of cnot gates, they scale with the number of parameters, and this is much larger than the number for genextreme qubit channels, which is only one. Therefore, we see there is a big difference between channels and superchannels.
The type-II circuits have an interesting connection with channels. For a general qubit superchannel, from the viewpoint of W , it can be viewed as the dilation of a channel from three-qubit to a qubit. From Ruskai's proof [33], it can be decomposed into a sum of two genextreme channels from three-qubit to a qubit, which only needs a qubit ancilla. Namely, the ancillary dimension vector is (4,2). Denote the set of these rank-8 channels as S 8 , and they are not gen-extreme qubit superchannels. This means they can be further decomposed into type-I or type-II circuits, which have ancillary dimension vector (4, 0) and (2, 2), respectively. However, we expect no exact decomposition of circuits of the form in S 8 . Therefore, in our numerical simulation algorithm we take an anonymous approach that employs the most general type-I circuits to represent genextreme qubit superchannels, and take arbitrary qubit superchannels as input.

C. Unital channels and superchannels
Besides the generic case, it is also important to see if there could be simpler scheme for special types of input superchannels. For the important class of unital ones, we find, different from channels, it is mostly not the case for superchannels. We first recall the case of unital channels. A channel is unital iff it preserves the identity operator 1. Unital channels form a convex set, denoted as S u d , and different from extreme channels, a product of unital channels is still unital. The extreme points for unital channels have also been characterized [31], which requires that {K † i K j K j K † i } is linearly independent. However, this condition leads to an upper bound of rank for gen-extreme unital channels even larger than the usual case. No circuit form for extreme unital channels has been known. Rather than using extreme unital channels, it turns out there is another scheme that is more suitable for simulation. It is proved that [53] a channel E is unital iff it is an affine combination of unitary channels with i a i = 1, and a i ∈ R. This is a quantum extension of the Birkhoff's theorem, which states that doubly stochastic matrices are in the convex hull of permutations. However, extreme unital channels can have rank higher than 1. Mixing-unitary (MU) channels only forms a subset of unital channels. Namely, we have For qubit channels, S mu 2 = S u 2 . For qutrit, there is a simple example of unital but not MU channel, whose Kraus operators are from the generators of SO(3) in the spin-1 representation [31]. Therefore, all unital channels can be simulated by sampling unitary operators, with possible negative probability. There is no need of a quantum ancilla, which means all initial quantum information in the initial state ρ is preserved. This implies that unital channel preserves or does not contain quantum information. So it is reasonable that unital channels preserve 1, i.e., the junk state.
For a non-unital channel E, it changes the input 1, namely, it can create quantum information E(1). It requires an ancilla, whose final state is It is a Gram matrix of the set of vectors {K i ⊗ 1|ω }.
The Gram matrix is of full rank if {K i } is linearly independent. Both E(1) and σ can be used to quantify the non-unital feature of a channel. It is interesting to observe that, physically, being unital is quite distinct from being extreme. The formula (26) above is better understood from its complementary channel, discussed in the Appendix A. Now we study unital superchannels. It turns out the structure is more complicated, and different generalizations of unital channels can be defined [28]. They can be completely specified by conditions on the Choi states of superchannels which take the form [3] ⊗ B i, [2] .
For simplicity, we use R as the symbol for a superchannel and R [i] to denote its marginal on the space i. A set of identity-preserving (IP) superchannels is defined with which is denoted as condition (a). A set of doublystochastic superchannels is defined with (a) and also which is condition (b) and requires V to be the dilation of a unital channel. A set of unital-channel preserving superchannels is defined with (b) and also  which is denoted as condition (c). A set of MU superchannels, which is the smallest amongst here, is defined when V and W are unitaries with a common source of random bits. Different from channels, there is a memory wire between V and W to realize superchannels, leading to the various cases above. For unital and even extreme unital superchannels, the memory wire cannot be replaced by classical bits in general. The condition for being extreme for each type of the unital superchannels presented above is unknown. As for the case of unital channels, we expect such conditions would not help directly for their circuit representations. However, we can a priori define a set of affine-mixing-unitary (AMU) superchannels as an extension of MU superchannels, with V and W replaced by affine-mixing of unitaries, possibly with different sources of mixing. The AMU set apparently belongs to the IP set above, therefore does not capture all extreme unital superchannels. At present, it is not clear if the quantum Birkhoff's theorem [53] can be extended to unital superchannels. For the purpose of quantum simulation, the wise strategy is to take the input superchannel anonymously, no matter it is unital or not, and then use optimization to find the convex sum of gen-extreme superchannels.

D. Convex decomposition algorithm of qubit superchannels
In this section, we present the numerical simulation results for the convex-decomposition of qubit superchannels. The simulation is summarized in Table I. Note the circuit cost is implicit based on the CSD formula and the number of parameters. For clarity, denote the set of qubit superchannels, gen-extreme, and rank-8 ones as S , S g , and S 8 , respectively, and elements of them asŜ,Ŝ g , andŜ 8 , respectively. A random input is generated by a pair of Haar-random unitary operators for a superchannel. For instance, a randomŜ is generated by Haar-random unitary V ∈ U(8) and W ∈ U(32). Then the Kraus operators {S a } for the superchannel is obtained from Eq. (20) with K m v = m|V |0 and K ma w = m|W |a for m ∈ [1,4], a ∈ [1,16]. The input Choi state ω (1) , which is a 16×16 matrix, is then obtained from {S a }. The optimization parameters are those rotations angles in the circuits for gen-extreme or rank-8 superchannels, expressed using the CSD formula. The optimization uses an optimization algorithm from Matlab, which is to perform a few local optimization first, then compare and find the best local one, and then generate a few new local samples, and iterate such a process [7][8][9].
We have simulated four decomposition tasks in our simulation as shown in the table. We pick uniform probability variables. For each task we simulated more than 20 random input instances following the standard routine in [7][8][9]. The first row is for the decomposition S = (Ŝ 8 1 +Ŝ 8 2 )/2, which can be exactly decomposed according to Ruskai's result [33]. Our simulation contains 656 parameters, coming from the CSD representation of the circuits. The simulation precision, i.e., the trace distance on Choi states, is of the order 10 −3 ∼ 10 −4 . This reveals the validity of our simulation program which is executed on a normal desktop computer. We believe that better precision can be obtained with more computational resources. The second row is for the decom-positionŜ = (Ŝ g 1 +Ŝ g 2 +Ŝ g 3 +Ŝ g 4 )/4, which uses 512 parameters and yields precision in the same order as the first row. This shows that our scheme is good enough for practical simulation of qubit superchannels, hence can be applied in experimental tasks which often involve noisy quantum gates. For instance, current noisy quantum devices [54][55][56][57] can achieve gate infidelity about of the order 10 −3 , often with single-qubit gate fidelity higher than the two-qubit one. With tens of gates in a circuit, the implementation error would dominant over the error from the numerical decomposition.
The third and forth rows in the table takes a ran-domŜ 8 as the input, and the goal is to check if it can be decomposed as a sum of twoŜ g as discussed in section III B. The simulation strongly suggests that this cannot be done numerically via our algorithm, and instead, a sum of fourŜ g is needed. Although the reverse direction, namely, a sum of twoŜ g can lead to rank-8 superchannels, this apparently only occupies a small portion of the whole setŜ 8 . Physically, this means that the last qubit ancilla for the post unitary operator in S 8 cannot be replaced by a random bit in general, and at least two random bits are needed. More generally, a low-rank (between 4 and 16) superchannel does not necessarily require fewer gen-extreme superchannels in a convex decomposition, as confirmed by the numerical simulation for both channels [8,9] and superchannels. This verifies that the rank of a channel, or completely positive map, is a 'global' feature of it, depending on its set of Kraus operators instead of individual ones. It also shows that it is not so straightforward to use the rank as a parameter to simplify our numerical decomposition. Instead, our numerical algorithm is effective to decompose arbitrary input superchannels regardless of their ranks.

IV. CONCLUSION
In this work, we developed the quantum circuit simulation algorithm of superchannels. It contains a classical optimization algorithm that designs the quantum circuits for the simulation of superchannels. This serves as a first study of superchannels and extreme superchannels from the viewpoint of quantum simulation and algorithm, laying the foundation for the near-term experimental implementation. Our algorithm applies to generic superchannels, and the convex-decomposition in terms of extreme ones proves to be effective, and can reduce the circuit cost compared with a direct dilation approach. The classical circuit-design algorithm may become intractable when the size of input channels gets larger, and it remains to see if better algorithms can be developed for these cases or cases of particular forms of input. We also studied the roles of special features or conditions in our algorithm, including being unital, extreme, and the rank of input superchannels. More specific ones can be studied case by case. For instance, there are superchannels that do not need a quantum ancilla, such as random-unitary superchannels, Pauli or depolarizing superchannels. The unitary twirling operation is such an example, which plays central roles for random benchmarking and error mitigation technique. Other analog of channels can be defined for superchannels, such as dephasing [58] and entanglement-breaking ones [59]. Their circuit simulations are quite straightforward given their concise forms.
We have seen that a key feature of superchannels is the memory (between the pre and post unitary operators). This feature actually has already been employed, even before the emergence of superchannels, for quantum teleportation and entanglement-assisted quantum communication [60]. Recently, we also pointed out the potential of using superchannels to design quantum superalgorithms and applications in quantum von Neumann architecture [30,51]. Given the growing importance of superchannels for quantum information science, it is therefore crucial to develop closer physical and al-gorithmic study of them. We believe our study can prompt the development of superchannels and their applications.
We emphasized that unital channels do not create nontrivial states from 1. Here we pursue the physical meaning of being extreme. So far we merely know that extreme points do not allow classical mixing, and the linearly independence condition is quite abstract. Below we study them from the viewpoint of complementary channels.
The complementary channel E c of a channel E is defined as If we treat ρ as the original quantum information, and it is in general impossible to recover it if we know E(ρ) and E. Only special types of channels can be recovered, and this actually leads to the error-correction condition tr(K † j K i ρ) = c ij [61]. However, if the final ancilla state σ f = E c (ρ) is known, then, together with the channel E, they can be used to recover the original state ρ if E is a full-rank extreme channel since {K i K † j } is a basis. That is, the ancilla aims to obtain the quantum information in ρ through a full-rank extreme channel. If the extreme channel is not of full rank, then only a partial information of ρ can be recovered. This agrees with an early scheme of using qubit extreme channels as an optimal asymmetric cloning machine [62], and also inspires potential extension of error-correction theory using extreme channels that shall be investigated separately.
There are extreme channels that are unital. These channels can be simulated by affine-mixing of unitaries, hence do not require a quantum ancilla. Therefore, we identify non-unital extreme (NUE) channels as genuinely quantum since they cannot be further decomposed. The amplitude-damping channel is such an example, which under composition will map any initial state into a fixed pure state, the ground state with the lowest energy. It has already been recognized as a resource for dissipative quantum state engineering and computing [13].
We expect that our understanding of extreme channels, especially non-unital ones, also applies to extreme superchannels. The complementary superchannel acts asŜ which extends the case for channels. If the set {S † b S a } is linearly independent and full-rank, then the original channel E can be recovered given the stateŜ c (ω E ) and the set {S a }. An additional feature is the memory wire between the pre V and post W , which forbids the reduction of a superchannel into two independent channels.