Efficient classical simulation of noisy random quantum circuits in one dimension

Understanding the computational power of noisy intermediate-scale quantum (NISQ) devices is of both fundamental and practical importance to quantum information science. Here, we address the question of whether error-uncorrected noisy quantum computers can provide computational advantage over classical computers. Specifically, we study noisy random circuit sampling in one dimension (or 1D noisy RCS) as a simple model for exploring the effects of noise on the computational power of a noisy quantum device. In particular, we simulate the real-time dynamics of 1D noisy random quantum circuits via matrix product operators (MPOs) and characterize the computational power of the 1D noisy quantum system by using a metric we call MPO entanglement entropy. The latter metric is chosen because it determines the cost of classical MPO simulation. We numerically demonstrate that for the two-qubit gate error rates we considered, there exists a length scale above which adding more qubits does not bring about an exponential growth of the cost of classical MPO simulation of 1D noisy systems. Specifically, we show that above the length scale, there is an optimal circuit depth, independent of the system size, where the MPO entanglement entropy is maximized. Most importantly, the maximum achievable MPO entanglement entropy is bounded by a constant that depends only on the gate error rate, not on the system size. We also provide a heuristic analysis to get the scaling of the maximum achievable MPO entanglement entropy as a function of the gate error rate. The obtained scaling suggests that although the cost of MPO simulation does not increase exponentially in the system size above a certain length scale, it does increase exponentially as the gate error rate decreases, possibly making classical simulation practically not feasible even with a state-of-the-art supercomputer.


I. INTRODUCTION
Quantum computers can provide significant computational advantage over classical computers as they can efficiently solve certain important problems that are believed to be not solvable in polynomial time with classical computers. Examples of such problems include integer factorization [1] and simulation of the real-time dynamics of large quantum systems [2]. While currently available quantum devices are not large and reliable enough to factor a large integer or simulate the dynamics of a large quantum system, it has been established over the past two decades that fault-tolerant quantum computing is in principle possible via quantum error correction [3][4][5][6][7]. However, despite the recent progress in reducing high resource overhead associated with the use of faulttolerant quantum computing schemes [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], large-scale and fault-tolerant quantum computing is not yet within reach of near-term quantum technologies.
Due to the lack of fault-tolerance, currently available noisy intermediate-scale quantum (NISQ) [24] devices are clearly not capable of realizing the full potential of quantum computing. Nevertheless, NISQ devices may be able to provide computational advantage over the best available classical computer in tackling certain computational tasks, whether or not solving them is practically useful. Various proposals for demonstrating such quantum computational advantage with NISQ devices have focused on sampling problems such as IQP [25], boson sampling [26,27], Fourier sampling [28], * noh827@gmail.com † liang.jiang@uchicago.edu ‡ wjf@uchicago.edu and random circuit sampling (RCS) [29]. In particular, various complexity-theoretic hardness results for these sampling problems have made them an appealing proposal for demonstrating quantum computational advantage.
Among these proposals, boson sampling was the first whose hardness was proven to be robust against adversarial total variation distance noise, under reasonable hardness assumptions from computational complexity theory [26]. Motivated by such a rigorous hardness result, experimental realizations of boson sampling followed shortly thereafter [30][31][32][33]. However, all the previous boson sampling experiments have been performed with a limited number of photons that is not large enough to make the system classically intractable. Moreover, various high-performing classical algorithms that are tailored to boson sampling have been developed [34][35][36][37]. These recent developments have thus made it much more challenging to demonstrate quantum computational advantage via boson sampling.
Over the past few years, RCS has risen as a promising candidate for achieving quantum computational advantage since it can be realized at scale in superconducting qubit systems [29]. While initially motivated by the experimental viability, RCS was also recently shown to have similar asymptotic hardness guarantees as boson sampling [38] and complementary hardness evidence was shown in Ref. [39], making RCS an even more compelling proposal. Notably, RCS was implemented in a superconducting qubit system which astonishingly consists of 53 qubits that are connected via two-qubit gates with very low gate error rates (p ∼ 0.006) in a planar architecture [40]. In particular, it was claimed in Ref. [40] that it would take about 10000 years for a state-of-the-art classical supercomputer to achieve a computational task that is equivalent to the one that their superconducting quantum device has arXiv:2003.13163v1 [quant-ph] 29 Mar 2020 Non-trivial quantum correlation

Shallow circuits
Deep circuits

