Limitations of probabilistic error cancellation for open dynamics beyond sampling overhead

Quantum simulation of dynamics is an important goal in the NISQ era, within which quantum error mitigation may be a viable path towards modifying or eliminating the effects of noise. Most studies on quantum error mitigation have been focused on the resource cost due to its exponential scaling in the circuit depth. Methods such as probabilistic error cancellation rely on discretizing the evolution into finite time steps and applying the mitigation layer after each time step, modifying only the noise part without any Hamiltonian-dependence. This may lead to Trotter-like errors in the simulation results even if the error mitigation is implemented ideally, which means that the number of samples is taken as infinite. Here we analyze the aforementioned errors which have been largely neglected before. We show that, they are determined by the commutating relations between the superoperators of the unitary part, the device noise part and the noise part of the open dynamics to be simulated. We include both digital quantum simulation and analog quantum simulation setups, and consider defining the ideal error mitigation map both by exactly inverting the noise channel and by approximating it to the first order in the time step. We take single-qubit toy models to numerically demonstrate our findings. Our results illustrate fundamental limitations of applying probabilistic error cancellation in a stepwise manner to continuous dynamics, thus motivating the investigations of truly time-continuous error cancellation methods.


I. INTRODUCTION
Dynamics of quantum systems are unavoidably noisy due to the interaction with the environment.Elements of quantum computers, therefore, are also subject to a variety of noises.The main techniques of reducing the effects of noise belong to two categories, quantum error correction [1] and quantum error mitigation [2].Quantum error correction encodes the quantum information in a much larger Hilbert space in order to achieve faulttolerance.Quantum error mitigation, on the contrary, makes use of a large number of copies of the quantum system and some knowledge of the noise structures in order to achieve a noise-free ensemble average result via statistical post-processing.Given the current and near-future error rates for a single quantum register [3], the implementations of various quantum error mitigation methods on actual quantum hardware [2] are considered to be an important milestone on the way to the long-term goal of fault-tolerant quantum computing.
One important application of quantum computing before full fault-tolerance can be achieved, is quantum simulation, especially the simulation of quantum dynamics [4].Indeed, different quantum error mitigation methods have been applied to different implementations of quantum computers to simulate unitary dynamics [5][6][7][8][9], with the possible potential of showing quantum advantages over classical computers.Recently, it has been proposed that quantum error mitigation can be applied to partially mitigate the noise of the quantum computer, so that open dynamics can be simulated, which have applications in many aspects in quantum chemistry and quantum biology [10].
Although widely used, methods of quantum error mitigation are known to suffer from an exponential scaling of the sampling overhead [11][12][13], and to what extent this can be optimized remains an open question.On the other hand, as quantum error mitigation aims at removing noise, how its effects change with different unitary parts of the dynamics is rarely explored, partially due to the exponential sampling cost mentioned above that poses significant numerical challenge.This is especially important for simulating the dynamics, as according to the theory of open dynamics [14,15], in general it is impossible to find a Hamiltonian-independent map that precisely converts a noisy time evolution to a noiseless one.
In this paper, we focus on the simulation errors of open and closed dynamics when probabilistic error cancellation, which is one of the quantum error mitigation methods, is applied to reshape the noise distribution of a noisy quantum simulation process.In particular, we assume that the number of samples is infinitely large, such that the probabilistic error cancellation is implemented in an ideal way.We describe a general framework to show that the simulation errors depend on the commutators between different superoperators governing the dynamics, and the contributing commutators are different depending on whether the simulation setup is digital or analog.We then illustrate our results numerically by taking a single-qubit toy model as an example, where analytical expressions are possible and the large sampling overhead for quantum error mitigation is not a limiting factor.Although the simulation errors we investigate here are unlikely to dominate over the numerical errors introduced by the sampling overhead for realistic arXiv:2308.01446v2[quant-ph] 30 Jan 2024 models, our work points out the intrinsic limitations of combining step-wise protocols to continuous dynamics, which can only be overcome by applying time-continuous protocols instead.

