Error-mitigated deep-circuit quantum simulation of open systems: steady state and relaxation rate problems

Deep-circuit quantum computation, like Shor's algorithm, is undermined by error accumulation, and near-future quantum techniques are far from adequate for full-fledged quantum error correction. Instead of resorting to shallow-circuit quantum algorithms, recent theoretical research suggests that digital quantum simulation (DQS) of closed quantum systems are robust against the accumulation of Trotter errors, as long as local observables are concerned. In this paper, we investigate digital quantum simulation of open quantum systems. First, we prove that the deviation in the steady state obtained from digital quantum simulation depends only on the error in a single Trotter step, which indicates that error accumulation may not be disastrous. By numerical simulation of the quantum circuits for the DQS of the dissipative XYZ model, we then show that the correct results can be recovered by quantum error mitigation as long as the error rate in the DQS is below a sharp threshold. We explain this threshold behavior by the existence of a dissipation-driven quantum phase transition. Finally, we propose a new error-mitigation technique based on the scaling behavior in the vicinity of the critical point of a quantum phase transition. Our results expand the territory of near-future available quantum algorithms and stimulate further theoretical and experimental efforts in practical quantum applications.


I. INTRODUCTION
Quantum computation [1] can solve physical problems, e.g., analyzing equilibrium and dynamical properties, of strongly-correlated quantum many-body systems more efficiently than its classical counterpart. One of the prominent computational frameworks is the so-called digital quantum simulation [2], where the unitary time evolution of local Hamiltonians is first discretized into small time steps, termed as Trotter steps, and then decomposed into a number of native gates, based on the Trotter-Suzuki formulas [3][4][5]. The quantum computational complexity, in terms of the gate count, is estimated to be polynomial in the number of constituents particles, while the classical simulation in general consumes an exponential amount of resources. With programmable quantum devices, this algorithm can be applied to a wide range of research fields, from condensed matter physics [6] to quantum chemistry [7,8]. As the quantum technologies improve, the number of qubits involved in digital quantum simulation has just increased to more than 16 in recent experiments [9,10].
Although powerful and flexible, digital quantum simulation always uses deep circuits, with the number of quantum gates increases proportionally to the number of qubits, the evolution time, and the desired simulation accuracy. It was thus largely maintained that this algorithm would eventually be corrupted by accumulated errors, coming from both the Trotter-Suzuki expansion and imperfect quantum devices. A recent theoretical paper [11], however, suggests that digital quantum simulation is robust against Trotterization errors, in a sense * zhangjn@baqis.ac.cn † yli@gscaep.ac.cn that the deviations of local observables are under control for relatively large Trotter steps. Moreover, it has been shown [11,12] that there is a sharp threshold behavior for the Trotter errors, which can be explained by the transition from quantum localization to many-body quantum chaos. These theoretical findings have strengthened our confidence in the prospect that digital quantum simulation will beat errors and show advantages over classical simulation with noisy intermediate-scale quantum devices [13]. Alternatively, we find that there are physical problems of open quantum systems [14] that are robust against noises. In other words, although using deep quantum circuits, the deviations in the final results of these deepcircuit quantum simulations still remain under control, and the corresponding ideal results can be extracted with certain error-mitigation techniques [15][16][17]. The classical simulation of open quantum systems is harder than that of closed quantum systems, since the storage and manipulation of density matrices generally require more computational resources [18], compared with those of state vectors. As a result, the universal quantum simulation of open quantum systems [19] covers diverse quantum algorithms that are appropriate for demonstrating the potential powerfulness of quantum computers. Here in this paper, we first prove, via perturbation theory, that the deviation in the steady state obtained via digital quantum simulation depends only on the error in a single Trotter step. Using the dissipative XYZ model as an example, we then show that there is a sharp threshold for the error rate below which the ideal steady-state properties and relaxation rates can be recovered by a specific error-mitigation technique [15], i.e., the so-called zero-noise extrapolation [16]. We explain this threshold behavior by the existence of a dissipation-driven quantum phase transition, which is supported by the mean-field theory. Finally, we propose a new error-mitigation technique based on the scaling behavior [20,21] in the vicinity of the critical point of a quantum phase transition, and provide numerical evidence with the mean-field calculation of the dissipative XYZ model.