Maximum achievable quantum correlation (focus of our work)
FIG. 1: Schematic plot of the degree of non-trivial quantum correlation as a function of the circuit depth. When the circuit depth is small, quantum correlation grows linearly in the circuit depth. On the other hand, when the circuit depth is large, the system converges to a depolarized state and thus the non-trivial quantum correlations are washed away. The focus of our work is to understand the optimal regime where the maximum non-trivial quantum correlation is achieved. See also Figs. 5 and 9.
achieved. In contrast, a recent work [41] has suggested that a refined simulation technique can bring down the required computing time to just a few days. In any case, what has become clear is that currently available superconducting qubit systems can tackle certain computational tasks that lie close to the borderline of what is achievable and not achievable with classical computing technologies. Going forward, an important thing to keep in mind is that currently available quantum devices are noisy. Thus, a crucial related question is how the classical computing time needed to simulate such noisy random quantum circuits would scale as a function of the system size and the gate error rate. Thanks to the rigorous complexity-theoretic results [25,26,38,42,43], it has been established that in the noiseless case, simulating the outputs of a random quantum circuit cannot be done classically in polynomial time in the system size. Moreover, for boson sampling and RCS, the classical intractability was shown to persist even in the presence of the total variation distance noise under suitable hardness assumptions [26,38,42,43].
However, modeling noise using only closeness in total variation distance does not suffice to address practically relevant settings such as the setting in which each gate is corrupted by an error channel with a non-zero gate error rate. This is because in realistic settings, the effects of noise dominate in the large circuit depth limit and the system eventually converges to a depolarized state. Thus, it is not immediately clear how much computational power can be gained by adding more qubits to noisy quantum systems. Addressing this question is thus essential for understanding the utility of near-term applications of NISQ technologies.
In this paper, we study noisy random circuit sampling in one dimension (i.e., 1D noisy RCS) as a simple model for exploring the effects of noise on the computational power of a noisy quantum device. Note that since noisy systems eventu-ally converge to a depolarized state, any non-trivial quantum correlations will be washed away in the large circuit depth limit, making the outputs of the system well approximated by a trivial uniform distribution. On the other hand, shallow circuits with a constant circuit depth can also be simulated efficiently via matrix product states (MPSs) [44] due to the limited growth of entanglement. Connecting these two extreme cases, we can expect that the degree of non-trivial quantum correlation will be peaked at a certain optimal circuit depth as illustrated in Fig. 1. The focus of our work is to understand this optimal regime where the maximum non-trivial quantum correlation is attained. In particular, we explore how hard it is to simulate the 1D noisy system at the optimal circuit depth.
Note that if one's goal is to approximately simulate ideal random quantum circuits with any non-zero fidelity, it may suffice to use MPSs with a constant bond dimension even for deep circuits (see, e.g., Ref. [45] and also the discussion in Section V). However, this is not our goal here and we instead aim to simulate noisy random quantum circuits to any desired accuracy. More specifically, we directly simulate the mixed state dynamics of 1D noisy random quantum circuits by using matrix product operators (MPOs) [46,47]. Note that this is a strictly more challenging task than sampling since any output probability (including marginal and conditional probabilities) can be computed efficiently from an MPO and thus sampling can be done straightforwardly.
The main contribution of our work is to characterize the computational power of 1D noisy quantum devices by using a metric we call MPO entanglement entropy. We choose the latter metric because it determines the cost of classical MPO simulation as well as the degree of non-trivial quantum correlation of a mixed state. We numerically demonstrate the maximum achievable MPO entanglement entropy is bounded by a constant that depends only on the gate error rate, not on the system size. In other words, the maximum achievable MPO entanglement entropy is saturated at a certain length scale and consequently the required MPO bond dimension does not increase exponentially in the system size above the length scale. Thus, our results indicate that there exists a length scale above which adding more qubits does not help increasing the cost of classically simulating a 1D noisy quantum device in an exponential way. We also provide a heuristic argument to get the scaling of the maximum achievable MPO entanglement entropy as a function of the gate error rate. The obtained scaling suggests that the cost of MPO simulation increases exponentially as the gate error rate decreases, possibly making classical simulation practically not feasible even with a stateof-the-art supercomputer.
Our paper is organized as follows. In Section II, we formulate the problem of 1D noisy RCS. In Section III, we briefly review matrix product operators (MPOs) and define MPO entanglement entropy that determines the cost of classical MPO simulation. In Section IV, we present the main numerical results on the MPO simulation of 1D noisy RCS. Moreover, in Subsection IV B, we provide a heuristic scaling analysis of the maximum achievable MPO entanglement entropy as a function of the gate error rate. In Section V, we discuss the relation of our results to previous results. We conclude the paper FIG. 2: Noisy random circuit sampling in one dimension. Each noisy two-qubit gate is given by a 4 × 4 Haar-random unitary operation followed by a two-qubit depolarization channel N 2 [p] with an error rate p. At the end of the circuit, all the qubits are measured in the computational basis. For simplicity, we only consider even number of qubits. Although the maximum circuit depth D is chosen to be even in the schematic illustration, we allow D to be odd as well.
by outlining several open questions in Section VI.

II. PROBLEM SETUP
In this section, we formulate the problem of noisy random circuit sampling in one dimension (i.e., 1D noisy RCS). We also introduce two-qubit depolarization error model which we assume to get the numerical results in Section IV.

A. Noisy random circuit sampling in one dimension
Consider n qubits laid out in a one-dimensional chain as shown in Fig. 2. For simplicity, we only consider even n. Initially, all the qubits are prepared in the computational zero state, i.e., In odd (or even) time steps, Haar-random two-qubit gates are applied to the l th and the l+1 th qubits for l ∈ {1, 3, · · · , n−1} (or l ∈ {2, 4, · · · , n − 2}). Eventually, we will assume that each Haar-random two-qubit gate is corrupted by a noisy twoqubit CPTP map [48] acting on the same sites. However, we assume that two-qubit gates are noiseless for now. The circuit depth D is defined as the number of time steps. Although only the case with an even D is shown in Fig. 2, D can also be odd. At the end of the circuit, all the qubits are measured in the computational basis {|0 , |1 }. Thus, we are left with an output n-bit string which is drawn from a probability distribution Here,Û C is the unitary operator associated with an instance of a depth-D random circuit C. Exact RCS is a sampling problem where the goal is to sample exactly from the ideal output distribution P C of a given quantum circuit C. One can also define an approximate version of RCS, i.e., approximate RCS, which is a sampling problem where the goal is to sample from a distribution P C that is -close to the ideal output distribution P C ) in total variation distance, i.e., In the ideal setting with noiseless two-qubit gates, it has been established that approximate RCS is hard in the average case [38,39,42,43]. In particular, the approximate hardness implies that the classical intractability of RCS persists even in the presence of the total variation distance noise given in Eq. (4). On the other hand, it is important to realize that realistic quantum devices are not able to sample, even approximately, from an ideal output distribution P C in the limit of large system size. In particular, all the gates in realistic quantum devices are corrupted by an error channel with a non-zero gate error rate and thus the noisy system eventually converges to a depolarized state. As a result, the fidelity between an actual output state obtained from a noisy quantum circuit C and an ideal output state obtained from a noiseless quantum circuit C decreases exponentially in the system size [40].
For these reasons, the adversarial total variation distance noise model is not directly relevant to realistic settings. Therefore, it is important to investigate noisy versions of RCS where each two-qubit gate is corrupted a noisy two-qubit CPTP map N . Specifically, we numerically study noisy RCS in 1D architecture by using matrix product operators (MPOs). By doing so, we explore the effects of noise on the computational power of a 1D noisy quantum device.