II. THE GENERAL FRAMEWORK
We consider the simulation of a Lindblad master equation that describes weak couplings between the system and a Markovian environment [14,16], where Ĥ is the Hamiltonian corresponding to the unitary dynamics of the system, Lk are Lindblad operators with error rate γ k , For simplicity, we have taken ℏ = 1, and have assumed that Ĥ and Lk are time-independent.We rewrite the master equation into the superoperator form by vectorizing the d×d density matrix ρ into the columnordered where L h is the superoperator describing the unitary part of the dynamics governed by the Hamiltonian, and L d is the superoperator describing the noise part of the open dynamics to be simulated.Both L h and L d are d 2 × d 2 matrices acting on the vectorized density matrix ρ from the left.As we have assumed that Ĥ and Lk are timeindependent, L h and L d are time-independent as well.
The solution of Eq. ( 2) is therefore Note that, according to the Baker-Campbell-Hausdorff formula [18], exp(L h t + L d t) = exp(L d t) exp(L h t) works only if L d and L h commute.Also note that the case of simulating closed dynamics can be included by simply taking L d = 0. Typically, a quantum computer aims at simulating the unitary dynamics L h , but unavoidable sources of noise imply that the quantum computer is in fact simulating some sort of open dynamics.We consider applying quantum error mitigation to attenuate the noise of the quantum computer, such that effectively we can simulate the dynamics governed by L h + L d rather than the intrinsic noise of the quantum computer.Specifically, we choose the method of probabilistic error cancellation [19].It aims at implementing a map that cancels the effect of the noise channel.As this map is not physical, it can only be realised statistically.A large number of trials are required.Within each trial, some unitary operators are sampled and applied to the state.Taking the statistical average over all the trials, together with classical postprocessing, will effectively implement the desired nonphysical map.We will explicitly illustrate these steps in Sec.III A with the example of a single qubit.As probabilistic error cancellation is designed to act in discrete time steps [6,10,20], we consider the time interval ∆t, such that the state of the system at time t + ∆t is given by MCρ(t), where M is the superoperator corresponding to the non-physical map to be realized by probabilistic error cancellation, and C is the superoperator describing the noisy implementation of the unitary dynamics during the time step ∆t.The form of C depends on whether the quantum computer is digital or analog, and will be discussed in the next subsections.The error mitigation superoperator M should be designed such that That is, we want to simulate open dynamics on a noisy quantum computer, which can be executed by mitigating the noisy unitary evolution such that effectively the remaining noise matches the noise coming from the interactions with the environment that are to be simulated.As a special case, this also includes the simulation of closed quantum systems.We are going to show that, if we require M to be independent of L h , which is a standard assumption as error mitigation protocols only aim at cancelling the noise, unless special conditions are satisfied, in general Eq. ( 6) cannot hold precisely.The errors resemble Trotter errors [18], but are in terms of the superoperators corresponding to the unitary part and noisy part of the dynamics.

A. digital quantum simulation
For a digital quantum computer, the noisy implementation of a unitary gate is usually modelled as a unitary layer followed by a layer of the noise channel [6,10].Specifically, where L n is the superoperator describing the device noise channel.
In some cases the device noise channel is defined by parameters that have already taken into account the exponential over ∆t.This will be discussed in detail with the numerical examples in Sec.III B 1. Assuming that exp(L n ∆t) is invertible, the error mitigation superoperator M should satisfy For the simulation of closed dynamics, L d = 0. We therefore require This is independent of the unitary part L h , and only cancels the device noise.On the other hand, for the simulation of open dynamics [10], L d ̸ = 0. Different forms of M are resulted in, depending on the commutator between the superoperators L d and L h .If [L d , L h ] = 0, we have The error mitigation superoperator M is still independent of L h .Note that, satisfying the relation [L d , L h ] = 0 may pose a requirement on the form of L h , but for special L d such as depolarizing noise [6], [L d , L h ] = 0 holds for any Hamiltonian.Finally, if [L d , L h ] ̸ = 0, according to the Zassenhaus formula [17], where we have neglected terms that contain high orders of ∆t in the exponential.The first term on the right hand side of Eq. (11) implies that M has to depend on the Hamiltonian Ĥ. Ignoring this term and the higher order terms will lead to Trotter errors, which depend on the superoperator commutator [L d , L h ] and can be reduced by decreasing the time step ∆t.However, the large sampling overhead of probabilistic error cancellation may put a lower bound on the implementable values of ∆t.

B. analog quantum simulation
For an analog quantum computer, it is no longer possible to separate the unitary dynamics from the device noise.Instead.both of them act on the system simultaneously [20], namely, we now have The error mitigation superoperator M is thus required to be Using Baker-Campbell-Hausdorff formula, we can rewrite M as where • • • represents the terms coming from commutating relations that have higher orders of ∆t in the exponential.Note that in this case, only if L d − L n commutes with L h will M be Hamiltonian-independent without Trotter error, i.e., M is given by Eq. (10).This is different from the case of a digital quantum computer where the device noise L n does not contribute to the commutating relation.Thus one important consequence for an analog quantum computer is that, for simulations of closed dynamics L d = 0, if [L n , L h ] ̸ = 0, it is impossible to precisely recover the unitary dynamics in a Hamiltonianindependent way.This no-go theorem applies to any protocol aiming at cancelling the noise in a dynamical process, including different methods of quantum error mitigation, such as the method described in Ref. [20].

III. TOY MODELS
In this section we take single-qubit toy models to demonstrate the various simulation errors we have discussed in Sec.II.Although oversimplified, single-qubit models have certain advantages, a few of which are the following.The unitary dynamics can be expressed exactly as exp(−i Ĥt) without Trotter errors, so that we can focus on the Trotter errors related to the noise superoperators.Due to the small Hilbert space dimension, analytical solutions can be kept track of.The number of samples for realizing probabilistic error cancellation can be very large as the numerical cost per sample is very low.
For simplicity, we focus on Pauli noise models for both the device noise and the noise part of the open dynamics, so that the error mitigation superoperator Eq. ( 10) has a closed form expression [6].Specifically, the Lindblad operators Lk are proportional to the single qubit Pauli operators (see Appendix C).Moreover, we consider both the exact implementation of Eq. ( 10) and approximate implementations where Eq. ( 10) is expanded to the first order in ∆t.Due to the Trotter errors that we have discussed, applying M according to Eq. ( 10) for each time step may not lead to the accurate simulation result.We therefore use the fidelity to quantify the simulation errors.