II. OPEN QUANTUM SYSTEM
Here we briefly recall the basic description of the dynamics of open quantum systems. Consider a quantum system weakly coupled to a Markovian environment. The quantum dynamics is described by the quantum master equation of Lindblad form, whereρ is the reduced density matrix of the system, and L 0 is the Lindbladian. The general form of the Lindbladian is whereĤ ≡ L l=1ĥ l is the system Hamiltonian and kÂ k is the dissipator of the system, operatorÂ k with non-negative coefficients γ i quantifying the dissipation strength. Given an initial stateρ(t = 0) ≡ρ in , the reduced density matrix at time t can be formally obtained aŝ

III. STEADY-STATE PROBLEM
In the theory of open quantum systems, the steadystate problem, i.e., obtaining properties of steady states of given physical models, is of particular relevance. The quantum master equation in Eq. (1) can be formally solved to obtain the reduced density matrix of the system at time t asρ(t) = e L0tρ in , withρ in being the initial density matrix. Intuitively,ρ(t) eventually converges to the steady stateρ 0 , which is by definition annihilated by the Lindbladian, i.e., L 0ρ0 = 0. In this paper, we consider only open system with single steady state, therefore the steady stateρ 0 = lim t→∞ e L0tρ in do not depend on the initial density matrix ρ in . To investigate the speed of the convergence, it is convenient to use the Hilbert-Schmidt space, where an arbitrary operatorÂ, which is a d×d matrix in the Hilbert space, is mapped to d 2 -entry column vectors |A , with d being the dimension of the Hilbert space and the inner product A|B defined as Tr A † B . In this representation, the Lindbladian L 0 in Eq. (2) becomes a d 2 ×d 2 matrix L 0 , and the steady state |ρ 0 is the right eigenvector of L 0 with the eigenvalue λ 0 = 0. Formally solving the Eq. (1) in the Hilbert-Schmidt space, we can prove that with an arbitrary norm. The lowest relaxation rate Γ = min {−Re (λ α ) |α = 0} determines the speed of the convergence, where λ α 's are the eigenvalues of L 0 . (see Appendix B). Obtaining the steady state by a classical computer needs to manipulate L 0 and storing |ρ . The classical computation costs exponential computational resources both in space and time with respect to the number of qubits N in the system, since L 0 and |ρ are respectively a 2 2N × 2 2N complex matrix and a 2 2N complex vector. Here we propose that the exponential cost can be alleviated by a near-term noisy quantum processor.

IV. QUANTUM ALGORITHMS
A universal quantum computer can prepare the steady state by directly evolving the quantum master equation for a long enough time T , such that the norm ρ (T ) −ρ 0 < δ, with δ > 0 is the predetermined accuracy goal. In contrast to the evolution of closed quantum systems, the basic operation block e L0τ , with τ being a small time interval to be determined by the error tolerance, cannot be represented as a unitary operator. Theoretically, a universal quantum computer can realize arbitrary quantum operations on N qubits with either 2N ancilla qubits [22] or a single ancilla qubit with 2N feed-back control cycles [23]. Together with Eq. (4), it is evident that to prepare the steady state, a universal quantum computer would require O(N ) qubits and O log δ −1 runtime.
One of the most promising ways that universal quantum computers simulate the time evolution of quantum systems is using the Trotter-Suzuki decomposition. Here we generalize the framework proposed in Ref. [2] to simulate the dynamics of open quantum systems. Specifically, we realize the following quantum operation, also called a Trotter step, on the universal quantum processer, with . . , M with M = K + L) are defined as the ideal native interactions being driven to implement e D k τ and U l (τ ).
For near-term noisy quantum processors, the real operations can be represented as the ideal operation slightly corrupted some noise process. Mathematically, we can express the real operation as e Gmτ with G m = G id m +r m E m , where r m > 0 is a small positive value that quantifies the strength of the noise for the m-th quantum operation, whose generator is formally denoted as E m (see Appendix F for the treatment of the depolarizing error channel). Using the Magnus expansion, the real effective Lindbladian on the near-term noisy quantum processor is where by definition L 0 = M m=1 G id m . Note that here we treat both r and τ as small quantities of the same scale. Intuitively, the deviation of L eff from L 0 is determined by both the noise strength r and the length of the Trotter step τ . Note that the high-order terms can be found in Appendix A. The effect of Trotter step τ has been widely studied and it has been shown that there is a sharp threshold for τ for closed system [11]. There is also such a threshold for τ in open systems, however, we choose a fixed τ = 0.01 below the threshold and focus on the unexploited field on the effect of noise strength r.
When the total evolution time T is sufficiently long, the state of the system converges to the steady state of the real effective Lindbladian, which satisfies L effρeff = 0. Using Eq. (4) and the triangle inequality, we find the difference between the final state on the noisy quantum processor and the ideal steady state satisfies the following inequality, where Γ eff is the lowest relaxation rate of the real effective Lindbladian L eff . The second term in the above equation can be arbitrarily suppressed by increasing the evolution time T , thus the error in the prepared steady state is dominated by the difference betweenρ eff andρ 0 . Next, we will show that the steady stateρ eff and the relaxation rates λ eff,α of the noisy effective Lindbladian L eff can be obtained from the corresponding ideal ones by perturbation theory, thus the deviations of the prepared steady state and the measured relaxation rates from the ideal ones are bounded by the error of a single Trotter step, even though the quantum algorithms consist deep quantum circuits.