B. Noise model: Two-qubit depolarization channel
To make the discussion concrete, we assume that the noise map N is given by a two-qubit depolarization channel N 2 [p] with an error rate p. The two-qubit depolarization channel N 2 [p] is defined as where E 2 ≡ {Î,X,Ŷ ,Ẑ} ⊗2 − {Î ⊗Î} is the set of 15 nontrivial two-qubit Pauli operators and {Î,X,Ŷ ,Ẑ} is the set of single-qubit Pauli operators. Thus, we incorporate the possibility of correlated two-qubit errors that occur during a twoqubit gate. Note also that the two-qubit depolarization channel N 2 [p] can also be expressed as whereÎ ⊗Î/4 is the maximally mixed (or completely depolarized) two-qubit state.
In the context of fault-tolerant quantum computing, twoqubit depolarization channels are used to model errors that happen during two-qubit gates, such as CNOT and CZ gates (i.e., to perform a detailed circuit-level noise analysis) [7,8,15,[18][19][20][21][22][23]. For instance, the fault-tolerance threshold of the surface-code p 0.01 [7] is obtained by applying the two-qubit depolarization channel N 2 [p] after each two-qubit gate. Since any two-qubit errors can be converted via a noise twirling [49] to a two-qubit depolarization channel, the latter serves as a pessimistic noise model for a detailed circuit-level fault-tolerance anaylsis.
In the context of noisy RCS, we consider the two-qubit depolarization model simply as a representative error model. We expect that different error models will give rise to the same conclusions we reach with the two-qubit depolarization model. Moreover, the MPO method is completely general and applies to any two-qubit error model.

III. MATRIX PRODUCT OPERATORS
In this section, we briefly review matrix product operators (MPOs) [46,47] which we use to analyze 1D noisy RCS. We also introduce MPO entanglement entropy that determines the cost of classical simulation based on MPOs. While we focus on qubits in later sections, we consider qudits here to make the review as general as possible.

B. Canonical update of MPOs
In 1D noisy RCS, we keep applying two-qudit CPTP maps on two neighboring qudits, say the l th and the l + 1 th qudits. Thus, to simulate 1D noisy RCS using MPOs, it is important to update an MPO in a canonical way after applying each twoqudit CPTP map, so that we are guaranteed to be left with an updated MPO in the canonical form. In the case of matrix product states (MPSs) that are pure states, a canonical update upon the action of a local unitary operator can be done locally [44]. For instance, if a single-qudit unitary operationÛ is applied to the l th qudit, one can simply update Γ where U i l j l ≡ i l |Û |j l . All the other parameters need not be updated because single-qudit unitary operations cannot affect the entanglement structure of the chain. Note that we used i l instead of I l = di l +ī l since we are dealing with pure states when considering MPSs. Similarly, if a two-qudit unitary operation is applied to the l th and the l + 1 th qudits, only need to be updated because the unitary operation can only affect the entanglement along the cut [1 · · · l] : [l + 1 · · · n].
In contrast, when we work with MPOs and CPTP maps, canonical update cannot be done locally any more because CPTP maps are generally not unitary. To elaborate more on this, let us get back to the case of pure states and consider an n-qubit GHZ state Note that in this case, the bond dimension χ = 2 suffices (i.e., α l ∈ {0, 1} for all l ∈ {1, · · · , n − 1}) and the canonical λ parameters (or singular values, defined similarly as in the case of MPOs) are given by This is because the Schmidt decomposition of the GHZ state is given by for any cut [1 · · · l] : [(l + 1) · · · n]. Then, as an example of local non-unitary action on the system, let us consider a nondestructive measurement of the l th qubit in the computational basis {|0 , |1 }. Note that we would get either |0 or |1 , each with 50% probability, as a measurement outcome. Then, conditioned on measuring |0 or |1 , the n-qubit GHZ state collapses to a product state |0 ⊗n or |1 ⊗n . Thus in any case, the updated singular values of the post-measurement state are given by This example clearly illustrates that even a local action can make a global impact on the entanglement structure of the chain when the action is not unitary. Note that in this example, the bond dimension needed to describe the output product output state is given by χ = 1, instead of χ = 2. Such a reduction in the required bond dimension is thanks to the fact that the non-destructive measurement decreased non-trivial quantum correlations between disjoint subregions in the chain that were present in the initial GHZ state.
With this observation in mind, let us now consider the case of mixed states (or MPOs) and two-qudit CPTP maps. Recall that in the 1D noisy RCS, we apply noisy Haar-random two-qubit gates that are corrupted by a CPTP noise map N (e.g., a two-qubit depolarization channel N 2 [p]). That is, each time we apply a noisy Haar-random two-qubit gate, we apply a two-qubit CPTP map where U is defined as U(ρ) ≡ÛρÛ † andÛ is a 4 × 4 Haarrandom unitary operator. Since the system starts from a completely uncorrelated product state |0 ⊗n , Haar-random unitary operations will initially make the system more correlated in a quantum way, and thus increase the required bond dimension needed to faithfully describe the system. On the other hand, the noise map N will generally tend to decrease non-trivial quantum correlations between disjoint regions across the entire chain, reducing the required bond dimension. In particular, we can expect that the effects of noise will eventually take over and the initially developed non-trivial quantum correlations will be washed away as the circuit depth increases. All these expected behaviors will be corroborated numerically in Section IV (see Figs. 5 and 9). To do so, we need to update MPOs such that they remain in the canonical form after apply-ing each noisy two-qubit gate. How this is done is explained in a great detail in Appendix A.

C. MPO entanglement entropy
We show in Appendix B that the time cost of MPO simulation of 1D noisy RCS (using the update method described in Appendix A) is given by Here, n is the number of qubits, D is the circuit depth, and χ is the bond dimension. Thus, the simulation cost is determined by the bond dimension χ for a given set of the system size n and the circuit depth D. A relevant quantity that determines the required bond dimension χ is the spectrum the singular values λ [l] α l of an MPO in the canonical form. In the case of pure states and MPSs, the singular values are directly related to the known entanglement measures of pure states such as the entanglement entropy [50].
In the case of mixed states and MPOs, although the singular values λ α l are not directly related to known entanglement measures for mixed states [51], they can still be used to characterize the degree of quantum correlations between two disjoint regions [1 · · · l] : [(l + 1) · · · n]. Specifically, we consider the following quantity, which we call MPO entanglement entropy, to measure the degree of quantum correlations in an MPO: Here, l ∈ {1, · · · , n − 1} and λ [l] α l are the singular values of an MPO |ρ in the canonical form. Note that we used the symbol S instead of S to distinguish the MPO entanglement entropy from the usual entanglement entropy. We also normalized the spectrum of the squared singular values (λ [l] α l ) 2 so that they sum up to unity. We took the squared singular values because then the MPO entanglement entropy is reduced to twice the usual entanglement entropy in the case of pure states, i.e., Here, S(ρ) = −Tr[ρ log 2ρ ] is the von Neumann entropy and is the reduced density matrix of the state |ψ with respect to the subsystem [1 · · · l]. The additional factor of 2 is due to the fact that density matrices are composed of both kets and bras whereas states vectors are described by kets only. More specifically, the MPO singular values of a pure state are simply given by the products of the two copies of the MPS singular values. We choose not to divide the MPO entanglement entropy by 2 because this additional factor reflects the fact that simulating density matrices is computationally more costly than simulating state vectors. We emphasize that for mixed states, the MPO entanglement entropy S l (ρ) does not correlate with the von Neu-mann entropy of the reduced density matrix S(ρ . For instance, for the completely and globally depolarized stateρ =Î ⊗n /2 n , the von Neumann entropy of the reduced density matrixρ [1···l] = 2 ⊗l /2 l is given by S(ρ [1···l] ) = l and thus grows extensively with the system size l. However, the completely and globally depolarized state does not have any non-trivial quantum nor classical correlations between disjoint subsystems. That is, the extensive von Neumann entropy only quantifies the extensive noise in the depolarized state, not its correlation. On the other hand, the MPO entanglement entropy of the completely and globally depolarized state vanishes for all subsystem size l ∈ {1, · · · , n − 1}, i.e., S l (ρ =Î ⊗n /2 n ) = 0. In this sense, the MPO entanglement entropy captures the non-trivial quantum correlations present in a mixed state, separating out the effects of noise.
In addition to quantifying the degree of non-trivial quantum correlation, the MPO entanglement entropy also determines the cost of classical MPO simulation. In particular, we define the maximum MPO entanglement entropy as follows The maximum MPO entanglement entropy can be used to estimate the required bond dimension χ needed to describe a mixed stateρ. Specifically, if the bond dimension χ satisfies taking only the largest χ singular values and discarding all the smaller ones will have a negligible effect on the accuracy of an MPO simulation. The cost of MPO simulation depends heavily on the required bond dimension χ (see Appendices A and B). Thus, the maximum MPO entanglement entropy can serve as a metric that characterizes the computational power of a 1D noisy system. In Section IV, we numerically demonstrate that in the 1D noisy RCS setting, the maximum MPO entanglement entropy is bounded by a constant independent of the system size (or the number of qubits n), implying that 1D noisy random circuits can be simulated efficiently using MPOs.

D. Efficient computation of the output probabilities
Recall that in RCS, all the qubits are measured in the computational basis at the end of the circuit. Here, we explain how the output probabilities can be efficiently computed from an MPO. Given a mixed stateρ, the probability of getting an output bit string x = (x 1 , · · · , x n ) from a computational-basis measurement is given by In the vectorized representation, the state | x x| is mapped to Thus, Pρ( x) is given by Plugging in the canonical MPO representation of a mixed statê ρ in Eq. (11), we find Note that this quantity can be efficiently computed via a sequential matrix multiplications of one 1 × χ matrix, 2n − 3 χ × χ matrices, and one χ × 1 matrix. Also, we can similarly compute the total probability Tr[ρ] as follows: Ideally, the total probability Tr[ρ] should be unity but it will be smaller due to the discarded small singular values. In our numerical demonstration below, we compute the total probability Tr[ρ] to characterize the truncation errors. Sampling from the output probability distribution can also be done efficiently. To do so, we first derive the marginal probability distribution of the first dit x 1 by computing for x 1 ∈ {0, · · · , d−1} and get a samplex 1 from the marginal distribution P [1] ρ (x 1 ). Then, we compute the conditional probability given the samplex 1 , i.e., where P and then get a samplex 2 from the conditional distribution P [2|1] ρ (x 2 |x 1 ). One can then sequentially get the remaining sample ditsx 3 , · · · ,x n from the conditional probabilities P which can be computed similarly as above.

IV. MAIN RESULTS
In this section, we present our main results on the MPO simulation of 1D noisy random quantum circuits. Note that in noisy RCS, it suffices to sample an output bit string x = (x 1 , · · · ,x n ) from the output distribution of a noisy quantum circuit. However, we aim to achieve a strictly more challenging task. That is, we will directly simulate the density matrix of the systemρ in real time using MPOs. As detailed in the previous section, sampling from the output distribution can be straightforwardly done if the output stateρ is available.
In the first time step (or circuit depth 1), we generate n/2 Haar-random two-qubit unitary operatorŝ by sampling n/2 Haar-random unitary matrices of size 4 × 4.
Then, we construct corresponding n/2 two-qubit CPTP maps by applying the two-qubit depolarization channel N 2 [p] with a two-qubit gate error rate p, i.e., where U are defined similarly as in the text below Eq. (21). We sequentially apply these two-qubit CPTP maps and update the MPO parameters in a canonical way as prescribed in Appendix A. After applying all n/2 two-qubit CPTP maps, we compute the MPO entanglement entropies S l (|ρ ) for all l ∈ {1, · · · , n − 1} and save them.
Similarly in the second time step (or circuit depth 2), we generate (n − 2)/2 Haar-random two-qubit unitary operatorŝ and sequentially apply the corresponding CPTP maps M and update the MPO parameters as described in Appendix A. Then, we compute and save the MPO entanglement entropy S l (|ρ ) for all l ∈ {1, · · · , n−1}.
We keep alternating between these two procedures and save the MPO entanglement entropies at the end of each time step. Since all the two-qubit unitary operators are chosen Haarrandomly, the obtained MPO entanglement entropies will vary across different circuit realizations. Thus, we repeat this entire procedure N s times and take the average of the obtained entanglement entropies over N s circuit realizations. Specifically, at each circuit depth and length of the subsystem l, we first average the MPO entanglement entropy S l (ρ) over N s circuit realizations. Then, we maximize the averaged MPO entanglement entropy over l to compute the maximum MPO entanglement entropy S max (ρ). In all the numerical simulations presented below, we choose N s = 24.
In Fig. 3, we consider the cases with n = 8 qubits subject to various two-qubit gate error rates 0 ≤ p ≤ 0.1. In particular, we choose the bond dimension to be χ = 2 8 = 256 so that there are no errors in the MPO simulation due to the bond dimension truncation. In the case of noiseless two-qubit gates (i.e., p = 0), the maximum MPO entanglement entropy is achieved in the middle of the chain (i.e., along the cut [1 · · · 4] : [5 · · · 8] or l = 4) in the large circuit depth limit and converges to S max = 6.56.
Note that this value is the same as twice the average entanglement entropy of the 8-qubit Haar-random states along the cut [1 · · · 4] : [5 · · · 8], which is given by Here, we used the formula established in Refs. [52][53][54][55], i.e., and [5 · · · 8]) consists of 4 qubits and thus 16 states. Such an agreement is consistent with the expectation that the system converges to a Haar-random state in the large circuit depth limit if all the two-qubit gates are noiseless.
On the other hand, if the two-qubit gates are noisy (i.e., p = 0), the system eventually converges to the completely and globally depolarized stateÎ ⊗8 /2 8 . As can be seen from Fig.  3, for any p = 0, the maximum MPO entanglement entropy indeed decreases exponentially as the circuit depth increases. This is consistent with the fact that the system converges to the completely and globally depolarized state, which does not possess any non-trivial quantum nor classical correlations.
In Fig. 4, to get more intuition on what happens in the n = 8 qubit cases considered in Fig. 3, we plot the output probability distribution P ( x) for all possible 2 8 = 256 outputs x ∈ {0, 1} 8 for a specific random circuit realization. Note that we ordered the outputs such that a heavier output with larger probability has a smaller label than that of a lighter output. In Fig. 4(a), we consider the noiseless case with the vanishing gate error rate p = 0. In this case, as discussed above, the system converges to the 8-qubit Haar-random states in the large circuit depth limit. Correspondingly, the sorted output probability distribution also converges to a fixed distribution in the large circuit depth limit. Note that the output distribution in the large circuit depth limit is far from being uniform. In other words, there are heavier outputs with larger probability that occur more often and lighter outputs that occur less frequently.
In Fig. 4(b), on the other hand, we consider noisy two-qubit gates with a gate error rate p = 0.1. In the noisy case, the system eventually converges to the completely and globally depolarized state in the large circuit depth limit. Indeed, as shown in Fig. 4(b), the output probability distribution converges to the trivial uniform distribution in the large circuit depth limit.
In Fig. 5, we consider the cases with a fixed two-qubit gate error rate p and vary the number of qubits n. In all the cases we consider, we observe that the MPO entanglement entropy is maximized at a certain optimal circuit depth independent of the system size n. Moreover, the maximum achievable MPO entanglement entropy at the optimal circuit depth does not depend on the system size n. Consequently, the minimum bond dimension χ needed to capture the majority of the total probability does not increase exponentially in the number of qubits n in the limit of large n.
For example, when the two-qubit gate error rate is given by p = 0.15 (see Fig. 5(a)), the MPO entanglement entropy is maximized at an optimal circuit depth D = 4 for all n ≥ 8. Also, the maximum achievable MPO entanglement at the optimal circuit depth is given by S max 2.75, which is independent of the system size n as long as n ≥ 8. Most importantly,    we observe that having more than 8 qubits does not help increasing the MPO entanglement entropy. As a result, the chosen bond dimension χ = 80 is at least ten times larger than 2 S max 6.7 and therefore is large enough to reliably describe the noisy system for any system size 4 ≤ n ≤ 18 we considered. Indeed, we numerically confirm that we accounted for at least 99.5% of the total probability (i.e., Tr[ρ] ≥ 0.995) on average with the bond dimension χ = 80.

Circuit depth
For a smaller gate error rate p = 0.1 (see Fig. 5(b)), the optimal circuit depth that maximizes the MPO entanglement entropy is given by D = 6 for all n ≥ 8. Also, the corresponding MPO entanglement entropy is given by S max 4 and is independent of the system size above the length scale n = 8. The chosen bond dimension χ = 150 is again about ten times larger than 2 S max 16 and we accounted for at least 99.7% of the total probability on average. Note that higher bond dimension is required in this case than in the case of p = 0.15, because higher MPO entanglement entropy can be achieved. Note also that the saturation of MPO entanglement entropy occurs with n = 8 qubits at the optimal circuit depth D = 6. As shown in Fig. 4(b) (see the red line), the sorted output probability distribution at this optimal circuit depth is far from being uniform and thus is still non-trivial. For an even smaller gate error rate p = 0.06 (see Fig. 5(c), we observe that the optimal circuit depth is given by D = 9 or D = 10 for all n ≥ 12 and the maximum achievable MPO entanglement entropy is given by S max 6.3. In this case, the chosen bond dimension χ = 400 is about five times larger than 2 S max 79 and the accounted total probability is at least 99.1% on average.
In Fig. 6, to get more understanding of why the saturation of the MPO entanglement entropy happens, we take the twoqubit gate error rate p = 0.1 (which was considered in Fig.  5(b)) and consider n = 32 qubits. We remark that simulating the 32-qubit system is not so costly because the constant bond dimension χ = 150 suffices even for 32 qubits. In particular, we zoom in to see a more fine-grained MPO entanglement structure and plot the MPO entanglement entropy S l as a function of the length of the subsystem l (i.e., with respect to the cut [1 · · · l] : [(l + 1) · · · n]) for various circuit depths. Note that the largest MPO entanglement entropy is achieved at an optimal circuit depth D = 6, which is consistent with the observation in Fig. 5(b). Notably, we can see that at the optimal circuit depth D = 6 (and far away from the boundaries, i.e., 4 ≤ l ≤ 28), the MPO entanglement entropy S l is independent of the subsystem size l. In other words, the MPO entanglement entropy follows an area law at the optimal circuit depth D = 6. This is because the optimal circuit depth is bounded by a constant independent of the system size due to noise and consequently qubits that are not contained within a finite causal cone cannot be correlated with the ones that lie in α l for the n = 32 qubits and the two-qubit gate error rate p = 0.1 at the optimal circuit depth D = 6. We took the subsystem size l = 14 where the maximum MPO entanglement entropy is achieved (see Fig.  6). Among the N s = 24 circuit instances, we only show two extreme instances with the largest and the smallest MPO entanglement entropy. the causal cone. We provide more discussions on the interplay between noise and the circuit depth in the following section.
Lastly in Fig. 7, to see the effects of the bond dimension truncation, we plot the spectrum of the MPO singular values λ [l] α l for the case with n = 32 qubit and the two-qubit gate error rate p = 0.1 considered in Fig. 6. In particular, we took the optimal circuit depth D = 6 and the subsystem size l = 14 where the MPO entanglement entropy is maximized. Note that the latter choice is somewhat arbitrary since the MPO entanglement entropy is nearly constant deep inside the bulk (i.e., for 4 ≤ l ≤ 28). Among the N s = 24 circuit realizations, we show the two extreme instances with the largest and the smallest MPO entanglement entropy. In any cases, we can clearly see from α l decrease exponentially as the bond index α l increases. Thus, the bond dimension truncation has a negligible effect as long as log 2 χ is much larger than the maximum achievable MPO entanglement entropy S max . In particular, due to the exponential decay of the singular values, the required bond dimension χ as well as the time cost of the MPO simulation would increase only poly-logarithmically in 1/ , where is the simulation error measured in total variation distance similarly as in Eq. (4) (see, e.g., Ref. [56]).

B. Maximum achievable MPO entanglement entropy
Recall that the time cost of MPO simulation of 1D noisy RCS is given by T = O(n 2 Dχ 3 ) (see Appendix B). Thus, the classical MPO simulation cost depends heavily on how large the bond dimension χ needs to be. As demonstrated above, it suffices in practice to choose the bond dimension χ such that for some constant c 1 (see also Eq. (26)). In Fig. 5, for instance, we chose at least c ≥ 5 and were able to capture more than 99.1% of the total probability on average. Combining the facts that the simulation cost increases cubically in the bond dimension χ, and the required bond dimension χ is proportional to 2 S max , we can infer that the simulation cost increases exponentially in the maximum achievable MPO entanglement entropy S max . Thus, the maximum achievable MPO entanglement entropy S max can be used as a measure for characterizing the computational power of a 1D noisy quantum system.
In Fig. 8, we plot the maximum achievable MPO entanglement entropy S max as a function of the number of qubits n for various two-qubit gate error rates 0.06 ≤ p ≤ 0.15. In all the cases, we observe that S max increases linearly in the system size n when the system size is small. In this case, the cost of MPO simulation increases exponentially in the number of qubits n. In other words, the computational power of a 1D noisy system increases exponentially in the system size n when n is small. However, as the system size n increases further, the maximum achievable MPO entanglement entropy S max converges to a constant value S max,∞ which is independent of the system size n. For p = 0.06, for example, S max converges to S max,∞ 6.3. In particular, the saturation occurs around the system size n = 12. This implies that for p = 0.06, once we have n = 12 qubits, adding more qubits does not bring about an exponential growth of the computational power because MPO entanglement entropy remains unchanged and thus the required bond dimension χ does not increase exponentially in the system size. Indeed, a constant bond dimension χ = 400 was sufficient to capture the majority (more than 99.1%) of the total probability for all n ≤ 18 in the p = 0.06 case. Consequently, above the length scale n = 12, the MPO simulation cost increases only quadratically, not exponentially, in the system size n since T = O(n 2 Dχ 3 ).
Note that as the two-qubit gate error rate p becomes larger, the saturation of the MPO entanglement entropy happens at a smaller system size and the saturated value S max,∞ becomes smaller. For example, when p = 0.1, the saturation occurs around the system size n = 8 and the corresponding MPO entanglement entropy is given by S max,∞ 4. In this case, the bond dimension χ = 150 is sufficient to capture more than 99.7% of the total probability on average for all n ≤ 18. Note that we need χ = 2 8 = 256 to exactly describe an 8-qubit system. Since a smaller bond dimension χ = 150 suffices for p = 0.1, we can infer that a 1D noisy system with a two-qubit gate error rate p = 0.1 cannot fully occupy the entire 8-qubit Hilbert space (see also Figs. 3 and 4). As a result, having more than 8 qubits does not bring about an exponential growth of the classical simulation cost.
We have demonstrated in Figs. 5 and 8 that there exists a length scale, determined by the gate error rate p, below which the classical MPO simulation cost increases exponentially in the system size n, but above which does so only polynomially in n. This is due to the saturation of the MPO entanglement entropy in the large system size limit. If the system size is large enough so that the MPO entanglement entropy is saturated, the saturated MPO entanglement entropy depends solely on the gate error rate p, i.e., S max,∞ = S max,∞ (p). Moreover, the latter becomes the key quantity that determines the required bond dimension χ and consequently the overall classical MPO simulation cost in the saturated regime.
Note that for the two-qubit gate error rates we considered (i.e., 0.06 ≤ p ≤ 0.15), the MPO entanglement entropy is sufficiently saturated with n = 16 qubits (see Fig. 8). In Fig.  9(a), we thus take a closer look into the n = 16 case and plot the maximum MPO entanglement entropy S max as a function of the circuit depth D for various two-qubit gate error rates. As can be seen from Fig. 9(a), the optimal circuit depth D (p) where the MPO entanglement entropy is maximized increases as the gate error rate p decreases. Moreover, the corresponding MPO entanglement entropy S max,∞ (p) increases as p decreases, as shown in Fig. 9(b). Numerically, we were not able to investigate the cases with low gate error rates (i.e., p ≤ 0.06) because the required bond dimension becomes larger than χ = 400 and thus the computing time gets too large. For comparison, note that with the bond dimension χ = 256 (χ = 512), we can perform an exact MPS simulation of any pure states of a 16-qubit (18qubit) system. Thus, instead of exploring the low gate error regime numerically, we provide a heuristic analysis of the optimal circuit depth D (p) and the maximum achievable MPO entanglement entropy S max,∞ (p) in the small p regime. Note that the MPO entanglement entropy is primarily determined by the qubits and gates that are contained within a causal cone. For a given circuit depth D, there are O(D) qubits and O(D 2 ) noisy two-qubit gates (with an error rate p per gate) that are enclosed in the causal cone in the bulk cases (e.g., l ∈ [n/4, 3n/4]). Thus, one can think of the following heuristic model for the MPO entanglement entropy S l : The saturated values of the maximum achievable MPO entanglement entropy S max,∞ (p) as a function of the two-qubit gate error rate p for 0.06 ≤ p ≤ 0.15. To find S max,∞ (p), we took the maximum MPO entanglement entropy at the optimal circuit depth D (p) for each p from the n = 16 data shown in (a). Note that for 0.06 ≤ p ≤ 0.15, the MPO entanglement entropy is saturated with n = 16 qubits (see Fig. 8).
is due to the linear growth of the MPO entanglement entropy in the absence of noise. The effects of noise are crudely modeled by the term (1−p) cD α p 1 − −− → De −cpD α . If all the O(D 2 ) gates contained in the causal contribute equally, the exponent α would be given by α = 2.
Since ∂ D S l (D, p) = (1 − αcpD α )e −cpD α = 0 implies αcpD α = 1, the optimal circuit depth D (p) for a given error rate p is given by Thus, we also find The scaling in Eq. (45) implies that the saturation of the MPO entanglement entropy observed in Fig. 8 will hold for any nonzero gate error rate p > 0. Thus, our heuristic analysis suggests that there does not exist a non-trivial threshold value of the gate error rate p below which efficient classical simulation is not possible. This is because the required bond dimension χ 2 S max does not grow exponentially in the system size n since the MPO entanglement entropy is saturated.
On the other hand, the heuristic analysis also suggests that the bond dimension χ increases exponentially in (1/p) 1/α . Hence, the classical MPO simulation cost increases exponentially as the gate error rate p decreases, possibly makes classical simulation practically not feasible. This is precisely the reason why we were not able to explore smaller gate error rates than p = 0.06 in our numerical simulation. In this regards, we remark that the heuristic analysis works in the regime where the optimal circuit depth D is large or equivalently when the gate error rate p is small. On the other hand, the gate error rates we considered 0.06 ≤ p ≤ 0.15 are not sufficiently small for us to extract the exponent α in a stable manner. Such a quantitative fitting will be reliable only in the low gate error regime. Numerically investigating the low gate error regime would require more advanced computing resources and techniques. We thus leave it as a future research direction.

V. RELATION TO PREVIOUS RESULTS
In this section, we compare our results with the related previous results. In essence, our work differs from the previous ones in that we directly simulate mixed states via MPOs that are corrupted by CPTP noise maps. That is, we take advantage of the fact that noise maps (e.g., depolarization channels) reduce non-trivial quantum correlations between disjoint subsystems, making it possible for us to reliably describe the corrupted mixed states with a bounded bond dimension. Indeed, we numerically demonstrate that for 1D noisy random quantum circuits, the MPO entanglement entropy, which characterizes the non-trivial quantum correlations, is bounded by a constant independent of the system size n (see Fig. 8). Moreover, the heuristic analysis given in the previous section suggests that the same qualitative behavior will hold for any non-zero gate error rate p > 0 (see Eq. (45)). Hence, our work suggests that there does not exist a non-zero threshold value of the gate error rate p below which efficient classical simulation is forbidden. In what follows, we explicitly compare these features with the ones observed in the previous works.
Recently, there have been numerous studies on 1D random quantum circuits subject to local projective (or weak) measurements [57][58][59][60][61][62][63][64][65][66][67][68]. In this model, local projective (or weak) measurements (which happens with a probability p m at each qubit site in a given time step) reduce non-trivial quantum entanglement between disjoint subsystems. In our model, on the other hand, the same role is played by CPTP noise maps. In both models, random unitary operations tend to increase non-trivial quantum correlations, and thus compete with the correlation-decreasing measurements or CPTP noise maps. In the case of projective (or weak) measurements, the 1D system always remains in a pure state. Also, it has been numerically observed that there is a threshold value of the measurement probability (p (th) m 0.16) above which area-law entanglement holds. Hence, above the threshold (p m > p (th) m ), dynamics of the random circuits can be efficiently simulated by using MPSs. However, volume-law entanglement holds below the threshold (p m < p (th) m ). Thus, the system cannot be efficiently simulated by using classical computers unless the random unitary operators have a special structure such as being Clifford operations which can be efficiently simulated regardless of the entanglement structure [69][70][71].
In contrast, since we consider CPTP noise maps which are realistic models of noise in NISQ devices, the 1D system is in a mixed state and hence we use MPOs to efficiently simulate the system. In particular, CPTP noise maps completely wash away non-trivial quantum correlations and the system eventually reaches the completely and globally depolarized statê I ⊗n /2 n in the large circuit depth limit. On the other hand, in the case of projective (or weak) measurements, there is always a non-zero constant entanglement that survives in the long time limit even in the area-law phase (except for the extreme case with p m = 1).
We remark that the projective (or weak) measurement models can be used as a sampling-based quantum-trajectory method to simulate 1D noisy RCS subject to CPTP noise maps. For instance, a single-qubit dephasing channel can be understood as a channel that results from performing a non-destructive measurement in the computational basis (or Z basis) with a measurement probability p m = 2p D and then forgetting about the measurement outcome, i.e., Note that (Î +Ẑ)/2 and (Î −Ẑ)/2 respectively correspond to the projection operators |0 0| and |1 1| that project the system onto a space associated with each measurement basis state. We remark that this method can be efficient only when the dephasing error probability p D is above a certain threshold value p D > p (th) m /2 0.08. However, we stress that the existence of a non-zero threshold value is specific to this quantum-trajectory method and is not intrinsic to the 1D noisy RCS model itself. Our results indicate that while each quantum trajectory in the measurement model yields a pure state with volume-law entanglement below the threshold (p D < p (th) m /2 0.08), the mixed state that results from putting together all the pure states in all the trajectories will have an area-law MPO entanglement entropy (see also the beginning of Section VIII in Ref. [72] for a related discussion). Hence, our method is effective as it directly simulates the mixed states via MPOs and thus maximally takes advantage of the reduction in quantum correlations due to CPTP noise maps. In contrast, such reduction in quantum correlations at the mixed state level is not exploited in the quantum-trajectory method based on the measurement models.
In another related recent work [45], an MPS method was used to simulate 1D and 2D random quantum circuits. Our work fundamentally differs from this work because in the latter, gate errors are introduced because of the truncation of small singular values when updating MPSs after each twoqubit gate, not because of the CPTP noise channels such as depolarization channels. In other words, Ref. [45] aims to approximately simulate an ideal random quantum circuit and is not concerned with what types of errors are introduced in the classical simulation as long as a non-zero fidelity is attained (see also the discussion on the difference between the total variation distance noise and the CPTP noise channels in Section II). For instance, this work has demonstrated that depth-20 2D random quantum circuits with 54 qubits can be efficiently simulated up to a global fidelity F ≥ 0.002. However, as was pointed out in Ref. [45], such a remarkable performance is specific to the setting that they considered where each two-qubit gate is fixed to be the CZ gate and only single qubit gates are chosen to be Haar-random. That is, at the quantitative level, the excellent performance is attributed to a simple grouping strategy available for CZ gates and it will be more costly to simulate a system of the same size if the two-qubit gates are chosen to be Haar-random. Moreover, the method in Ref. [45] is shown to be efficient only above a certain non-zero threshold value of the gate error rate, i.e., for > ∞ ∼ 0.01, where is the error rate per gate.
In contrast, while we restricted ourselves to a 1D setting, we specifically addressed gate errors that are introduced due to a practically relevant CPTP noise map (hence we use MPOs instead of MPSs), not due to the bond-dimension truncation in MPO simulations. Moreover, we considered general Haarrandom two-qubit gates as opposed to CZ gates. Also in regard to the bond dimension truncation, we demonstrate that a constant bond dimension suffices and the errors associated with the bond-dimension truncation are insignificant in our case. This is again thanks to the fact that we directly simulate mixed states with MPOs and maximally take advantage of the correlation reduction caused by CPTP noise maps. Moreover, our numerical results and heuristic analysis suggest that there does not exist a non-zero threshold value of the gate error rate.
Lastly in the context of IQP, Ref. [73] showed that most IQP circuits can be classically simulated approximately if they are subject to depolarization errors with a non-zero error rate. In particular, the simulation efficiency comes from the fact that the output probability distribution P ( x) of an IQP circuit becomes sparse in the Fourier-transformed basis, since highorder Fourier coefficients are suppressed exponentially due to the depolarization errors. A similar technique was used in the context of RCS in Refs. [74,75] to claim efficient classical simulability of noisy random quantum circuits (see also Ref. [76] and the supplementary material of Ref. [40] for a related discussion).
We remark that the methods in Refs. [74,75] are applicable only in the large circuit depth limit. In contrast, while our MPO simulation is only applicable to 1D settings, it works for any circuit depth D. Most importantly, although our MPO simulation is applicable to deep circuits, we are not primarily concerned with the large circuit depth limit. Instead, as illustrated in Fig. 1, we have focused on identifying the optimal circuit depth where the maximum non-trivial quantum correlation is attained and understanding how hard (or easy) it is to simulate such an optimal regime.

VI. SUMMARY AND OPEN QUESTIONS
In this work, we have numerically investigated the computational power of 1D noisy quantum systems by using MPOs. The key observation is that the maximum achievable MPO entanglement entropy S max is bounded by a constant S max,∞ (p) that is independent of the system size, but depends only on gate error rates (see Fig. 8). Our numerical results thus imply that the classical simulation cost of a 1D noisy system increases exponentially in the system size n only until the system size n reaches a certain length scale. Above the length scale, which is determined solely by the gate error rate p, adding more qubits brings about only polynomial increase in the classical simulation cost. We have also provided a heuristic argument which suggests that the maximum achievable MPO entanglement entropy would scale as S max,∞ (p) ∝ p −1/α in the small p limit for some α > 0. This scaling relation implies that the classical simulation cost of a 1D noisy system would increases exponentially as the gate error rate p decreases. Thus, decreasing the gate error rate p can make the classical MPO simulation practically impossible, given that the system size n reached a certain length scale. For this reason, we were not able to numerically investigate the low gate error rate regime with p ≤ 0.06.
An immediate future research direction is thus to numerically explore the low gate error regime with an advanced computing resource and see if the same behavior in Fig. 8 and the scaling in Eq. (45) hold. Moreover, it would be ideal to make the heuristic analysis in Subsection IV B more rigorous. These studies will allow us to understand the scaling of the length scale where the saturation of MPO entanglement entropy occurs as a function of the gate error rate p for the currently available gate error rates p ∼ 10 −3 − 10 −2 and the cost of classically simulating such 1D noisy quantum systems via MPOs.
Another important open question is whether the same conclusions hold also in the case of 2D noisy RCS, which is more relevant to the currently deployed state-of-the-art superconducting qubit systems. A natural way to extend our results to the 2D cases would be to use the projected entangled pair operators (PEPOs), which are a mixed state generalization of the projected entangled pair states (PEPSs) [77,78]. Note that generalizing the heuristic arguments given in Subsection IV B to the 2D cases, we find that the PEPO entanglement entropy (defined similarly as in the case of the MPO entanglement entropy) of a region A will scale as where |∂A| is the length of the boundary of the region A, D is the circuit depth, p is the gate error rate, and c is a constant. The prefactor D 2 is due to the fact that O(D 2 ) qubits are contained in a causal cone formed by the circuit depth D. Also, if all the O(D 3 ) gates in the causal cone contribute equally, the exponent α would be given by α = 3. Eq. (48) then implies that the PEPO entanglement entropy S A (D, p) is maximized at an optimal circuit depth and thus the maximum achievable PEPO entanglement entropy scales as Numerically investigating whether these scaling relations hold will be very helpful for understanding the computation power of the state-of-the-art superconducting qubit systems in a planar architecture.
If the maximum achievable PEPO entanglement S A (p) follows an area law as suggested in Eq. (50), the corresponding mixed states may be efficiently described by a PEPO with a constant bond dimension χ that depends only on the gate error rate p. However, unlike in the case of MPO, exactly computing an observable from a PEPO is not feasible because exactly contracting PEPSs is #P complete in the worst case [79] and in the average case [80]. Nevertheless, these hardness results do not immediately rule out the possibility of efficient and approximate simulation of 2D noisy RCS because it may be possible to efficiently contract the output PEPOs approximately.
We remark that the computational complexity of 2D RCS with a constant circuit depth has recently been studied in Ref. [81]. In particular, Ref. [81] suggests that all but superpolynomially small fraction of constant-depth 2D random quantum circuits can be simulated approximately and provides numerical evidence supporting the claim. While the two-qubit gates are assumed to be noiseless in this work, the results of Ref.
[81] on constant-depth 2D circuits are very relevant to understanding the complexity of 2D noisy RCS. This is because in the presence noise, maximum non-trivial quantum correlations may be achieved at a constant circuit depth. Moreover, since the simulation algorithms for 2D systems given in Ref. [81] are based on simulating 1D systems using MPSs, it would be interesting to see if one can integrate the methods given in Ref. [81] with the MPO approach presented in our work.
We finally remark that the efficient approximate simulability of typical 2D noisy random circuits (which is subject to future studies) is not in contradiction with the feasibility of faulttolerant quantum computing. This is because circuits that are used for fault-tolerant quantum computing are far from being random since they are carefully designed so that gate errors are propagated in a restricted manner. If it turns out that typical instances of noisy 2D RCS with Haar-random two-qubit gates can be efficiently simulated (which is again to be explored), another important related question is whether there exists a special class of robust random circuits that are more feasible than fault-tolerant quantum computing, but are structured enough so that gate errors are propagated in a restricted way and thus are harder to simulate classically.
Define M I l I l+1 ,J l J l+1 as where I l = di l + i l , · · · , J l+1 = dj l+1 + j l+1 . Then, upon the action of the two-qudit CPTP map M, we have as an output state  α l as follows: Note that the index β is replaced by α l . In the case of unitary two-qudit gates, all the other singular values are unchanged and we update Γ [l]I l α l−1 α l and Γ [l+1]I l+1 α l α l+1 accordingly, leaving all the other Γ parameters unchanged [44]. However, for general two-qudit CPTP maps, this is not the case any more and we should update the singular values and the Γ parameters globally.
To further update the MPO parameters (on the left hand side), note that where B [l−1,l] α l−1 ,I l α l is defined as Similarly as above, applying SVD to B [l−1,l] , we get where V [l−1] is a χ×χ unitary matrix and W [l] is a d 2 χ×d 2 χ unitary matrix. Note that the summation index β goes from 0 to χ − 1 since B [l−1,l] is a χ × d 2 χ matrix. Thus, we find the following Schmidt decomposition of the output state with respect to the cut [1 · · · (l − 1)] : [l · · · n]: where |Φ α l−2 as follows: Repeating the same procedure, Γ α l−3 , · · · , λ [1] α1 can be updated in the same way. In the very last step, we should update Γ [1]I1 α1 as follows: (A19) The remaining MPO parameters on the right hand side are also updated in the same way starting from I l+1 α l ,α l+1 is defined as α l+2 as follows: α1 as follows: I1α1 , λ [1] α1 ← λ [1] α1 .
Right edge case: Lastly, let us consider another boundary case with l = n − 1. Note that an input MPO in the canonical form is given by  αn−1 , · · · , λ [1] α1 can be updated in the same way as in the bulk case, following the procedure described in Eqs. (A9)-(A19).