A. probabilistic error cancellation
The quantum error mitigation method of probabilistic error cancellation [19] uses multiple copies of the quan-tum system to stochastically implement a non-physical map.For a single qubit as described above, the error mitigation superoperator Eq. ( 10) is in the following form, where q 0 +q 1 +q 2 +q 3 = 1 and q 0 > 0 (see Appendix C for details).As M might convert a time evolution to a less noisy version, it might be non-physical, corresponding to negative values of q 1 , q 2 and q 3 .In this case, Eq. ( 16) can only be implemented stochastically, in combination with post processing.To be specific, we first define sampling probabilities, where | • | refers to taking the absolute value.For each sample, a unitary operator is applied to the qubit, and a classical prefactor is multiplied to the density matrix.With probability 1 − µ 1 − µ 2 − µ 3 , the unitary operator is the identity operator 1 and the prefactor is ), the unitary operator is the Pauli operator Pi ( P1 = X, P2 = Ŷ , P3 = Ẑ) and the prefactor is sign where sign(•) refers to taking the sign.Averaging over an infinite number of samples corresponds to the ideal implementation of Eq. (16).For practical numerical simulations, however, there are finite errors that decrease with an increasing sample size.In Appendix B, we use numerical examples to show that deviations from the ideal sampling probabilities Eqs. ( 17)-( 19) may lead to non-physical simulation results, as the post-processing is not a physical process.

B. closed dynamics
The desired dynamics to be simulated are determined by the Hamiltonian Eq. ( 15).The time evolution of the population of the excited state can be expressed analytically independent of the parameter β, It represents an oscillation between 0 and 1 with period π/ω, as shown in the green dashed line in Fig. 1(a).

digital quantum simulation
We first consider a noisy digital quantum simulation of the closed dynamics.As discussed in Sec.II A, the noise  20).The orange dot-dashed line is the analytical result for the ideal implementation of the error mitigation maps, i.e., in the limit of an infinitely large sample size.(a) Exact probabilistic error cancellation.The parameters q0, q1, q2, q3 in Eq. ( 16) follow the exact inversion of the noise channel, as given in Eqs. ( 22)- (24).Here the green dashed line also corresponds to the ideal implementation of probabilistic error cancellation, therefore there is no systematic error.(b) Approximate probabilistic error cancellation.The parameters q0, q1, q2, q3 are the approximate ones following Eq.( 25).The orange dot-dashed line follows Eq. ( 26).It implies that even in the limit of an infinitely large number of samples, the noise cannot be fully cancelled.
is typically modelled by a separate layer of noise channel following the unitary layer within each time step.We assume that the noise channel is expressed as where the coefficients satisfy λ 1 , λ 2 , λ 3 ≥ 0, λ 1 +λ 2 +λ 3 ≤ 1, so that N (ρ) is a physical map.The error mitigation layer aims at cancelling the effect of N (ρ).Rigorous calculation steps involve finding L n expressed via λ 1 , λ 2 , λ 3 , and then calculating the error mitigation superoperator following Eq.( 9).This will be described in details in Sec.III C 1 when we consider the digital simulation of open dynamics.For the case here, we can directly take the inverse matrix of the superoperator form of N , and get M expressed in the form of Eq. ( 16) with coefficients These define the exact probabilistic error cancellation map.Note that, however, the calculation requires inverting a large matrix, thus not scalable to more qubits.Moreover, more advanced methods are required to generalize to noise models that are not Pauli [19,21].We may therefore also consider an approximate implementation of Eq. ( 9) if λ 1 , λ 2 , λ 3 ≪ 1, corresponding to negating the coefficients [10] in Eq. (21).
It also corresponds to expanding Eqs. ( 22)- (24) to the first order in λ 1 , λ 2 , λ 3 .Coefficients in Eq. ( 25) thus define the approximate probabilistic error cancellation map, which does not involve intensive calculations of matrix inversion.
As a numerical example, we consider β = 0 in the Hamiltonian Eq. ( 15), and take ω = 1.We choose the depolarizing noise channel, λ 1 = λ 2 = λ 3 = 0.05 in Eq. (21).We take the time step ∆t = 0.5 and 20 layers of error mitigation.The results are shown in Fig. 1.The numerical data come from sampling according to the probabilistic error cancellation procedure described in Sec.III A. The lines corresponding to the ideal implementations of probabilistic error cancellation, on the other hand, come from directly applying the superoperator Eq. ( 16) to the vectorized density matrix ρ, therefore equivalent to an infinite sample size.As discussed in Sec.II A, for the digital simulation of closed dynamics, Eq. ( 9) fully recovers the unitary dynamics without any Trotter error, as Trotter errors come from the commutator [L d , L h ] and for closed dynamics L d = 0. Thus in Fig. 1(a), the ideal implementation of the exact probabilistic error cancellation coincides with the actual closed dynamics, both of which are described by Eq. (20).The approximate probabilistic error cancellation plotted in Fig. 1(b), on the other hand, is based on a linear approximation to Eq. ( 9).This approximation induces error in the simulation, which persists even in the limit of an infinite number of samples.In fact, the orange dot-dashed line in Fig. 1(b) has an analytical closed-form expression due to the simplicity of the depolarizing noise channel we have chosen, where the time t = N ∆t is only defined at an integer multiple of ∆t, N = 0, 1, 2, • • • .It is interesting to look at the limit of taking an infinitesimal time step, ∆t → 0. We can consider two distinct situations.For the first situation, we assume that λ 1 is independent of ∆t.Taking the limit brings the right hand side of Eq. ( 26) to 1/2.Features in the dynamics are completely lost as we have introduced a fixed amount of noise per step and we have an infinite number of steps.For the second situation, we assume that λ 1 is proportional to ∆t.Taking the limit brings the right hand side of Eq. ( 26) to Eq. ( 20), which is the closed dynamics.Details are included in Appendix A.