V. ZERO-NOISE EXTRAPOLATION
To use the perturbation theory, we first write L eff = L 0 + L and treat L as the perturbation. The steady state ofρ eff can be expressed as a power series of superoperators working on the ideal steady stateρ 0 , where L −1 0 is the superoperator corresponding to the generalized inverse of L 0 in the Hilbert-Schmidt space. Therefore, the expectation of any observableÔ can be written as Similarly, the relaxation rates, i.e., the real parts of the eigenvalues of the Lindbladian, can also be expressed as The power series ofρ eff will converge within a few orders under the condition L −1 0 L 1, with L −1 0 ∼ 1/Γ and L ∼ r. Thus it is anticipated that if the noise strength r is smaller than the lowest relaxation rate Γ, the steady stateρ 0 of the ideal target Lindbladian L 0 can be inferred by extrapolation with linear or low-order polynomial functions, given noisy data obtained with the intrinsic error rate r 0 and boosted error rates cr 0 with c > 1. (See Appendix G for details.) In the following, we will give a concrete example of using the error mitigation technique to obtain the physical properties of the ideal steady state, with a programmable universal quantum processor.

VI. DISSIPATIVE XYZ MODEL
We consider the dissipative version of the anisotropic spin-1 2 Heisenberg model [24], also called the XYZ model, on a two-dimensional square lattice as an example. The ideal Lindbladian L 0 takes the general form in Eq. (2), with the Hamiltonian for the spin lattice beinĝ with the row and column parts defined as, withσ α i,j (α = x, y, z) being Pauli matrices on the qubit at the i-th row and the j-th column. Here the anisotropic Heisenberg interaction strength on the nearest-neighbor sites are quantified by J α (α = x, y, z). Besides the above Hamiltonian, each spin is also subject to a dissipation process governed by the following dissipator, where γ is the intrinsic single-site dissipation rate. The steady stateρ 0 can be obtained from an arbitrary initial state ρ in after the following Trotterized evolution: where quantum gate G m is either generator of local column (row) This model plays an important role in the field of unconventional magnetism and nonequilibrium phase transition [24]. In contrast to the equilibrium case, each spin in the system keeps precessing around the effective magnetic field stemming from its surrounding spins. Thus the phase transition can be understood from the aspect of the competition between the precession and the dissipation, which is quantified by the dimensionless transversal aspect ratio Specifically, in the thermodynamic limit, the critical point is g cri = 0.062 for the dissipative XYZ model with J z = γ [25,26], and the steady state lies in the paramagnetic (ferromagnetic) phase when g < g cri (g > g cri ).

VII. EXPERIMENT PROTOCOL
There are many theoretical efforts [27][28][29][30] been made to simulate Moarkovian open quantum systems with unitary operations. Recently, experimental simulation of a two-level open quantum system is demonstrated in a superconducting circuit [31]. It is not hard to generalize these ideas to the dissipative XYZ model. The steady state of the dissipative XYZ model on a L × L square lattice can be prepared with a programmable quantum processor consisting of L × 2L qubits. The quantum circuit for the preparation of the steady state on a 2 × 2 spin lattice is shown in Fig. 1 as an example, where the construction only requires local interactions and can be straightforwardly generalized to L > 2. In the Trotterized evolution (13), the order of quantum gate G m can be chosen arbitrarily. In numerical calculation, we choose to group gates corresponding to column (row) Hamiltonian together as e G c(r) τ s and gates corresponding to dissipators together as e G d τ s as shown in Fig. 1.

VIII. SMALL-ANGLE GATES AND STANDARDIZED CIRCUITS
Having designed the circuit architecture, we can implement each module in the circuit by small-angle gates or a standard universal gate set. Mainstream physical platforms for quantum information processing, like the superconducting circuit system and the trapped ion system, have their own native gates with continuously tunable parameters. Here we take the controlled-phase (CP) gate, defined as CP(θ) = exp (−iθ|11 11|), in the superconducting circuit system as a representative example. With simple algebra, it is clear that CP(θ) is equivalent to exp − iθ 4σ z 1σ z 2 and single-qubit rotations. With single-qubit π/2-rotations on the transverse axis and the CP gate, we can realize Ising-type interaction alone all three axis. (See Appendix D for details and realization of other gates.) There is a column of ancillas besides each column of qubits. In this layout, the column HamiltonianĤc can be directly implemented by a universal gate set among qubits. (c) Evolution governed by the row Hamiltonian. Since neighboring qubits are interleaved by ancillas in this layout, the terms in the row HamiltonianĤr need to be implemented by using the ancillas as mediators. This can be achieved either by swap-like gates on the ancilla-qubit pairs or using the ancillas as tunable couplers [32]. (d) Evolution governed by on-site dissipators. A recent experiment [33] shows that with a single ancilla, it is possible to implement arbitrary quantum operations on a qubit.

IX. NOISE MODEL
We consider the depolarizing error model for the noisy quantum processor. Specifically, in the numerical simulation of the quantum circuit for the preparation of the steady states of the dissipative XYZ model, every quantum gate is followed by a quantum operation representing the depolarizing error, e Gmτ = e rEmτ e G id m τ , with the following error generator where 1 1 Am is the identity operator on the Hilbert space of the qubit set A m , which are involved in the m-th quantum operation. Note that the number of qubits in the set |A m | = 1 (2) if the m-th quantum operation involves one (two) qubit(s). For simplicity, we assume a homogeneous strength for all quantum operations, i.e., r m = r. We believe that introducing different noise strengths for each quantum operation would not cause quantitative changes in our conclusion.
Before presenting the numerical results, we wouldd like to make some more comments on the choice of choosing depolarizing as the representative noise model. Although being idealized, the depolarizing error model is widely used in the quantitative analysis of the performance of quantum algorithms. The reason for this choice is twofold. First, the depolarizing error model introduces all possible Pauli errors, thus provides an unbiased criterion for the algorithmic robustness. Second, by Pauli twirling [34], it is possible to convert generic noise into the Pauli channel, of which the depolarizing model is a special case. Finally, we mention that though we choose depolarizing error as the representative, our conclusions which will be shown below hold for various noise models. In Appendices F and G, we present other noise models and corresponding numerical results, from which we see that the details of error mitigation for different noise models depend on their symmetries.