analog quantum simulation
Now we consider a noisy analog quantum simulation of the closed dynamics.In this case, the device noise acts simultaneously with the unitary dynamics to be simulated.As discussed in Sec.II B, the Hamiltonian-independent error mitigation superoperator M is defined according to Eq. ( 9).However, there is no Trotter error only if [L h , L n ] = 0. We are going to consider two noise models.The first noise model is a depolarizing noise model, whose superoperator commutes with the superoperator of the Hamiltonian.The second noise model only contains Pauli X operator, therefore only commutes with L h if the parameter β in the Hamiltonian Eq. ( 15) is taken as β = π/2.Depolarizing noise.Under a depolarizing noise model, the quantum state evolves according to the master equation where the Hamiltonian Ĥ is given in Eq. (15).The solution can be found to be representing a Rabi oscillation with damping.In order to apply probabilistic error cancellation, we discretize the continuous dynamics into time step ∆t, and after each time step, apply the error mitigation step following Eq.( 9) and the procedure described in Sec.III A. Note samples, and the errorbars correspond to the standard deviation.The orange dot-dashed line is the analytical result for the ideal implementation of the error mitigation maps, i.e., in the limit of an infinitely large sample size.(a) Exact probabilistic error cancellation.The parameters q0, q1, q2, q3 in Eq. ( 16) correspond to the exact probabilistic error cancellation map, as given in Eq. (30).Here the orange dot-dashed line also corresponds to the actual closed dynamics to be simulated, as [Ln, L h ] = 0 implies no Trotter error.(b) Approximate probabilistic error cancellation.The parameters q0, q1, q2, q3 are the approximate ones following Eq.( 31).The orange dot-dashed line follows Eq. (32).It implies that even in the limit of an infinitely large number of samples, the result is different from the actual closed dynamics.
that the device noise superoperator, commutes with L h , therefore there is no Trotter error.As there are only Pauli errors, exp(±L n ∆t) have simple closed form expressions [6].The error mitigation operator M is therefore expressed in the form of Eq. ( 16) with actual dynamics popula�on FIG. 3. Combining step-wise probabilistic error cancellation with a noisy analog quantum computer to simulate closed dynamics, assuming a Pauli-X noise channel.We have chosen the parameters as ω = 1, ∆t = 0.5, κ = 0.3.The numerical data points are the ensemble averages over 5 × 10 6 samples, and the errorbars correspond to the standard deviations.The green dashed line is the actual closed dynamics to be simulated, as in Eq. ( 20).The parameters q0, q1, q2, q3 in Eq. ( 16) correspond to the exact probabilistic error cancellation map, as given in Eq. ( 34).The brown square data points for β = π/2 correspond to the case without Trotter error.Decreasing β to π/4 (purple triangle) and 0 (red circle) leads to larger Trotter errors resulting from the non-zero [Ln, L h ]. coefficients In addition to this exact map, we can also define an approximate error mitigation map by keeping only the linear term in ∆t, although taking the approximation at different places lead to different expressions.Here we expand exp(L n ∆t) to the first order in ∆t, and define the approximate M as the inverse of the expansion of exp(L n ∆t), resulting in We are going to introduce another way of linear expansion when we deal with the simulation of open dynamics in Sec.III C 2. In Fig. 2(a) and (b), we show the numerical simulations following the coefficients in Eq. (30) and Eq.(31), respectively.Note that the results are independent of the parameter β in the Hamiltonian.While the exact coefficients Eq. (30) result in recovering the actual closed dynamics to be simulated in the limit of an infinite number of samples, the approximate coefficients Eq. (31) lead to errors that increase in time.Note that these errors are not the Trotter error mentioned before.In fact, the dotdashed orange line in Fig. 2(b) has an analytical form, popula�on actual dynamics FIG. 4. Combining step-wise probabilistic error cancellation with a noisy analog quantum computer to simulate closed dynamics, assuming a Pauli-X noise channel.We have chosen the parameters as ω = 1, ∆t = 0.5, κ = 0.3.The numerical data points are the ensemble averages over 5 × 10 6 samples, and the errorbars correspond to the standard deviations.The green dashed line is the actual closed dynamics to be simulated, as in Eq. ( 20).The parameters q0, q1, q2, q3 in Eq. ( 16) correspond to the approximate probabilistic error cancellation map, as given in Eq. ( 35).Increasing β from 0 (red circle) to π/4 (purple triangle) to π/2 (brown square) leads to increasing simulation errors.Instead of Trotter errors, these errors are from the approximation of linear expansion, and can lead to non-physical results.
defined for t = N ∆t with N = 0, 1, 2, • • • .This corresponds to the solution of a master equation, Eq. ( 27) replacing κ with Note that this is smaller than 0, indicating a non-physical result.This is because the process of probabilistic error cancellation as described in Sec.III A involves post processing that does not satisfy physical rules.Pauli-X noise.Next we consider a noise model that does not commute with the Hamiltonian Eq. ( 15) for general β.The evolution of the system follows In this case, L n = κ( X ⊗ X − 1 ⊗ 1) does not commute with L h unless β = π/2.This also implies that a succinct form of ⟨1|ρ(t)|1⟩ no longer exists.For the exact probabilistic error cancellation map, we still define M according to Eq. ( 9).This corresponds to the coefficients in Eq. ( 16) as Numerical simulations are shown in Fig. 3.Although not large, the deviation of the error mitigated results We have chosen the parameters as ω = 1, ∆t = 0.5, λ1 = 0.16, λ2 = 0.12, λ3 = 0.2, γ = 0.3.The fidelity is defined by comparing the exact solution of Eq. ( 36) with the result from the ideal implementation (i.e., infinite sample size) of the exact error mitigation map with coefficients Eq. ( 39).The red solid line with circle marks is for β = 0.The purple dashed line with triangle marks is for β = π/4.The brown dot-dashed line with square marks is for β = π/2, in which case there is no Trotter error, so that the fidelity remains 1.
compared with the actual closed dynamics increases as β decreases from π/2 to 0. This deviation comes from the non-commuting L h and L n .We can also define the approximate probabilistic error cancellation map, using the linear expansion method described above in the part of depolarizing noise.The resulting coefficients are Numerical results are shown in Fig. 4. The deviation of the error mitigated results compared with the actual closed dynamics is much larger than the case of Fig. 3.This is because of the linear expansion we made, which causes much larger errors than the Trotter errors.