X. NUMERICAL RESULTS
We numerically simulate the process of preparing the steady state and obtaining the values of observables of the dissipative XYZ model by a noisy quantum processor. In Fig. (2), we show the performance of the zero-noise extrapolation technique with the measured steady-state magnetization, on a 3 × 3 qubit lattice. Due to the finite-size effect, the magnetization curve smoothly increases, instead of abruptly jumping, from −1 towards 0, as the aspect ratio g changes. Specifically, we show in the Fig. 2(a) the steady-state magnetization for the exact and the noisy prepared steady statesρ 0 andρ eff , and the steady-state magnetization obtained by zero-noise extrapolation at r 0 = 0.01. To make a closer investigation, we show the discrepancies of the noisy and the extrapolated data from the true values in Fig. 2(b). We find that the discrepancy between M 0 and M eff can be greatly decreased by the zero-noise extrapolation results M ex . The zero-noise extrapolation technique works quite well in the regimes far from the critical point, where the discrepancies are substantially suppressed. While in the vicinity of the critical point, the discrepancies become larger after extrapolation. To find the reason for the failure of the extrapolation, we choose three parameters, specified by stars in Fig. 2(b), and plot the magnetization as a function of the modular error rate r in Figs. 2(c)-(e), respectively. For parameters that lie deep in the paramagnetic and ferromagnetic regime, the steady-state magnetization exhibits monotonic behavior as the modular error rate r increases. However, in the vicinity of the phase transition, the curve first goes down and then bends upward, which results in the failure of the extrapolation when the raw modular error rate is close to or larger than the inflection point. The inflection point occurs at r = 0.03 in Fig. 2(d), therefore in order for the zero-noise extrapolation technique to work, we need to choose an intrinsic error rate r 0 such that 2r 0 < 0.03. The threshold below which the correct result of observables can be recovered by the quantum error mitigation is determined the minimum r at which the inflection point occurs as g varies. We find that an acceptable threshold is r 0 = 0.01. Then we use the mean-field method [24] to reveal the physical insights behind the distinct behavior in different parameter regimes. Based on the effective Lindblad quantum master equation, the equations of motion for the expectation values of the local spin components can be written as follows, where σ x i and σ y i serve as order parameters characterizing the transition between the paramagnetic and ferromagnetic phases. We then introduce the lowest-order of the mean-field approximation to express the correlation function σ α the modular error rate r, and present the numerical results in Fig. 3. Starting from the vicinity of the critical regime, with g = 0.1 as marked by the red asterisk in Fig. 2(b), it is clear that there is a noise-driven phase transition as r increases. Moreover, the critical point r cri predicted by the mean-field calculation is quantitatively consistent with the inflection point in Fig. 2(d).
Thus we believe the noise-driven phase transition well accounts for the failure of the zero-noise extrapolation. Also inspired by this finding, we propose to use the scaling behavior to mitigate errors. The order parameter of the dissipative XYZ model is In the vicinity of the critical point g cri , the order parameter m(g) can be described by a power-law function, m (g) ∝ |g − g cri | β , with β being the critical exponent. Consider the case that the noises introduced by the realistic quantum processor possess the same symmetry as the model being simulated, for example, the dissipative XYZ model and the depolarizing error channel which both have the Z 2 symmetry. (See Appendix F for other noise examples.) Then the error-perturbed order parameter will still retain the power-law dependency on not only the system parameter g but also the error rate r, i.e., where the r-dependency comes from the perturbation to the critical point and the critical exponent. To mitigate the effect of the errors, we can repeat the measurement of the order parameter under several different values of the noise strength r, e.g., the intrinsic noise strength r = r 0 and boosted ones r = r ≡ cr 0 with c > 1, and fit the data with power-law functions to obtain g cri (r) and β(r), based on which we can infer the true values g cri and β by extrapolation to the zero-noise limit (r = 0). After establishing the error-mitigation technique based on the scaling behavior, we continue to investigate its performance in predicting true values of physical quantities in the whole parameter regime. We first accumulate data for the order parameters, i.e., the expectation values of the local spin components σ x and σ z , for the intrinsic noise strength r 0 and a boosted noise strength r = 2r 0 using the above mentioned mean-field method. By utilizing both error-mitigation techniques, i.e., the direct extrapolation and the extrapolation based on the scaling behavior, we obtain the zero-noise results as shown in Fig. 4, and find the scaling-behavior-based extrapolation correctly recovers the true position of the critical point, while the direct extrapolation performs well on predicting the exact values of the order parameter when the parameter g goes deeply into the ferromagnetic phase. (It also predicts values correctly in paramagnetic phase where the order parameter vanishes.) For comparison, we also put the true values of the order parameters. We can clearly see that the effect of the symmetry-preserving noise on the order parameter is that it moves the critical point towards the ferromagnetic phase. The failure of the scaling-behavior-based extrapolation outside the critical regime is straightforwardly connected with the invalidation of the scaling behavior outside the critical regime. The scaling-behavior-based extrapolation is helpful as it can be used to estimate the critical point which is of vital importance in many physical problems. We show in Fig. 5 the critical points g cri of the dissipative XYZ model under the influence of various noises and evaluated via different extrapolations. (See Appendix G for details of noises and extrapolations.) The estimation of the critical point predicted by the linear extrapolation is quite accurate when r 0 is below the threshold 0.01. Note the deviation of the estimated critical point can be made smaller above the threshold if we use quadratic extrapolation. It seems that the scaling-behavior-based extrapolation performs well on estimating the critical point even though the value of the order parameter itself is not so well predicted when the intrinsic error rate r 0 exceeds the threshold.