C. open dynamics
Now we will demonstrate the quantum simulation of open dynamics, where the Trotter errors nontrivially depend on the superoperators of the noise part of the open dynamics to be simulated (L d ), the unitary part (L h ), and the device noise (L n ), as argued in Sec.II.Specifically, we consider simulating the following master equation, describing Rabi oscillations with Pauli-X error continuous in time.Same as before, the Hamiltonian Ĥ is given by Eq. ( 15).As L d = γ( X ⊗ X − 1 ⊗ 1) only commutes  36).The orange dot-dashed line is the analytical result (i.e., in the limit of an infinitely large sample size, we have also called it "ideal implementation") for the approximate error mitigation maps.The parameters q0, q1, q2, q3 in Eq. ( 16) follow Eq. ( 42).The difference between the green dashed line and the orange dot-dashed line implies that due to the approximation to the first order of γ∆t, λ1, λ2 and λ3, even in the limit of an infinitely large number of samples, the noise cannot be fully cancelled.(a) with L h for β = π/2, in general the exact solution as expressed by Eq. ( 5) does not have a simple analytical form.Instead, we plot the exact solutions in the green dashed lines in Fig. 6, for three values of β.The time evolution of the exited state population is a damped Rabi oscillation, while β = 0 corresponds to the weakest damping and β = π/2 corresponds to the strongest damping.The numerical data points are the ensemble averages over 5 × 10 6 samples, and the errorbars correspond to the standard deviations.The parameters q0, q1, q2, q3 in Eq. ( 16) correspond to the exact probabilistic error cancellation map, as given in Eq. (30).The brown square are data points for β = π/2, the purple triangles are for β = π/4, and the red circles are for β = 0.

digital quantum simulation
We have shown in Sec.II A that, if [L d , L h ] ̸ = 0, there will be Trotter errors in the quantum simulation.Importantly, this criterion is independent of the device noise superoperator L n .We choose the device noise channel to follow the same form as in Eq. ( 21).Its vectorized form, can be rewritten in the form of exp(L n ∆t) in order to match Eq. ( 7), with and coefficients Here we have assumed that the device noise channel is weak, e.g., λ 1 +λ 2 +λ 3 ≤ 1/2.The exact error mitigation superoperator, as defined in Eq. (10), is in the form of Eq. ( 16) with coefficients where k = 1, 2, 3 and δ i,j is the Kronecker delta.In Fig. 5, we show numerical examples of how the commutator [L d , L h ] affects the Trotter errors in the quantum simulation.As we aim at demonstrating the Trotter errors, which are rather small (as clearly seen in the case of Fig. 3), we only show the ideal implementation of the error mitigation map, i.e., considering an infinite number of samples.We quantify how close the simulation result is to the exact open dynamics via fidelity, which for qubits is expressed as [22] where det(•) refers to the determinant.As seen in Fig. 5, only the case of β = π/2 corresponds to simulation results without any error, as for this β, the Hamiltonian only contains the Pauli X operator.We can also define the approximate error mitigation superoperator by expanding Eq. (38) to the first order in λ 1 , λ 2 , λ 3 , The approximate error mitigation superoperator thus has coefficients Note that this is the approach taken by Ref. [10].In Fig. 6, we show numerical examples of how the approximation described above affect the simulation results.Although getting the approximated expression of FIG. 9. Combining step-wise probabilistic error cancellation with a noisy analog quantum computer to simulate open dynamics, assuming a biased Pauli-X noise channel.We have chosen the parameters as ω = 1, ∆t = 0.5, κ = 0.1, γ = 0.3.The fidelity is defined by comparing the exact solution of Eq. ( 36) with the result the ideal implementation of the approximate error mitigation map with coefficients Eq. ( 48).The red solid line with circle marks is for β = 0.The purple dashed line with triangle marks is for β = π/4.The brown dot-dashed line with square marks is for β = π/2.the error mitigation map M following Eq.( 42) does not require inverting any matrix, the resulting simulation error turns out to be significantly larger than the Trotter error.Most oscillatory features of the damped Rabi oscillation cannot be resolved.We will look into two noise models.The first one is a depolarizing noise model, which will induce Trotter errors.The second one is a biased Pauli-X noise model, which has more Pauli-X noise than a depolarizing noise model.As the open dynamics we aim at simulating only has Pauli-X noise, this noise model will not have Trotter error.Depolarizing noise.For the depolarizing noise channel, the system evolves following the master equation Eq. ( 27).Step-wise probabilistic error cancellation implements the map Eq. (10).In the form of Eq. ( 16), the coefficients are given by (43) Interestingly, the simulation results in this case are exactly the same as those shown in Fig. 5.To prove this, we focus on the superoperator applied to the vectorized state at each time step.For the case in Fig. 5, the superoperator is e (L d −Ln)∆t e Ln∆t e L h ∆t = e L d ∆t e L h ∆t . (44) Here the effective L n is defined via Eq.(37) and Eq. ( 38), and we have used the fact that both L d and L n only contain Pauli terms, therefore they commute with each other.For the analog simulation case considered here, the superoperator at each time step is Here L ′ n is the depolarizing noise channel as written in Eq. ( 29), and we have used that the depolarizing noise commutes both with L d and with L h .As the open dynamics to be simulated corresponds to the step-wise superoperator exp(L d ∆t + L h ∆t), Trotter error exists, as illustrated in Fig. 5.
We can also define the approximate error mitigation map, by expanding Eq. ( 43) to the first order in κ∆t and γ∆t.This results in the coefficients The numerical results are shown in Fig. 7.Although the fidelities are close to 1, the errors are not just Trotter errors, but also the errors from the linear expansion of q 1 , q 2 and q 3 .Compared with the case in Fig. 6, the errors here are much smaller, as there is no step similar to the approximation in Eq. ( 41).Interestingly, the errors here are even smaller than the case of the exact error mitigation map, Fig. 5, possibly because the Trotter errors and the linear approximation errors cancel with each other.Biased Pauli-X noise.Now we consider a biased Pauli-X noise model, such that [L h , L d − L n ] = 0 is satisfied.We assume that the system evolves following the master equation The exact step-wise probabilistic error cancellation superoperator Eq. (10), written in the form of Eq. ( 16), has coefficients that are same as Eq.(30).Numerical results are shown in Fig. 8.As there is no Trotter error for all values of β, the numerical data match the analytical actual dynamics in Fig. 6.
For the approximate error mitigation map, we can expand Eq. (30) to the first order in κ∆t and get The numerical results are shown in Fig. 9.The small deviations of the fidelity from 1 are results of the linear approximations we have made, rather than the Trotter errors.