XI. SPECTROSCOPY
Finally, we show by numerical simulation that we can obtain spectroscopic information about the ideal Lindbladian L 0 with a noisy quantum processor. Starting from a random initial density operator, we first calculate the nonequilibrium dynamics of expectation values of local physical observables, e.g., the local spin component σ z , evolving towards the steady stateρ eff (r) on a noisy quantum process with the intrinsic and boosted noise rates r 0 and cr 0 . These noisy nonequilibrium dynamics can be fitted with a multimodal decaying exponential function, i.e., where A α ∈ R, θ α ∈ [0, 2π] and λ α ∈ C are the parameters of the model, and we use the positive integer p to control the complexity of the model. Note that Re (λ α ) ≤ 0 for a physical Lindbladian. We use the matrix pencil method [35] to extract λ α (r) from the nonequilibrium dynamics governed by L eff (r) with r = r 0 and cr 0 . The error-mitigated estimation of the true eigenvalues of the target Lindbladian L 0 can be obtained by linear extrapolation to the zero-noise point with r = 0. By choosing different initial density operators and physical quantities in the nonequilibrium dynamics, we can extract most of eigenvalues of the target Lindbladian with small-magnitude real parts, while eigenvalues with largemagnitude real parts are difficult to extract due to their fast decay. As shown in Fig. 5(b), the extrapolated eigenvalues is well consistent with those of the true Lindbladian, i.e. lim r→0 λ α,eff (r) λ α .

XII. CONCLUSIONS
Starting with an initial state, we can simulate its time evolution via DQS in principle and extract any relevant physical quantity we want. The availability of only quantum computers constituting of noisy intermediate-scale quantum devices seems to prevent us from achieving this ambition. Due to the accumulation of errors, the simulation of long-time evolution is unreliable. The question is if we can find cases where physical quantities are not suffered from the breakdown of long-time evolution. In this paper, we show that the steady-state problem of open systems meet our demand. It is robut against error accumulation and we can extracted accurate values of observables via error mitigation.
Our paper consists of two main parts. In the first part, we propose a perturbation theory, by which we prove that the deviation of the steady state depends only on the error in a single Troter step, regardless of the evolution time of DQS. This result has a simple physical interpretation. Error accumulation may corrupt the simulation of long-time evolution, but the steady state, as the goal of evolution, is not affected by it. The deviation of the noisy steady state from the exact one depends only on the difference between the noisy Lindbladian L eff and the exact one L 0 , which has nothing to do with the evolution time. Further, we can express the noisy steady state as a power series of the single-step error. Therefore, the exact steady state and any relevant observables can be obtained by error mitigation. In the second part, we take the dissipative XYZ model as an example to convince us of the theoretical results. When the strength of error is below a threshold, the perturbation theory works. By error mitigation and a technique based on scaling behavior, we can recover the exact values of magnetization, relaxation rate, and the location of the phase transition point.

XIII. ACKNOWLEDGEMENT
We acknowledge the support of the National Natural Science Foundation of China (Grant No. 11875050, 12088101, and 92065205) and NSAF (Grant No. U1930403).