IV. CONCLUSION AND DISCUSSION
We have analyzed the limitations of combining stepwise probabilistic error cancellation with noisy quantum simulations of continuous dynamics beyond the exponentially large sampling overhead.We have pointed out the Trotter errors originating from the superoperators governing the unitary dynamics L h , the device noise L n , and the noise part L d of the open dynamics to be simulated.These Trotter errors exist even in the limit of an infinite sample size.For a digital quantum simulation, a nonzero commutator [L d , L h ] leads to the Trotter error.For an analog quantum simulation, a non-zero commutator [L d − L n , L h ] leads to the Trotter error.Importantly, even for simulating closed dynamics, L d = 0, this Trotter error can still exist for an analog quantum simulation.We have also pointed out that, the commonly used coefficients in the error mitigation map, which do not involve inverting a large matrix, are actually linear approximations of the exact error mitigation map, therefore may induce additional simulation errors.Although dominated by the numerical errors in real quantum devices for now, these errors that we have investigated here put fundamental limitations on the constructions of the simulation methods.It can be expected that as the qualities of the quantum hardware improve, these errors are going to become limiting factors for quantum simulations.
We note that, however, a systematic study as what we have done here is very difficult to go beyond few-qubit toy models and simple noise channels such as Pauli noise channels.How to reduce the Trotter errors in a scalable way remains an open question.One possible way is to design the simulation setup such that, the Trotter error from the non-commuting Hamiltonians of a multi-qubit system partially cancels the Trotter error between the unitary part and the noise part of the dynamics to be simulated.Another possible way may be to remove the time discretization as a middle step, and instead find methods to map one dissipative dynamics to another in a matter that is continuous in time [23][24][25], that can be implemented under noisy quantum operations.ing with the noise intrinsic to the quantum computer.As a simple example, we suppose that the noise in the quantum computer effectively applies a depolarizing channel after each quantum gate is applied [10], and the strength of the depolarizing noise channel is proportional to the time-step ∆t.To be specific, in each step the density matrix is transformed in the following way, Here ρ 1 denotes the resulting density matrix for this digital simulation, and κ∆t ≪ 1 corresponds to the error rate in the depolarizing channel.Note that we have included the factor ∆t in the definition of the error rate to facilitate a direct comparison in the dissipation rate with Eq. ( 27).Starting from the initial state ⟨1|ρ(0)|1⟩ = 1, the simulation result is derived as (A2) This corresponds to the solution of the master equation ( 27) with a time-step-dependent dissipation rate In the limit ∆t → 0, Eq. (A1) is equivalent to Eq. ( 27), which is also an algebraic way [26] of deriving Eq. (28).However, for finite ∆t, Eq. (A1) always corresponds to a stronger dissipation than Eq. ( 27).As shown in Fig. 10, although the steady state is not changed by this difference, the time evolution before reaching the steady state is changed.samples, and the errorbars correspond to the standard deviation.The orange lines are the analytical results in the limit of an infinitely large sample size, Eq. (B2).(a) In each error mitigation step, the probability µ ′ of sampling the X (or Ŷ or Ẑ) is smaller than the ideal probability µ1, µ ′ = 0.97µ1.This corresponds to κ ′ ≈ 0.0041 and ξ ≈ 1.01, meaning that the trace is increasing with the number of steps.(b) The probability µ ′ of sampling the X (or Ŷ or Ẑ) is larger than the ideal probability µ1, µ ′ = 1.03µ1.This corresponds to κ ′ ≈ −0.0042 and ξ ≈ 0.99, meaning that the trace is decreasing with the number of steps, but instead of a decay, the oscillation amplitude increase in time.

Appendix B: imprecise sampling probabilities
In principle, adding probabilistic error cancellation after every time-step in a digital quantum simulation can remove the noise of the quantum computer in a statistical way.However, it has been known that the major disadvantage of quantum error mitigation is that the sampling cost scales exponentially as the number of quantum gates to be applied [2,11].Here we will point out another issue: When randomly sampling Pauli gates to be applied to the quantum state, deviation in the sampling probabilities may lead to non-physical simulation results.
The map in Eq. ( 16) is non-physical if not all q 1 , q 2 and q 3 are non-negative, as shown in the examples in Sec.III.
Therefore M can only be implemented stochastically, together with post-processing, following the description in Sec.III A. We consider a case where the sampling probabilities are deviated from the ideal values µ 1 , µ 2 , µ 3 , but the prefactors are not modified.This can be due to an unknown bias in the sampling process.We take the example of Eq. (A1), and aim at removing the noise characterized by κ∆t.In order to simplify the expressions, we assume that the probability of applying the X (or Ŷ or Ẑ) gate is µ ′ , which is different from µ 1 (= µ 2 = µ 3 ).
The implemented map M(ρ) can be not trace-preserving as we have Applying the error mitigation map M(ρ) after each step of the noisy digital quantum simulation Eq. (A1) results in the statistically averaged simulation result ⟨1|ρ(N ∆t)|1⟩ = ξ N 1 2 1 + e −4κ ′ N ∆t cos(2ωN ∆t) , (B2) with In addition to the change in the trace described by ξ, another non-physical feature is that κ ′ can be negative.The reason of simulation results being non-physical is that, even though the process of randomly sampling a Pauli matrix to apply it to the state is a physical operation that preserves the trace and makes sure that the diagonal elements of the density matrix are non-negative, the post-processing of multiplying the density matrix with prefactor is not a physical process.Numerical examples are shown in Fig. 11.Although ξ only slightly deviates from 1 and if κ ′ < 0 it is very close to 0, the exponential in N means that the non-physical features grow exponentially in number of steps, i.e., circuit depth.

FIG. 1 .
FIG.1.Combining step-wise probabilistic error cancellation with a noisy digital quantum computer to simulate closed dynamics.We have chosen the parameters as ω = 1, ∆t = 0.5, λ1 = λ2 = λ3 = 0.05.The blue numerical data points are the ensemble averages over 2 × 10 7 samples, and the errorbars correspond to the standard deviation.The green dashed line is the actual closed dynamics to be simulated, following Eq.(20).The orange dot-dashed line is the analytical result for the ideal implementation of the error mitigation maps, i.e., in the limit of an infinitely large sample size.(a) Exact probabilistic error cancellation.The parameters q0, q1, q2, q3 in Eq. (16) follow the exact inversion of the noise channel, as given in Eqs.(22)-(24).Here the green dashed line also corresponds to the ideal implementation of probabilistic error cancellation, therefore there is no systematic error.(b) Approximate probabilistic error cancellation.The parameters q0, q1, q2, q3 are the approximate ones following Eq.(25).The orange dot-dashed line follows Eq. (26).It implies that even in the limit of an infinitely large number of samples, the noise cannot be fully cancelled.

FIG. 2 .
FIG. 2. Combining step-wise probabilistic error cancellationwith a noisy analog quantum computer to simulate closed dynamics, assuming a depolarizing noise channel.We have chosen the parameters as ω = 1, ∆t = 0.5, κ = 0.1.The blue numerical data points are the ensemble averages over 5 × 10 6 samples, and the errorbars correspond to the standard deviation.The orange dot-dashed line is the analytical result for the ideal implementation of the error mitigation maps, i.e., in the limit of an infinitely large sample size.(a) Exact probabilistic error cancellation.The parameters q0, q1, q2, q3 in Eq. (16) correspond to the exact probabilistic error cancellation map, as given in Eq. (30).Here the orange dot-dashed line also corresponds to the actual closed dynamics to be simulated, as [Ln, L h ] = 0 implies no Trotter error.(b) Approximate probabilistic error cancellation.The parameters q0, q1, q2, q3 are the approximate ones following Eq.(31).The orange dot-dashed line follows Eq. (32).It implies that even in the limit of an infinitely large number of samples, the result is different from the actual closed dynamics.

FIG. 5 .
FIG.5.Combining step-wise probabilistic error cancellation with a noisy digital quantum computer to simulate open dynamics.We have chosen the parameters as ω = 1, ∆t = 0.5, λ1 = 0.16, λ2 = 0.12, λ3 = 0.2, γ = 0.3.The fidelity is defined by comparing the exact solution of Eq. (36) with the result from the ideal implementation (i.e., infinite sample size) of the exact error mitigation map with coefficients Eq. (39).The red solid line with circle marks is for β = 0.The purple dashed line with triangle marks is for β = π/4.The brown dot-dashed line with square marks is for β = π/2, in which case there is no Trotter error, so that the fidelity remains 1.

FIG. 6 .
FIG.6.Combining step-wise probabilistic error cancellation with a noisy digital quantum computer to simulate open dynamics.We have chosen the parameters as ω = 1, ∆t = 0.5, λ1 = 0.16, λ2 = 0.12, λ3 = 0.2, γ = 0.3.The green dashed line is the actual open dynamics to be simulated, which is the solution of Eq. (36).The orange dot-dashed line is the analytical result (i.e., in the limit of an infinitely large sample size, we have also called it "ideal implementation") for the approximate error mitigation maps.The parameters q0, q1, q2, q3 in Eq. (16) follow Eq.(42).The difference between the green dashed line and the orange dot-dashed line implies that due to the approximation to the first order of γ∆t, λ1, λ2 and λ3, even in the limit of an infinitely large number of samples, the noise cannot be fully cancelled.(a) β = 0. (b) β = π/4.(c) β = π/2.

FIG. 7 .FIG. 8 .
FIG.7.Combining step-wise probabilistic error cancellation with a noisy analog quantum computer to simulate open dynamics, assuming a depolarizing noise channel.We have chosen the parameters as ω = 1, ∆t = 0.5, κ = 0.1, γ = 0.3.The fidelity is defined by comparing the exact solution of Eq. (36) with the result from the ideal implementation of the approximate error mitigation map with coefficients Eq. (46).The red solid line with circle marks is for β = 0.The purple dashed line with triangle marks is for β = π/4.The brown dot-dashed line with square marks is for β = π/2.

2 .
analog quantum simulation Next we focus on noisy analog quantum simulations of the open dynamics.As derived in Sec.II B, the Trotter error now depends on the commutator [L h , L d − L n ].

FIG. 11 .
FIG.11.Using probabilistic error cancellation to remove the noise, considering imprecise sampling probabilities.We have chosen the parameters as ω = 1, κ = 0.1, ∆t = 0.5.The blue numerical data points are the ensemble averages over 2 × 10 7 samples, and the errorbars correspond to the standard deviation.The orange lines are the analytical results in the limit of an infinitely large sample size, Eq. (B2).(a) In each error mitigation step, the probability µ ′ of sampling the X (or Ŷ