Appendix B: Jordan normal form of Lindbladian and time evolution superoperator
For convenience, we express the Lindbladian in the Hilbert-Schmidt space [36][37][38][39]. We use |A to denote the d 2 -dimensional column vector corresponding to each d × d matrix A, such that for two square matrices Here, 1 1 d is the d-dimensional identity matrix. Using the Jordan normal form, we can write L 0 in the block diagonal form where each α denotes a Jordan block, λ α is the corresponding eigenvalue, d α is the dimension of the block, and α d α = d 2 . Here, {|R α,m } and { L α,m |} are dimensionless vectors satisfying L α,m |R α ,m = δ α,α δ m,m and α,m |R α,m L α,m | = 1 1 d 2 .
The evolution superoperator of L 0 is V 0 (t) = e L0t . Using the Jordan normal form, the corresponding matrix representation is Therefore, the state evolves to the steady state as e L0t (ρ in ) −ρ 0 = O e −Γt poly(t) , where Γ = min {−Re (λ α ) |α = 0} determines the speed of the convergence, where λ α 's are the eigenvalues of L 0 .

Appendix C: Perturbation theory
In this section, we derive the perturbation theory of Lindbladian based on the modified Jordan normal form.

Steady state perturbation theory
We express the steady state asρ eff = ∞ k=0ρ (k) , wherê ρ (k) is the k-th-order correction to the steady state. Substituting this expression ofρ eff into L effρeff = 0 gives (L 0 + L )( ∞ k=0ρ (k) ) = 0. Retain terms of the k-th order, we have L 0 (ρ (k) ) + L (ρ (k−1) ) = 0. Here,ρ (−1) = 0 andρ (0) =ρ 0 . Now, we switch to the Hilbert-Schmidt space, and the equation of each order becomes where L 0 is not invertible because of the zero eigenvalue λ 1 . We need thatρ eff is normalized for each order of the correction, i.e. Tr(ρ (k) ) = 1 Then, Therefore, We remark that this generalized inverse makes sure that ρ (k) is traceless, i.e. the steady state with a truncation at any order is normalized.

Non-degenerate perturbation theory
We consider a non-degenerate eigenvalue λ α , i.e., d α = 1. Let λ (k) α and |R (k) α,1 be the k-th order corrections to the eigenvalue and right eigenvector, respectively. Then, where λ α,1 . Retain terms of the first order, we have Multiply this equation by L α,1 | from the left, we obtain the first-order correction to the eigenvalue Multiple by L β,m | with β = α, we have where L β,d β +1 | = 0. Then, Similar to the perturbation theory of Hamiltonian, we have L α,1 |R (1) α,1 = 0, such that the eigenvector after correction is still normalized, i.e., L eff,α,1 |R eff,α,1 = 1, where |R eff,α,1 and |L eff,α,1 are right and left eigenvectors with the corrections. Then, we obtain the first order correction to the right eigenvector We can obtain the correction to the left eigenvector in a similar way. Retain terms of the second order, we have Multiple this equation by L α,1 | from the left, we obtain the second order correction to the eigenvalue

Appendix D: Circuits
In this appendix we show the circuits for implementing the quantum gates for the dissipative XYZ model in Fig. 7. Here and below, we omit the details of the qubit-ancilla permutation in L × 2L lattice which has been discussed in the main part and consider only the realization of two-qubit gate acting on qubit k and l and single-qubit acting on qubit k. We consider two setups. In the first setup, operations are realized using smallangle gates R αα (θ) = exp −i θ 2 σ α ⊗ σ α . We remark that θ α = 2J α τ and cos φ = e −γτ /2 , and θ α , φ 1 in the Trotterization algorithm. In the second setup, operations are realized using controlled-NOT gate and single-qubit gates. The single-qubit rotation gate reads R z (θ) = exp −i θ 2 σ z .

Appendix E: Error twirling and boosting
The Pauli twirling of controlled-NOT gate uses four Pauli gates σ t,q = σ I , σ x , σ y , σ z , where q = a, b, c, d, as shown in Fig. 8(a). We randomly choose Pauli gates σ t,a and σ t,c according to the uniform distribution, and then  we take σ t,c ⊗σ t,d = Λ x σ t,a ⊗σ t,b Λ x , where Λ x is the unitary operator of controlled-NOT gate. By applying the Pauli twirling, we can convert a general error model into the Pauli error model in the form N [Λ x ], where [Λ x ] denotes the completely-positive map of the ideal controlled-NOT gate, and N = α,β=I,x,y,z p α,β σ α ⊗ σ β denotes Pauli errors. Here, p I,I is the fidelity of the gate, and p α,β [(α, β) = (I, I)] is the probability of the Pauli error σ α ⊗ σ β . To implement the error extrapolation, we need to measure error probabilities p α,β using methods such as quantum process tomography.
Suppose we want to effectively realize the error model N = α,β=I,x,y,z p α,β σ α ⊗ σ β , we can compute N e = N −1 N . Note that the map N is always invertible for high-fidelity gate. N e is always a Pauli-error map in the form N e = α,β=I,x,y,z q α,β σ α ⊗ σ β . We remark that q α,β is not always positive. To effectively realize the depolarizing map, we take p I,I = 1 − and p α,β = /15 for (α, β) = (I, I). We can always find a finite such that all q α,β are positive. Intuitively, we can take = 15 max{p α,β | (α, β) = (I, I)}. Similarly, by taking N as the error model with increased error (e.g. corresponding to 2r raw ), we can boost the error.
Pauli twirling circuit of the operation e −i[H k,l ,•]τ is shown in Fig. 8(b). This operation can be realised us-ing either small-angle gates or controlled-NOT gate, and the twirling circuit works for both cases. To implement the twirling, we apply Pauli gates σ s = σ I , σ x , σ y , σ z on both qubits before and after the operation, according to the uniform distribution. It is similar for the operation e D l τ . Pauli twirling circuit of the operation e D l τ is shown in Fig. 8(c). This operation can be realized using either small-angle gates or controlled-NOT gate, and the twirling circuit works for both cases. To implement the twirling, we apply Pauli gates σ s = σ I , σ z before and after the operation, according to the uniform distribution. Suppose For depolarizing error on controlled-NOT gate, because the map of two-qubit depolarizing error commute with single-qubit gates and two-qubit gates on the same two qubits, the error model of whole circuits in Figs. 7(b) and (d) is depolarizing. In this case, further twirling operations using circuits in Figs. 8(b) and (c) are unnecessary.

Appendix F: Noise models
In this section, we show how we deal with the errors rE in the quantum operation. Besides a description of depolarizing error which is an equal superposition of Pauli errors, we also include a general case of random Pauli error and transverse damping error which has no Z 2 symmetry. At last, we make some comments on the non-Markovian noise.

Depolarizing error
As described in the main part, the real operation used in the universal quantum processer is e Gmτ , with G m = G id m + r m E m and the generator of depolarizing error is where we use a superscript below "d" to denote the depolarizing error, and "r" for random Pauli error, "t" for transverse damping error in the next subsection. In each Trotter step of quantum algorithm, we apply the real operation on the system. We choose r m = r for both single-qubit gates and two-qubits gates without loss of generality.
According to the perturbation theory, the steady state can be expressed aŝ Therefore, M (r) is also a power series of the noise strength r, and the linear ansatz is the simplest approximation. In Fig. 9, we show the zero extrapolation and scaling extrapolation with the quadratic ansatz M (r) = M 0 + br + ar 2 .
In Fig. 10, we show the zero-noise extrapolation and scaling extrapolation of random Pauli error and transverse damping error. It can be seen that the error without Z 2 symmetry will destroy the phase transition. Therefore we modify Eq. (19) as m (g, r) ∝ |g − g cri (r)| β(r) + a(r) + b(r)r, (G5) where a(r) accounts for the deviation of "order parameter" from zero. There is a slight positive slope b(r) which may not be directly seen from the plot. In Fig. 11, we show some lowest eigenvalues obtained by the mean-field method. Some of the eigenvalues differ much from that obtained on the 3 × 3 qubit lattice.