Quantum Error Mitigation via Quantum-Noise-Effect Circuit Groups

Near-term quantum computers have been built as intermediate-scale quantum devices and are fragile against quantum noise effects, namely, NISQ devices. Traditional quantum-error-correcting codes are not implemented on such devices and to perform quantum computation in good accuracy with these machines we need to develop alternative approaches for mitigating quantum computational errors. In this work, we propose quantum error mitigation (QEM) scheme for quantum computational errors which occur due to couplings with environments during gate operations, i.e., decoherence. To establish our QEM scheme, first we estimate the quantum noise effects on single-qubit states and represent them as groups of quantum circuits, namely, quantum-noise-effect circuit groups. Then our QEM scheme is conducted by subtracting expectation values generated by the quantum-noise-effect circuit groups from that obtained by the quantum circuits for the quantum algorithms under consideration. As a result, the quantum noise effects are reduced, and we obtain approximately the ideal expectation values via the quantum-noise-effect circuit groups and the numbers of elementary quantum circuits composing them scale polynomial with respect to the products of the depths of quantum algorithms and the numbers of register bits. To numerically demonstrate the validity of our QEM scheme, we run noisy quantum simulations of qubits under amplitude damping effects for four types of quantum algorithms. Furthermore, we implement our QEM scheme on IBM Q Experience processors and examine its efficacy. Consequently, the validity of our scheme is verified via both the quantum simulations and the quantum computations on the real quantum devices.

While the above successful results of the research and development of quantum computers have been reported, nearterm quantum computers based on circuit models have been built as intermediate-scale quantum devices yet and are fragile against quantum noise effects: they are called noisy intermediate-scale quantum (NISQ) devices [10,12,38].Quantum noise effects (decoherence) are major obstacles for performing quantum computation and historically many great efforts have been made on reducing such effects [39,40].One of the traditional and representative schemes for this is the quantum-error-correction (QEC) coding [5,8,10,11,17,[41][42][43][44].Another important one is the dynamical decoupling which plays fundamental role in extending coherence times of qubits [9,12,17,43,[45][46][47][48].The QEC codes are, however, not implemented on NISQ devices and to obtain quantum computational results in good accuracy with NISQ devices we need to search for alternative approaches for mitigating quantum noise effects.This research field is called quantum error mitigation (QEM), and these days, it is one of the important themes of the research and development of quantum computation [26,27,.The difficulty of the treatment of quantum noises (e.g., amplitude damping, phase damping (dephasing), depolarizing channel) is that we cannot directly construct their inverse processes by quantum gates due to their non-unitarity.On the other hand, it is possible to formulate quantum noise effects as quantum circuits by using ancilla bits and measurements on them [5,[78][79][80][81][82][83][84][85].By utilizing the quantum circuits representing the quantum noise effects under consideration, we expect that we can establish QEM schemes for reducing such effects.If this is established, we become able to mitigate the quantum noise effects by the gate operations and measurements, i.e., QEM conducted by all-quantum-computational operations.In other words, we become able to programmably run quantum algorithms with mitigating the quantum noise effects solely by the quantum computational operations and realize high-accurate quantum computation.
In this work, we propose our QEM schemes for quantum computational errors which occur owing to couplings with environments (decoherence) during gate operations: errors of state preparation (initialization) and measurement, imperfections of quantum gates, and cross talks among qubits are not taken into account.In particular, we make detailed analysis on quantum computational errors generated by amplitude damping (AD) of single-qubit states.We show the schematic representation of our QEM scheme in Fig. 1 and it consists of two elements, the quantum circuit for the quantum algorithm under consideration (original circuit) represented by the blue rectangle and the ensemble of quantum circuits which yields the theoretical value of the quantum computational error due to the quantum noise effect, namely, quantum-noise-effect circuit group and is represented by the orange rectangles.By utilizing the quantum-noise-effect circuit groups, we formulate our QEM scheme as a perturbation theory with respect to a strength(s) of quantum noise(s) and perform it by subtracting the expectation values given by the quantum-noise-effect circuit groups from those generated by the quantum circuits for the quantum algorithms under consideration as expressed by the formula in the green rectangle; see also the right-hand side of the first line in Eq. (8).As a result, the quantum noise effects are mitigated and we approximately obtain the ideal expectation values.Then, we discuss the numbers of elementary quantum circuits which compose the quantum-noiseeffect circuit groups and show that they scale polynomial (linear) with respect to the products of the numbers of register bits and the depths of the quantum algorithms (circuit depths or the numbers of unitary gates composing the quantum algorithm).Finally, we numerically demonstrate the validity of our QEM scheme by running noisy quantum simulators of qubits under the AD effects for four types of quantum algorithms in the linear-order perturbation regime.The detailed explanation on how to extend our QEM scheme to other kinds of quantum noise effects including phase damping, generalized amplitude damping (thermalization), and depolarizing channel, and extension of our QEM scheme to higher-order quantum noise effects are given in Supplementary Material.
The structure of this paper is given as follows.It begins by Sec.III with our modeling of the quantum computation under the influence of the quantum noise effect.After then we explain the formalisms of our QEM scheme.In Sec.IV, which presents our main results, we demonstrate numerically our QEM schemes for the noisy quantum simulations for four types of quantum algorithms.These simulations are done by both our original numerical code and Qiskit [86].In Sec.V, we discuss our quantum computation results for our QEM scheme run on the IBM Q Experience processors [87].In Sec.VI, we make comparisons between our scheme and other QEM methods.Sec.VII is devoted to the conclusion of this paper.

III. QEM SCHEMES A. Modeling and Formulation
Let us explain our modeling of quantum computation under the influence of quantum noise effects [26,27,49,60,61,86].In the following, we focus on the amplitude damping (AD) effect: generalized-amplitude-damping (GAD) effect at zero temperature.As discussed later, it is straightforward to generalize the argument for the AD effect to other quantum noise effects such as phase damping (PD) and stochastic Pauli noises.We schematically represent such a circumstance as a quantum circuit and show it in Fig. 2.There are N q register bits and the quantum algorithm to be run is represented by the unitary transformation U QC .It is comprised of d unitary transformations described by The unitary transformation U k (k = 1, 2, . . ., d) is composed of single-and two-qubit gates.We assume that the duration time (gate operation time) of the unitary transformation U k is ∆t for any k.During the time interval ∆t, the register bits are influenced by the AD effects due to couplings with environments, e.g., electromagnetic field in the vacuum, phonons in solids, etc.The quantum master equation describing the AD process in the interaction picture is given by [5,[88][89][90] where ρ(t) is the density matrix of the N q register bits at a time t and γ is the decay rate.The symbol L AD denotes the Lindblad superoperator of the AD process and the operators σ± j = Xj ∓iYj 2 are the ladder (raising and lowering) operators acting on the register bit Q rj .X j and Y j are X and Y gates acting on Q rj , respectively.A, B is the anti-commutator between the operators A and B. In our model, we assume that the N q register bits experience homogeneously the AD effect of single-qubit state given by the decay rate γ.At the initial time t = 0, all the register bits are in |0⟩ state (ground state), namely, ρ(0) = |0⟩⟨0| ⊗Nq .Let us write the total amount of quantum computational time (running time of the quantum algorithm under consideration) by T (= d • ∆t) while we introduce the dimensionless time τ = γ∆t.By assuming τ ≪ 1, in the following let us evaluate the density matrix at the time T , ρ(T ), by using the quantum master equation (1) and express it as a perturbation series with respect to τ given by Here • U QC † describes the noise-free (ideal) quantum state of the register bits.In other words, it is the ideal output quantum state generated by the quantum algorithm given by U QC .The quantity .On the other hand, the quantum-noise-effect circuit group, which is represented by the orange rectangles (right side), is constructed from the original circuit by inserting an additional operation between U k and U k+1 (gray box).It yields the theoretically-estimated quantum computational error ⟨ Ô⟩ (∆ AD 1 ρ d...1 ) real .By using these two expectation values, we obtain the equation for our QEM scheme in the green rectangle.Here we have taken d = 3.The symbol E AD Q rj expresses the occurrence of AD effect on the register bit Qrj (the j-th register bit).
the theoretically-evaluated p-th-order AD effect.Let us focus on the first-order AD effect ∆ AD 1 ρ d•••1 which has the form where with In the above equation we have used describes the projection onto the quantum state |1⟩ j with 1 j and Z j denoting the identity operator and the Z gate acting on Q rj , respectively: On the other hand, the projection operator of the quantum state |0⟩ j is given by P .

B. QEM Scheme
Since we have evaluated the single-qubit-state AD effect, next we discuss our quantum error mitigation (QEM) scheme.We denote the operator of which we are aiming to take an expectation value by Ô.When we implement the quantum state ρ on a real device what we actually obtain is a quantum state which is different from ρ due to quantum noise effects: Note again that hereinafter we only consider the AD effect.Let us write it by ρ real .We represent the density matrix ρ real in terms of ρ (ideal state) as ρ real = ρ + δ AD ρ, where δ AD ρ represents the deviation from ρ owing to the AD effect on a real device.Note that we use the symbol δ AD to describe the AD effect on a real device while we use ∆ AD to describe the theoreticallyestimated AD effect like Eq. ( 2).Namely, a quantum computational error occurs due to the deviation δ AD ρ.QEM is a prescription for mitigating the error coming from the deviation δ AD ρ.Mathematically, this is a task to make the value of Tr( Ôδ AD ρ) as small as possible.In our scheme, we mitigate the error Tr( Ôδ AD ρ) by perturbatively treating the deviation δ AD ρ with respect to τ and using the theoretically-estimated AD effect ∆ AD p ρ.In the following we show such a perturbative analysis up to the first order in τ .The extension of QEM scheme to higher-order AD effect is discussed in Sec.I in the Supplementary Material.The key procedure of our QEM scheme is to construct quantum circuits for computing the quantity Tr( Ô∆ AD 1 ρ d•••1 ), which describes the theoreticallyestimated quantum computational error of the expectation value Tr( Ôρ) in the first order of τ .For doing this, there are two difficulties: (i) the generation of the anti-commutator term P 1 j , ρ k•••1 in Eq. ( 4) and (ii) the implementation of the non-unitary operators σ− j and P 1 j .Let us discuss from our solution to the difficulty (i).We denote some sort of quantumcomputational operation (gate operation or measurement) by A. When the operation A acts on the quantum state , in contrast, is not represented in this way, and thus, it is not clear how to generate such a term by the quantum-computational operations.We solve this in the following way.To make our argument simple, here let focus on the single-register-bit system (N q = 1); the generalization to N q ≥ 2 is straightforward and is discussed later.First, we rewrite ρAD k•••1 in Eq. ( 4) as In the above way, all the four terms in Eq. ( 5) are written in the form Aρ k•••1 A † , and thus, we have solved the difficulty (i).Let us analyze the mathematical structure of the righthand side of Eq. ( 5).The quantum circuit for creating the first term is straightforward because it is obtained by the quantum circuit composed of U QC (the quantum algorithm under consideration).The implementation of the quantum circuit for the second term    is also straightforward because we just apply the Z gate after the operation of U k .The unclear part is to find ways to construct the quantum circuits for generating the third and fourth terms given by the non-unitary operators σ− and P 1 and this is nothing but the difficulty (ii).We solve this by using an ancilla bit and a measurement on it [5].For the creating the operation σ− , we use the quantum circuit presented in Fig. 3(a) (AD-effect circuit A) while for the operation of P 1 we use the one in Fig. 3(b) (AD-effect circuit B).In both quantum circuits, we regard the ancilla bit Q a as the environment which induces the AD effect on the register bit Q r0 .The interactions between these two qubits are represented by the controlled-rotational gate describes the rotation about y axis by the rotational angle ϑ and it is composed of the control bit Q r0 and the target bit Q a .On the other hand, for the gate operation U CX [Q a ; Q r0 ] the ancilla bit Q a is the control bit while the register bit Q r0 is the target bit.We have used the notation such that the control bit(s) comes before the semicolon while the target bit(s) comes after.Let us explain the output states generated by the AD-effect circuits A and B. For both quantum circuits, the initial quantum states of Q r0 and Q a are the same and it is ρ Qr0Qa (0) = ρ Qr0 (0) ⊗ ρ Qa (0) with ρ Qr0 (0) = |0⟩ Qr0 ⟨0| and ρ Qa (0) = |0⟩ Qa ⟨0|.The AD-effect circuit A is given by the unitary operation , where 1 2×2 is the two by two identity operator.Owing to these unitary operations, the quantum state generated by the AD-effect circuit A is given by At the end, we measure the ancilla bit Q a .Then the quantum states of Q r0 (reduced density matrices) are described by the Kraus operators [5,91].
For the AD-effect circuit A (B) the Kraus operators ) acts on the register bit Q r0 when the measurement outcome of the quantum state of the ancilla bit ) operates when the measurement outcome is |1⟩ Qa .When we average these two outcomes, the quantum state of Q r0 created by the AD-effect circuit A is given byρ , where the symbol Tr Qa denotes the partial trace with respect to Q a degrees of freedom.Similarly, for the AD-effect circuit B we have In particular, for the AD-effect circuit A when we take ϑ to be ϑ t such that cos 2 ϑt 2 = e −γt [5], the matrix representation of ρ ADA,Qr0 (ϑ t ) is given by The matrix element [ρ ADA,Qr0 (0)] nn ′ (n, n ′ = 0, 1) denotes the (n, n ′ )-element of ρ ADA,Qr0 (0).The reduced density matrix ρ ADA,Qr0 (ϑ t ) in Eq. ( 7) is nothing but the solution of the quantum master equation (1).Further, when we take ϑ t → π, the Kraus operators in Eq. ( 6) becomes Therefore, for the case of the AD-effect circuit A by using the measurement result of Q a such that we post-select the output state of Q a to be |1⟩ Qa we can realize the operation of σ− on Q r0 .On the other hand, for the AD-effect circuit B by post-selecting the output state of Q a to be |0⟩ Qa we realize the operation of P 1 on Q r0 .To show the above things concretely, let us present the examples of the quantum circuits for the generation of 3, and we show them in Figs.4(a) and 4(b), respectively.As a result, by using the AD-effect circuits A and B we can perform the actions of σ− and P 1 as described by the third and four terms in Eq. ( 5), and thus, we have solved the second difficulty (ii).Since the difficulties (i) and (ii) have been solved, we are now ready to establish our QEM scheme.To compute the quantity Tr( we need four types of quantum circuits: the original circuit given by U QC , the quantum circuit where the additional Z gate is performed, and the AD-effect circuits A and B. The latter three quantum-circuit ensembles composed of the additional Z-gate, σ− , and P 1 operations form the quantum-noise-effect circuit group for the AD effect, i.e, AD-effect circuit group.Hereinafter, let us write the trace of the product between the operator Ô and ρ by Tr( Ôρ) = ⟨ Ô⟩ ρ .We perturbatively express the quantum states . Furthermore, we write the quantity which is obtained by the implementation of . Then, by using Eqs.( 2), ( 3), (4), and ( 5) we obtain the quantum-errormitigated expectation value of Ô given by The idea of our QEM scheme is clearly represented in the second line of Eq. ( 8).The first term is the ideal expectation value while the second term represents the conduction of our QEM scheme.It is described as the subtraction between (the quantum computational error occurring on a real device) and ⟨ Ô⟩ ∆ AD 1 ρ d•••1 (theoretically-evaluated quantum computational error), which is computed by the ADeffect circuit group.The heart of the idea for doing this is that we have considered that the real noise effect δ AD . When the second term in the second line of Eq. ( 8) becomes small enough, we consider that we have accomplished in mitigating the error of quantum computation on a real device.Note that the quantum noise effects coming from δ AD (∆ AD ρ d•••1 ) are suppressed by multiplying ⟨ Ô⟩ (∆ AD ρ d•••1 ) real by τ (the second term in the right-hand side of the first line of Eq. ( 8)).This is because the lowestorder error AD effect on the implementation of Consequently, in the first-order perturbation theory with respect to τ we have established our QEM scheme by the usage of the AD-effect circuit group and is expressed by the formula given by Eq. (8).The above argument on QEM-scheme derivation can be straightforwardly generalized to register-bit systems for N q ≥ 2. In this case, 3) is represented as We can apply our QEM scheme to the N q register-bit system in the following way.We prepare N q register bits and one ancilla bit Q r0 , Q r1 , . . ., Q rNq−1 , Q a .Then we create an ensemble of quantum circuits composed of the j-th register bit Q rj (j = 0, 1, . . ., N q − 1) and the ancilla bit Q a which describes that Q rj is subject to the AD effect induced by the ancilla bit Q a .Namely, we create the ensemble of four types of quantum circuits composed of Q rj and Q a , the original quantum circuits given by U QC , the quantum circuits with additional Z-gate operations, and the AD-effect circuits A and B. By summing up all these quantum circuits, we obtain the AD-effect circuit group which enables us to perform QEM for N q -register-bit system under the AD effect.The total number of quantum circuits which compose the AD-effect circuit group is 3dN q + 1, and thus, it scales polynomial in dN q , which is not so high-cost computational performance.Before ending this section, let us explain two ways to extend our QEM formalism.Firstly, we can extend into cases of other quantum noise channels including generalized amplitude damping (GAD), phase damping (PD), their composite channels, and stochastic Pauli noise models such as bit flip, phase flip, bit-phase flip, and depolarizing channel.Secondly, we can create quantum-noise-effect circuit groups which en-able us to perform QEM for higher-order quantum noise effects.We present arguments on these two cases in Sec.I in the Supplementary Material.

IV. NUMERICAL SIMULATIONS
In this section, we numerically demonstrate our QEM scheme for four types of algorithms.For the quantum noise effect we choose AD effect.In Sec.IV A, as a preliminary of our QEM demonstration, we present the results of two algorithms: the algorithm composed of the initial X-gate operation and the repetition of the Hadamard gate acting on a single register bit and that composed of the initial X ⊗ X operation and the controlled-Hadamard gate acting on two register bits.In Sec.IV B, we show the results of QEM for a long-term quantum algorithm and here we choose quantum amplitude amplification (QAA) algorithm (Grover's search algorithm).In Sec.IV C. we show the results of recently developed NISQ quantum algorithms, Quantum Approximate Optimization Algorithm (QAOA).
In the following, let us explain the formalism of our noisy quantum simulations (numerical simulations of running quantum algorithms with real quantum devices performed by classical computers).As shown in Fig. 2, every time we apply an unitary (gate) operation, the N q register bits experience AD effects.Suppose that at a time t 0 the quantum states of the register bits were given by the density matrix ρ(t 0 ).According to the quantum master equation (1), when the unitary gate U has been applied to the register bits within the time interval ∆t the quantum state of the register bits at t = t 0 + ∆t is expressed by the density matrix where E AD [• • • ] is the superoperator which describes the AD effect on the N q register bits and it is given by [26,27,49,60,61,86] where n Qr0 , . . ., n Q rNq −1 = 0, 1, and cos 2 ϑτ 2 = e −τ .The operators M AD 0 and M AD 1 are the Kraus operators acting on single-qubit states and describe the influence of the AD effect on a single register bit during the time interval ∆t.Here we assume that the N q register bits homogeneously experience the single-qubit-state AD effect as described in Eq. (11).
For later convenience, let us introduce the notation which describes the operation of the unitary operator U on the quantum state ρ by T [ρ, U ](= U ρU † ).Let us write the quantum state generated by the unitary transformation U QC under the AD effect by ρ AD d•••1 .By using the superoperator Eqs. (11) and (12) are the basic equations of our noisy quantum simulations.Namely, we perform our noisy quantum simulations by identifying the quantum state ρ AD d•••1 in the above equation with ρ real d•••1 , which is the AD-affected quantum state generated by the unitary transformation U QC on real devices.We conduct QEM described by Eq. ( 8) for various values of τ by tuning the value of ϑ τ .When we execute the AD circutis A and B, we run quantum simulations of the N q + 1 qubit systems.The right-hand side of Eq. (11) becomes Note that the ancilla bit Q a is treated similarly as the other N q register bits such that Q a is subject to the same quantum noise effect described by the Kraus operators in Eq. (11) as the other N q register bits do.Correspondingly, we mitigate the quantum noise effects influencing both Q a and the N q register bits via Eq.( 8).
For performing our noisy quantum simulations, we have created two types of numerical codes.The first one is our original numerical code and the second one is the numerical code created by Qiskit [86].In our numerical code, we simply compute the matrix products of the density matrices, the unitary operations U k , the Kraus operators, and the additional operations such as Z, σ± , P 1 .Furthermore, we compute the trace operations between the density matrices and the physical operators O. Namely, our original numerical code is a density-matrix simulator.On the other hand, the Qiskit code is programmed by the two numbers, N QC and N samp .To understand them concretely, first let us discuss an example of quantum computing of a single-qubit system.We execute the given quantum circuit and we obtain an output state which is either |0⟩ or |1⟩.We repeat this process N QC times and say that we obtained |0⟩ for N |0⟩ times as the output state while we obtained |1⟩ for N |1⟩ times.Then, the probability weight of |0⟩ is NQC .Namely, the number N QC is the repetition number of quantum computing executed by the given quantum circuit, and the numerical simulations in our original code describes the quantum simulations in the limit N QC → ∞.Next let us explain what N samp is.It is the repetition number of "the execution of the quantum computation for N QC times under the given (fixed) quantum circuit".By introducing such a number, we effectively perform the quantum computation under the given quantum circuit with the repetition number N QC × N samp in our Qiskit code.In contrast to N QC , the repetition number N samp is not coming from foundations of quantum mechanics and we have introduced this quantity owing to the following two reasons.First, due to our survey there is an upper limit on N QC in the Qiskit program and is 10 7 .By introducing N samp , we become able to effectively perform quantum computations with repetition numbers greater than the upper limit of N QC (= 10 7 ) in the Qiskit program.Second, we consider that it is not enough to just show a single data point of "the quantum computational result obtained by the fixed quantum circuit and N QC " to see how largely it deviates from the true quantum computational result (quantum computation in the limit N QC → ∞).To present how largely the quantum computational results for fixed and finite N QC deviate (finite-size effects of N QC ) from the true ones, we introduce N samp and show visually such deviations.In order to describe the second reason more mathematically, let write a binary which labels a quantum state of qubit Q α (α = 0, . . ., N q − 1) by n Qα : |n Qα ⟩ with n Qα = 0, 1 and an output state is described by the computational basis states Let us say that we focus on the specific quantum state |n Q Nq −1 • • • nQ0 ⟩ and consider how many times it is obtained as the output state for given N QC .By saying that |n Q Nq −1 • • • nQ0 ⟩ has been obtained as the output state for

NQC
. Furthermore, we write the probability such that |n We redescribe such a circumstance as the binomial distribution denoted by B(p nQ Nq −1 •••n Q 0 , N QC ) and introduce the random variable X i such that in the ith round of quantum computation we take X i = 1 when the output state is In other words, the binomial distribution B(p, N QC ) approaches to the Gaussian distribution function given by the mean p nQ Nq −1 •••n Q 0 and the standard deviation Let us say that we perform quantum computing with accuracy ϵ.From the standard deviation of we can estimate the lower bound of N QC for doing this and is equal to . Note that in order to perform QEM with the accuracy ϵ the total repetition number of quantum computing gets larger with respect to d and N q due to the quantum-noise-effect circuit group and this issue is addressed later.The number N QC describes the repetition number of quantum computation owing to the given quantum circuit whereas the number N samp describes how many times you conduct the sampling for the expectation values of physical operators obtained by the N QC -repeated quantum computation.Owing to this sampling, the repetitive number of the quantum computations effectively becomes N QC × N samp , and our simulation results become more trustable.In our simulations we take N QC = 2 10 and N samp = 100 for both the original quantum circuit and each elementary circuit of quantum-noise-effect circuit group.Our original numerical code, on the other hand, is the code for a noisy quantum simulation in the limit N QC → ∞, and basically, it performs pure linear algebraical computations such as matrix-product and trace operations.We note that when we create the numerical codes with Qiskit, we need to be careful with how controlled-unitary operators are implemented.We have examined that on Qiskit program the controlled-unitary operators are implemented as the decomposition of U CX gates and single-qubit unitary gates.For example, the control-R y gate Therefore, when we simulate QAA for three-qubit systems and QAOA with our original code we implement U CRy(ϑ) [Q r0 ; Q r1 ] in the same way as Qiskit program does.
In order to quantitatively describe the validity of our QEM scheme we introduce the measure defined by The numerator of RT QEM in Eq. ( 13) describes the absolute of the difference between the expectation value owing to the noisy quantum simulation On the other hand, the denominator represents the absolute of the difference between the expectation value obtained by our QEM scheme 8)) and the ideal expectation value.In other words, the measure RT QEM in Eq. ( 13) is the ratio between the absolute of the error without QEM and the one with QEM.Thus, when RT QEM > 1 the expectation value , which means that our QEM scheme is working.In addition to the ratio RT QEM in Eq. ( 13), we display the results of the expectation values obtained by the ideal simulations, the noisy simulations without QEM, and the noisy simulations with QEM, and show explicitly the validity of our QEM scheme.Note that for our original code we take the noise-strength parameter to be ϑ τ = i o × 0.01 with i o = 0, 1, . . ., 50 while for Qiskit code we take ϑ τ = i Q ×0.05 with i Q = 0, 1, . . ., 10.For computing the expectation values we include the case i o = 0 and i Q = 0 whereas for computing the ratio RT QEM we omit i o = 0 and i Q = 0.This is because in this case we have and we encounter in the indefinite 0 0 .In the following, we create a subsection for each quantum algorithm and discuss the results in detail.

A. Preliminary
As a preliminary, let us conduct the noisy quantum simulations for two simple algorithms.The first algorithm is given by the unitary operation U QC pre1 = H ⊗d−1 • X, where H denotes the Hadamard gate.The second one is given by is the controlled-Hadamard gate composed of the control bit Q r0 and the target bit Q r1 .Here we take the circuit depth d to be d = 1 + 2 n with n being positive integers.Let us show the quantum circuits for these two algorithms in Fig. 5.By changing the values of d and τ , we numerically analyze how well our QEM scheme works in terms of these two parameters.Let us discuss from the results of noisy quantum simulations conducted by the quantum circuit in Fig. 5(a) and show them in Fig. 6.We have taken the physical operators Ô as Ô = X, Z.The horizontal axis represents ϑ τ , which can be regarded as the strength of AD effect.On the other hand, the vertical axis in Figs.
not as the expectation values of X with respect to the quantum states generated by U QC pre1 but as those of Z with respect to the quantum states generated by U QC pre1 and the subsequent operation of H, i.e., we have switched the basis vectors for the measurement from the computational basis to {|+⟩, |−⟩}, where |+⟩ = H|0⟩, |−⟩ = H|1⟩.Let us analyze our simulation results by comparing the behaviors of the expectation values and the ratio RT QEM in Eq. ( 13) as functions of ϑ τ .In this way, we can clearly see whether our QEM scheme is working or not, and for this purpose in the following we rewrite RT QEM as RT QEM (ϑ τ ) to emphasize that they are the functions of ϑ τ .Furthermore, we introduce the angle ϑ c τ such that RT QEM (ϑ τ = ϑ c τ ) = 1, which indicates that the point ϑ τ = ϑ c τ is the critical point of our QEM to become failed.Let us look from the results shown in Figs.6(a) and 6(c) by focusing on how the behaviors of , and RT QEM change by increasing ϑ τ .In Fig. 6(a), as the def- , and correspondingly, in Fig. 6(c) we see RT QEM (ϑ τ ) ≥ 1.As we increase the value of ϑ τ from ϑ c τ , the absolute and correspondingly, the ratio RT QEM (ϑ τ ) decreases monotonically from one.For the region ϑ τ ≥ ϑ c τ to improve the quality of our QEM we need take into account higher-order AD effects and establish QEM schemes for mitigating them and we expect the value of ϑ c τ to become larger.Next, let us analyze how the quality of our QEM becomes when we vary the circuit depth d.We see that for every | become bigger and RT QEM (ϑ τ ) decreases as we increase d.This is reasonable because when d gets larger the amount of error gets bigger.For the results in Figs.6(b) and 6(d), basically we see that both the expectation values of Z and RT QEM (ϑ τ ) show the similar behaviors as those for Ô = X: (I) the validity of QEM (RT QEM (ϑ τ ) ≥ 1) in the range 0 < ϑ τ ≤ ϑ c τ and monotonic decrease of RT QEM (ϑ τ ) for ϑ τ > ϑ c τ , and (II) worsening of the quality of our QEM for large d.In contrast to the above characteristics of ⟨X⟩ and ⟨Z⟩, we have numerically verified that the expectation value of Y takes zero for any ϑ τ .This is because when the density matrix ρ is a real matrix the expectation value ⟨Y ⟩ is zero.Since both the quantum algorithm and the AD effect are described by real numbers (see also Eq. ( 7) or the Kraus operators in Eq. ( 11)) the density matrix generated by these two things is real and we have ⟨Y ⟩ = 0.
Let us discuss the results in Fig. 7.They are the noisy simulation results of the quantum algorithm given by U QC pre2 (see the quantum circuit in Fig. 5(b)) and we have taken d = 1 + 2 n as in the case of simulations for U QC pre1 .Here we have simulated RT QEM (ϑ τ ) for the expectation values of the operators Ô = ZX, ZZ.As similar to the computa- , we have computed not as the expectation values of ZX with respect to the quantum states generated by U QC pre2 but as those of ZZ with respect to the quantum states generated by U QC pre2 and the subsequent operation of H on Q r1 .Overall, we see the same characteristics with the cases of Ô = X, Z: the characteristics (I) and (II) mentioned above.For any ϑ τ , the ratio RT QEM (ϑ τ ) for the noisy simulations of U QC pre2 are smaller than those of U QC pre1 .This is because U QC pre1 is solely comprised of the single-qubit gates (X and H) while U QC pre2 is constructed by n-operation of the controlled-Hadamard gate (two-qubit gate), and thus, the bigger amount of errors are accumulated in the latter case.The difference between the characteristics of noisy simulations for U QC pre1 and those for U QC pre2 , although it is not an essential point for the validity of our QEM, is that we see both one minima and one maxima in each plot for RT QEM (ϑ τ ) in Fig. 7(c) while only one maxima appears in Fig. 7(d).Let us denote the point where RT QEM (ϑ τ ) takes the minimum (maximum) by ϑ min τ (ϑ max τ ): note that these values depend on d.We can understand why these points emerge by looking at Figs. 7(a) and 7(b).Let us explain from the origins of the minima and the maxima in Fig. 7(c) by looking at the plots in Fig. 7(a).In the range 0 < ϑ τ ≤ ϑ min As a result, the minima appears at ϑ τ = ϑ min τ whereas the maxima emerges at ϑ τ = ϑ max τ in Fig. 7(c).The origin of the maxima in Fig. 7(d) can be similarly explained by looking at the plots in Fig. 7(b).In the range 0 and consequently, the maxima appears at ϑ τ = ϑ max τ .Like the case of the noisy simulations of U QC pre1 , the density matrices are generated as real matrices (the unitary transformation U QC pre2 as well as the AD effects are described by real numbers), and the expectation values of the Pauli operators, IY, XY, Y I, Y X, Y Z, ZY are zero for both ideal and noisy simulations.Here we have rewritten 1 2×2 as I for convenience.Note that the expectation value of the identity operator (= 1 4×4 : four by four identity operator) is one for any quantum state including noise-affected quantum states since the trace of density matrix is one for any quantum state.In other words, it is unnecessary to do QEM for the expectation value of the identity operator.Note, however, that when leakage occurs the trace preservation is not held anymore and we need to consider QEM for the error induced by the leakage.

B. Quantum Amplitude Amplification
By taking account of the previous analysis, let us apply our QEM scheme to Quantum Amplitude Amplification (QAA) [33] for three-qubit systems: two-register bits and one oracle bit.One of the important application of QAA is the database retrieval and the quantum algorithms for this is called the Grover's search algorithm [33,[92][93][94][95]. Let us denote the (classical) oracle function by f and a binary by x which takes "00", "01", "10", "11".We consider that we have only one solution of f and write it by x * , which satisfy f (x * ) = 1, and assume x * = "11": for x = "00", "01", "10" we have f (x) = 0.The oracle operator O QAA can be implemented on a quantum circuit by using one oracle bit 8. Quantum circuits for U QAA in Eq. ( 14).(−1) , where the superposition state is created by applying H • X on the oracle-bit state |0⟩ Qo .In our case, (−1) f (x) = −1 when |x⟩ Qr0Qr1 = |11⟩ Qr0Qr1 , and the oracle operator O QAA is equivalent to the Toffoli gate comprised of the two controlled bits Q r0 and Q r1 and the target bit Q o [33], and write it by To construct QAA we need one more unitary transformation and that is QAA is given by the unitary operation [33] We show the quantum circuit for the unitary operation U QAA in Fig. 8 [33]: the quantum circuit for the oracle operator FIG.11.Histogram of the probability distribution of the computational basis states for QAA simulations given by U QAA in Eq. ( 14).Here we have set ϑτ = 0.2.
II in the Supplementary Material.We note that on the quantum circuit in Fig. 8, what is actually implemented is −U G = U ψ • O QAA and the global phase factor (−1) does not affect our result.The unitary operation U G is called the Grover operator and k is the repetitive number number of its operation and here we have k = 1.After the operation of U QAA in Eq. ( 14), ideally both of the probability weights of |110⟩ and |111⟩ are 1  2 , and thus, the probability of obtaining the quantum state |11⟩ as the output state is 1  2 × 2 = 1, which implies the success of searching the solution x * = "11".By taking account of the above theoretical framework, we examine whether our QEM scheme works or not for QAA given by U QAA in Eq. ( 14) by computing the probability weights of |110⟩ and |111⟩ which we name as P 110 and P 111 , respectively, and show these results in Fig. 9. Solid lines in Figs.9(a by P Qiskit noisy and P Qiskit QEM , respectively.Similarly, in Figs.9(c) and 9(d), we have denoted the ratio RT QEM (ϑ τ ) calculated by our Qiskit code by RT Qiskit QEM : for the results obtained by our original code we have just used the notation RT QEM for describing them.Let us look from the simulation results of P 110 and the associated ratio RT QEM (ϑ τ ) given by Figs.9(a) and 9(c), respectively.In the range 0 ≤ ϑ τ ≤ 0.5, overall the simulation results with QEM are numerically closer to the ideal values than the noisy simulation results without QEM, and correspondingly, we have RT QEM (ϑ τ ) > 1.The similar characteristics can be seen in Figs.9(b) (simulation results of P 111 ) and 9(d) (RT QEM (ϑ τ ) for P 111 ).For the results obtained by our Qiskit code, the deviation between RT QEM and RT Qiskit  QEM becomes prominent in the small ϑ τ region and we consider this as follows.When noise strength ϑ τ is weak enough, on the Qiskit code the difference between the noisy value and the ideal value is very tiny such that our QEM becomes invalid and RT Qiskit  QEM gets lower than one.On the other hand, we see that some red points are above RT QEM (ϑ τ ) = 1.We consider that by greatly increasing N QC , we expect that RT Qiskit  QEM approaches to RT QEM .As a result, our QEM works for the noisy simulations of both P 110 and P 111 .To show clearly that the simulation results obtained by our Qiskit code approach to those obtained by our original code by increasing N QC , in Figs.10(a) and 10(b) we show the N QC dependencies of P 110 and P 111 , respectively.In these plots the horizontal axes represent N QC whereas the vertical axis in Fig. 10(a) represents the numerical values of P 110 and that in Fig. 10(b In addition to these two variances, we also plot the quantity defined by . Such a quantity describes the deviations of ⟨O⟩ Qiskit QEM (i, N QC ) from ⟨O⟩ QEM owing to the finite effect of N QC and plays the role of variance.Like we take it as the linear function of N QC given by the coefficient α Qiskit QEM (O, N samp , N q , d).Note that the dependency of the coefficient α Qiskit QEM (O, N samp , N q , d) in terms of N q and d originates in the size of the quantum-noiseeffect circuit group 3N q d.Owing to our simulation, we obtain α Qiskit QEM (O, N samp , N q , d) = 3. , however, is about 1.34 for both cases which implies that the broadening of the deviations is not so big.We leave the detailed mathematical analysis on as well as the coefficient α QEM (O, N samp , N q , d) as our future work.In addition to N samp = 100, we also perform the simulations for As a result, by increasing N samp both ⟨O⟩ Qiskit noisy (i, N QC ) and ⟨O⟩ Qiskit QEM (i, N QC ) approach to ⟨O⟩ noisy and ⟨O⟩ QEM , respectively, owing to the law of large numbers.Let us also show the simulation results of the rest of the six probabilities of the computational basis states for ϑ τ = 0.2 as the histogram in Fig. 11, which also includes P 110 and P 111 .The ideal values of these six probabilities are all zero and we see that the noisy simulation results with QEM are numerically close to them compared with those without QEM, which indicates that our QEM scheme also works for the other six probabilities.
In addition to the above simulation, let us present the simplified version of QAA for the two-qubit systems [96].In this problem setting, we consider the effective two-dimensional space spanned by the two Bell states , and write their superposition state by Our goal is to amplify the probability amplitude c Φ + .Here we take the initialization operator to be while we take the Grover operator to be k is the number of V G to be applied.In total, the unitary operation for running QAA is given by V QAA = (V G ) k • V init .We present the corresponding quantum circuit in Fig. 12.As similar to the above case, on the quantum circuit we implement Ṽs instead of V s .In the following, let us take a look at the meanings of the three unitary operations V init , V s , and V ω .First, the unitary operation V init generates the superposition of |Φ + ⟩ and Second, the unitary operation V ω is the oracle operator and when it is applied to the initial state |s⟩ we have , the oracle operator V ω is the operator such that it reverses the sign of the Bell state |Φ + ⟩.Third, V s is the reflection of the vector V ω |s⟩ with respect to the vector |s⟩.When we operate V G on |s⟩ for k times we have and since ϑ V = π/3 we have k = 1.We can understand the geometrical meaning of the operation V G as follows.Let us consider the effective three-dimensional space spanned by the two Bell states |Φ + ⟩ and |Ψ + ⟩ and the vector |ξ⟩ which is perpendicular to both |Φ + ⟩ and |Ψ + ⟩.Furthermore, we call the axis which is parallel to |ξ⟩ (perpendicular to the twodimensional plane spanned by |Φ + ⟩ and |Ψ + ⟩) as ξ-axis.The unitary operation V G is the rotation about ξ-axis by the angle ϑ V in this effective three-dimensional space.We can verify whether |Φ + ⟩ has been generated as the output state or not by measuring the probabilities of |00⟩ state and |11⟩ state, which are denoted by p 00 and p 11 , respectively.In other words, the probabilities p 00 and p 11 are the expectation values of the projection operators P 00 and P 11 , respectively, where P 00 (P 11 ) is the projection operator of |00⟩ (|11⟩) state.Namely, we run noisy quantum simulation for Ô = P 00 , P 11 and perform QEM on them: Note that in the case of the ideal simulation we obtain p 00 = p 11 = 1 2 .We plot the numerical results of p 11 and p 00 for the range 0 ≤ ϑ τ ≤ 0.5 in Figs.13(a) and 13(b), respectively, and in Figs.13(c) and 13(d) we plot RT QEM (ϑ τ ) for p 11 and p 00 , respectively.All these results shown here are obtained by our original numerical code and we have denoted ⟨P 00 (11) by P ideal , P real , and P QEM , respectively.Let us look from the results of the probability p 11 .In Fig. 13(a) we see that for any ϑ τ the absolute of the deviation and correspondingly, as we see in Fig. 13(c) the ratio RT QEM (ϑ τ ) is greater than one.Therefore, our QEM scheme works well for the noisy simulation of p 11 .In contrast, in Figs.13(b) and 13(d) we see that the probability p 00 shows essentially a different behavior.That is the expectation value without QEM ⟨P 00 , and correspondingly, the ratio RT QEM (ϑ τ ) is lower than one.Such a characteristic is understood as follows.Firstly, we have analytically examined that the expectation value ⟨P 00 ⟩ ρ real d•••1 does not include first-order term in τ , i.e., ⟨P 00 ⟩ ∆ AD 1 ρ d•••1 = 0.The lowest-order term included in the numerator of RT QEM (ϑ τ ) is O(τ 2 ).Secondly, due to our QEM the lowest order of the denominator of RT QEM (ϑ τ ) is also O(τ 2 ).As a result, the ratio RT QEM (ϑ τ ) becomes lower than one, which indicates that it is not appropriate to adopt our QEM scheme.We consider that this is because our QEM scheme described by Eq. ( 8) is the scheme for mitigating the first-order AD effect.In the limit of ϑ τ → 0, the ratio RT QEM (ϑ τ ) takes finite value, and analytically it is the ratio between the absolute of the coefficient of . We note that we have performed two types of simulations with our original code.In the first one we have directly implemented U CZ [Q r1 ; Q r0 ] gate while in the second one we have implemented it via the decomposition According to the results of these two simulations, we have analyzed that on the Qiskit code U CZ [Q r1 ; Q r0 ] gate is automatically implemented by the above decomposition since the result of the second simulation with our original code shows better matching with that obtained by our Qiskit code.In such a case, the first-order term in τ appears for ⟨P 00 ⟩ ρ real d•••1 and our QEM works well.Besides p 00 and p 11 , let us briefly discuss the noisy simulation results of the probability weights of |01⟩ and |10⟩ and write them by p 01 and p 10 , respectively.We show them in the his-togram in Fig. 14 which describes the probability distribution of the computational basis states of the two register qubits Q r0 and Q r1 .Here we have taken ϑ τ = 0.2.We see that like the probability weight of |11⟩, our QEM scheme works for the probability weights of |01⟩ and |10⟩.
Let us end this subsection by giving the following comment.In the previous subsection, we have seen that our QEM scheme becomes meaningless in the cases when the expectation values of the ideal simulations are equivalent to those of noisy simulations such as the simulation for the expectation value ⟨Y ⟩.Besides these cases, our QEM scheme represented by the formula in Eq. ( 8) does not work when noisy expectation values do not include the first-order term in τ like the noisy simulation for the probability p 00 discussed above.In other words, if we construct the QEM formula which describes the mitigation for a higher-order quantum noise effect, which is discussed in Sec.IA in the Supplementary Material, by using it we become able to accomplish the noisy quantum simulation obtaining RT QEM (ϑ τ ) > 1.In practice, however, when we run quantum algorithms on real quantum devices we cannot compute RT QEM (ϑ τ ) since we cannot compute ideal expectation values.We can check whether noisy expectation values include the first-order terms in τ or not, for example, in the following way.We perform two types of QEM, QEM of both the first-and second order quantum noise effects and the one of only the second-order effect.Let us denote the density matrices obtained by the former QEM and the latter one by ρ QEM 2nd and ρ QEM 2nd,only , respectively.Next, we introduce the measure where O is a physical operator.If the first-order QEM fails (noisy expectation values do not include the first-order terms in τ ) then the measure M 1st/2nd is O(τ 2 ).On the other hand, if the firstorder QEM succeeds (noisy expectation values include the first-order terms in τ and is mitigated) then M 1st/2nd is O (1).By extending this approach we can examine whether noisy expectation values include higher-order quantum noise effects or not.We consider, however, that the failure of the first-order QEM for the noisy quantum computation when the linear order in τ does not appear is not so crucial compared with the case when we have failed in mitigating the linear-order quantum noise effects included in noisy expectation values, and indeed we can see this by looking at Fig. 14.For the probability weight of |00> state, the numerical differences among the three expectation values, , are small.On the other hand, for the probability weight of |11> state, As a final example, let us apply our QEM scheme to the noisy simulation of the variational quantum algorithm called Quantum Approximate Optimization Algorithm (QAOA) [27][28][29][30][31][32][33].In the following, we analyze QAOA for a max-cut problem which is to divide vertices (nodes) of a given graph into two groups so that the number of edges connecting two vertices belonging to the different groups is maximized and is a NP (Non-deterministic Polynomial time)-hard problem.
First, we discuss from a theoretical framework of a classical approximate optimization.We express the given graph G as G = (V, E), where V = {v 0 , . . ., v i , . . ., v NV−1 } is the set of vertices with N V denoting their total number and for each vertex v i the binary value is the set of the edges with ⟨v i v i+1 ⟩ denoting the edge connected by the vertices v i and v i+1 .The quantity C ij (i, j = 0, . . ., N V − 1 with i ̸ = j) is the adjacency matrix element (weight) for the edge ⟨v i v j ⟩ which is semi-positive.Let us write the N V strings of z i by z = (z 0 , . . ., z NV−1 ).The goal of a classical approximate optimization is to minimize the cost function or equivalently to maximize the ratio r CAO (≤ 1) which satisfies where C min is the minimum value of C(z).In our simulation, as illustrated in Fig. 15 we adopt the square graph given by the four vertices v 0 , v 1 , v 2 , and v 3 , and the edges are ⟨v 0 v 1 ⟩, ⟨v 1 v 2 ⟩, ⟨v 2 v 3 ⟩, and ⟨v 3 v 0 ⟩, and take C ij = 1 for any edge ⟨v i v j ⟩.Next, we discuss the theoretical framework for QAOA.The four vertices v 0 , v 1 , v 2 , and v 3 are encoded in four qubits Q r0 , Q r1 , Q r2 , and Q r3 , respectively, and the values z i z j in the expectation values of the operators Z Qi ⊗ Z Qj .The cost function C(z) in Eq. ( 15) is given by the expectation FIG. 15.
of the Hamiltonian where the symbol (i, j) (i, j = 0, 1, 2, 3) denotes the summation for the edges connected by the qubits Q ri and Q rj under the square-graph structure in Fig. 15.In this simulation the physical operator Ô is the Hamiltonian H C in Eq. ( 17).The unitary operation for running QAOA, which we denote by U QAOA , consists of three elements.The first one is the unitary operation for creating the reference state and is given by the Hadamard-gate operation on all four qubits, U int = H ⊗4 .The other two are the unitary operations U C (ϑ QAOA j ) and U X (φ QAOA j ) which are generated by the Hamiltonian H C in Eq. ( 17) with the angle ϑ QAOA j and the term H X = i X Qri with the angle φ QAOA j , respectively: in QAOA H X is called the transverse-field (mixing or driving) term.The two types of angles ϑ QAOA j and φ QAOA j (j = 1, . . ., p) are the variational parameters and p is the repetition number of applying the unitary operation ) (the number of iteration), which determines the accuracy of QAOA.In total, U QAOA is given by The quantum circuit for U QAOA in Eq. ( 18) is presented in Fig. 16.In our simulation we set p = 2 and the circuit depth is d = 15.The unitary operation U X (φ QAOA j ) is implemented by the R x gate (rotation about x axis) with the angle 2φ QAOA j .Meanwhile, the quantum circuit for U C (ϑ QAOA j ) is composed of the sets of the quantum gates where R z denotes the rotational gate about z axis and the associated angle is 2ϑ QAOA j .Corresponding to Eq. ( 16), the goal of QAOA simulation is to compute and minimize the expectation value   ).Let us also call the expectation value C ϑ QAOA , φ QAOA in Eq. ( 19) as the cost function and its minimization is equivalent to the optimization of the variational parameters ϑ QAOA j , φ QAOA j , and write them by ϑ QAOA,opt j , φ QAOA,opt j .We show their values in Sec.III in the Supplementary Material.Let us now discuss our simulation results shown in Figs.17-19.All these results are obtained by computing the probability distribution of the computational basis states since the Hamiltonian H C in Eq. ( 17) is given by the Z gate operations.First, let us take a look at the simulation results in Fig. 17.In Fig. 17 ).The dashed black line, the blue solid line, and the orange solid line are C ideal , C noisy , and C QEM , respectively, and they are all obtained by our original code.On the other hand, the blue and orange circles are , respectively, and they have been calculated by our Qiskit code.We have denoted by C Qiskit noisy and C Qiskit QEM , respectively.We have plotted 100 circles for each ϑ τ .We see that in the range 0.1 ≤ ϑ τ ≤ 0.5 our QEM works well.In Fig. 17(b), we have plotted the ratio RT QEM (ϑ τ ).The black curve is the result obtained by our original code while the red circles are those obtained by our Qiskit code and 100 circles are plotted for each ϑ τ .Corresponding to the result shown in Fig. 17(a), the ratio satisfies RT QEM (ϑ τ ) > 1.Like the results in Fig. 10, we present the N QC dependencies of the probabilities P 0101 and P 1010 in Figs.18(a) and 18(b), respectively.We take N QC = 2 nQC with n QC = 8, 9, . . ., 14, N samp = 100, and ϑ τ = 0.2.We see that as we increase  QEM (i, N QC ) from ⟨O⟩ QEM due to our QEM method is not so large (not so high computational cost).
We have also performed our simulations for N samp = 1000 and we obtain Finally, let us explain the results in Fig. 19.Here we have presented the histogram of the probability distribution of the computational basis states for the ideal case (P ideal ), the noisy case (P noisy , P Qiskit noisy ), and the case with QEM being performed (P QEM , P Qiskit QEM ), with setting ϑ τ = 0.2.We see that ideally the probabilities of the two quantum states |0101⟩ and |1010⟩ are both equal to 0.5.This implies that under the optimized variational parameters the cost function C ϑ QAOA , φ QAOA becomes minimized such that the four qubits Q ri are partitioned into the two groups Q r0 , Q r2 and Q r1 , Q r3 , and it implies that all the edges of the square are to be cut, i.e., the maximum number of edges to be cut is four.Correspondingly, as shown in Fig. 17

V. QEM SCHEME IMPLEMENTATION
In this section, we demonstrate our QEM scheme using two IBM Q Experience processors, ibmq_belem and ibm_perth [87].In Fig. 4 in Sec.V of the Supplementary Material, we show schematics of spatial configurations for qubits in these machines: Fig. 4(a) is the illustration for ibmq_belem whereas Fig. 4(b) is that for ibm_perth.As similar to the quantum simulations demonstrated in Sec.IV, we examine the efficacy of our QEM scheme for the real quantum devices by varying the depths of the quantum algorithms to be run.We do this for two quantum algorithms, U QC pre1 in Fig. 5(a) and a quantum algorithm for a two-qubit system defined by , where n rep is an integer.For the usage of ibmq_belem, we choose qubits Q 1 and Q 3 as register qubits whereas we use Q 2 (Q 4 ) as an ancilla bit for mitigating the quantum noise effect on Q r1 (Q r3 ).On the other hand, we use Q 1 and Q 3 as register qubits while we use an ancilla bit Q 2 (Q 5 ) for QEM of the quantum noise effect on Q 1 (Q 3 ) for the usage of ibm_perth.In the following, we rewrite Q 1 and Q 3 as Q r1 and Q r3 , respectively, to emphasize that they are used as register bits.Our QEM scheme can be applied provided that noise parameters (strengths) are given a priori.Thus, we start from the acquisition of the quantum noise strengths of these devices (noise characterization).Then, we discuss the results of our QEM method obtained with these machines (implementation).

A. Noise Characterization
Let us first discuss from how to obtain the AD strength τ (= γ∆t) or the T 1 time.The relation between the decay rate γ and the T 1 time is T 1 = 1 γ(2n+1) , where n is the Bose-Einstein distribution function.For superconducting qubit systems, qubit frequencies and temperatures are about 5.0 GHz (see also Table II, Table IV, and Table V in the Supplementary Material) and 10 mK [9,17], respectively.As a result, the Bose-Einstein distribution function n is estimated to be n ∼ 10 −11 and we approximate n to be zero.The value of T 1 can be obtained as follows.First, we prepare a single qubit and apply the X gate so that the qubit is in the |1⟩ state (excited state).Next, we let the qubit relax until certain time t relax a (a = 1, 2, . . ., N m and 0 < t 1 < t 2 < • • • < t Nm ), measure the output state of the qubit which is either |0⟩ or |1⟩, and repeat this process to obtain the probability weights of the |0⟩ and |1⟩ states.The total exposure time of the qubit is t tot,T1 a = ∆t X + t relax a , where ∆t X denotes a X-gate operation time.This is equivalent to measuring the dynamics of the expectation value ⟨P 1 ⟩(t tot,T1 a ) (the dynamics of the probability such that |1⟩ is to be measured as the output state), and by doing an exponential fitting on the ⟨P 1 ⟩(t tot,T1 a ) plots we can extract the value of T 1 .We can also extract a T 2 time and here let us consider a T * 2 time which is obtained by the Ramsey measurement explained in the following and hereinafter we just denote T * 2 time as T 2 time.Basically, what we do is first we apply the Hadamard gate, let the qubit relax for t relax a , apply the second Hadamard gate, measure the output state, and repeat this process: The actual gate operations implemented on these experiments are not the two Hadamard gates owing to the issue of transmon qubit systems and its details is described in Sec.IV in the Supplementary Material.The above procedure is essentially equivalent to measuring the dynamics of ⟨P 1 ⟩(t tot,T2 a ), where t tot,T2 a = 2∆t H + t relax a with ∆t H denoting Hadamard-gate operation time.From the ⟨P 1 ⟩(t tot,T2 a ) curve we obtain a T 2 time.In the Supplementary Material IV, we explain how to estimate the T 1 and T 2 times with presenting the experimental data plots of the expectation values: Fig. 2 and Fig. 3 show the data plots of the T 1 and T 2 times, respectively, and Table I lists the values of them for ibmq_belem.
In Sec.V in the Supplementary Material, we present the physical properties (single-qubit and two-qubit gate times, qubit frequencies) of ibmq_belem and ibm_perth: Table II and Table III . By using the data in Sec.IV and Sec.V in the Supplementary Material we are now able to implement our QEM scheme.Before showing the results, let us comment five things about our experiment.Firstly, the values of T 1 and T 2 times differ from qubit to qubit on a real device due to imperfection (inhomogeneities of T 1 and T 2 ).To take these inhomogeneities into account and extract T 1 and T 2 times, we need to perform the above two procedures independently on each qubit.Secondly, a gate time ∆t is actually gate dependent as shown in Table II-Table VII in the Supplementary Material.By taking the inhomogeneities of both the T 1 and T 2 times and the gate times into account, instead of τ , which does not depend on both the qubits and the quantum gates, we perform our QEM scheme by using a perturbative parameter τ jk = γ j ∆t k , where j(= 0, . . ., N q − 1) is the index for qubit numbering whereas k(= 1, . . ., d) is that for labeling the gates.In the following, we denote the T 1 and T 2 times of the qubit Q j by T 1,j and T 2,j , respectively.Thirdly, T 1 and T 2 times fluctuate temporally on a real device.To run our QEM scheme, it is necessary to record the data of T 1 and T 2 times day-by-day and use these values.Thus, we indicate the time and the date in the coordinated universal time (UCT) when we list the data of T 1 and T 2 times.Fourthly, the CZ gate operation is decomposed into the form , and according to the transpilation of IBM Q Experience programming , and U QC imp,2 has been executed as Note that 1 denotes the identity operator.Fifthly, the time when the quantum computation has been conducted (see the captions of Fig. 20, Fig. 21, and Fig. 22) indicates the time when the original circuit and the quantum-noise-effect-circuit group have been executed.Finally, the measured (experimental) value of τ , which we denote by τ exp , can be differ from the true value of τ due to, for instance, imperfection of experimental apparatuses.Although that is the case, the conduction of our QEM method is said to be a success provided that the condition RT QEM > 1 is satisfied.

B. Implementation
As mentioned previously, in real quantum devices both T 1 and T 2 times and gate operation times are inhomogeneous, and moreover, both AD and PD effects exist.By taking these elements into account, we perform our QEM scheme by improving the formulas given by Eqs. ( 3) and ( 8) as

⟨ Ô⟩
We note that the quantum circuits such that the additional operations {Z, σ− , P 1 } are inserted after the virtual Z gates are not needed for (or do not contribute to) the calculation of Eq. ( 20) because the operation time of the virtual Z gate is zero.
First, let us discuss from the results for U QC pre1 which are shown in Fig. 20.Here all the horizontal axes represent d = n rep + 1, where the repetition number n rep is taken to be 16,32, and 64.Fig. 20(a) presents a quantum computation result of expectation values ⟨Z 1 ⟩ obtained by ibm_perth whereas the plots in Fig. 20(b) is a quantum simulation result of ⟨Z 1 ⟩.On the other hand, Fig. 20(c) displays a quantum computation result of RT QEM with ibm_perth and Fig. 20(d) plots a quantum simulation result of RT QEM .In Fig. 20(a), we observe that both the expectation values with and without QEM get farther from the ideal expectation value (black dashed line) rapidly as we increase n rep , and correspondingly, the monotonic decreasing of RT QEM is exhibited in Fig. 20(c), which are similar to the characteristics in Figs.6(a) and 6(c).Compared with these results, however, RT QEM decreases gradually since not only the noisy expectation values but also the expectation values with QEM increase rapidly.In contrast, in Fig. 20(b), we see that while the expectation values without QEM or noisy expectation values (blue open circles) increase gradually as d (or n rep ) gets larger the expectation values with QEM (orange open circles) take almost the same value.In Fig. 20(d), the monotonic increasing of RT QEM is exhibited.In the experiment of the implementation of U QC pre1 , both the expectation values and RT QEM obtained by the quantum computation and those by the quantum simulation show qualitatively different behaviors.On the other hand, all the ratios RT QEM are greater than one.As a result, our QEM scheme has worked, however, not effectively.
= 1 • 1: see also the fourth comment in Sec.V A. Note that the number of the implementation of the virtual Z gate is not included in the depths of the quantum algorithms.Let us explain from the results in Figs.21(a) and 21(b).We observe that both of these results show qualitatively the same behavior, i.e., while the noisy expectation values (blue open circles) gradually get farther from the ideal expectation value (black dashed line) as n rep increases, the expectation values with QEM (orange open circles) approximately take the same value.On the other hand, in Figs.21(c) and 21(d) the ratios RT QEM exhibit the monotonic increasing, and moreover the ratio RT QEM is greater than one for every d: RT QEM shown in Fig. 21(c) are larger than those in Fig. 20(c) for every d.In addition to the quantum computation with ibm_perth, we have also conducted a quantum computation using ibmq_belem and a quantum simulation for  I in the Supplementary Material.The time when the latter data has been taken (12:33) is closer to the time when the quantum computation has been conducted (11:43), and thus we consider that the latter RT QEM are greater than the former ones: see also the third comment in Sec.V A. As a result, we observe RT QEM > 1 for every quantum computation in our experiment, and we consider that our QEM works for the IBM Q Experience processors.
To summarize our experiment of the implementation of our QEM scheme, the quantum computations for U QC pre1 and U QC imp,2 show the different behaviors although the same machine has been used: the former case exhibits the different behavior from the simulation result while the latter case shows qualitatively the similar behavior.Moreover, the ratios RT QEM for U QC imp,2 are larger than those for U QC pre1 , which are in contrast to the results in Fig. 6 and Fig. 7: On the whole, RT QEM for U QC pre2 are basically smaller than those for U QC pre1 .One way to interpret the characteristics in Fig. 21(c) and Fig. 22(c), the increasing behavior of RT QEM with respect to n rep , is as follows.When the noise strength is too small the noisy expectation values become sufficiently close to the ideal ones and in such circumstances the conduction of our QEM scheme can give rise to negative effects since computers can treat up to certain digits.Indeed, such a characteristic has been observed in the quantum simulation results in Figs.9(c) and 9(d) and Fig. 17(b) in the range 0 < ϑ τ ≤ 0.1: for a superconducting qubit system T 1 ≈ 100 µsec and ∆t ≈ 100 nsec and ϑ τ is estimated to be ϑ τ ≈ 0.063.Thus, in order to utilize our QEM method effectively we need to use it under circumstances with moderate quantum noise strengths or for running moderately long quantum algorithms under weak quantum noise effects: such a way of interpretation, however, cannot be adapted to understand the characteristics in Fig. 20(c).Consequently, the interpretation of the discrepancy between the quantum simulation result and the quantum computational result for U QC pre1 and the discrepancy between the characteristic of RT QEM for U QC imp,2 and that for U QC pre1 remain unresolved for this experiment.
We note that our experiment has been conducted under restricted conditions such as the time for which we could have used the IBM Q Experience processors and the machines which have been available.Provided that we have no such restrictions, let us give several comments on how to improve our results.The first way is to increase the value of N samp .As indicated in Fig. 10 and Fig. 18, by increasing the value of N samp we expect that the expectation values with QEM approach to unique values and we become clearer to see whether our QEM scheme is working or not.The second way is to mitigate other types of errors.By combining our QEM scheme with QEM methods for other errors such as state preparation and measurement errors or errors according to other quantum noise channels such as crosstalk [97], we anticipate that the value of RT QEM becomes larger.
In addition to the above discussion, let us consider the effectiveness of the implementation of our QEM scheme on different quantum hardware and here we choose ion trap qubit systems.Ion trap qubit systems are engineered, for instance, as linear chains and a two-qubit gate operation can be exploited such that it can be performed on any pair of qubits [98]: In contrast, in order to implement a two-qubit gate on two separated qubits, say Q a and Q b , on the superconducting quantum devices which have been used in this experiment we need to insert SWAP operations acting on qubits which locate between Q a and Q b [87].Therefore, all quantum algorithms as well as quantum-noise-effect circuit groups are able to be implemented as indicated by the quantum circuits for them.In other words, we can harness our QEM scheme on ion trap qubit devices without reformulating the quantum circuits.Next, let us discuss quantum noise in ion trap qubit systems.The quantum noise occur in, for instance, hyperfine-state type ion-trap qubit systems are considered to be the phase damping and T 2 times are about 10 sec while single-qubit gate times are around 10 microseconds and two-qubit gate times are about 100 µsec [7,12].By setting T 2 = (2γ p ) −1 , we have τ p = γ p ∆t ≈ 5 × 10 −6 , where we have set ∆t = 100 µsec.τ p is sufficiently small and therefore, we expect that our QEM scheme also works for ion-trap NISQ devices.From these ingredients, we expect that the implementation our QEM scheme works more effectively and is more suitable for ion trap qubit systems compared with the superconducting qubit systems.

VI. COMPARISON WITH OTHER METHODS
We make comparisons between our method and other methods.Although many types of QEM methods have been proposed up to now [26,27,, here we focus on the following three methods, probabilistic error cancellation (PEC) [26,27,49,52,57,60,65,69,72,74,75,77], zero noise extrapolation (ZNE) [26, 27, 49-52, 73, 75, 77], and error suppression by dearangements (ESD) [67,76,77] and virtual distillation (VD) [66,77].In Table I, we summarize and present the comparison between our method and the others in terms of the number of ancilla bits N a and the number of additional circuits N AQC .We choose PEC from these three methods and numerically compare with our method and the reason we do this is the following.Both methods are theoretically similar such that they evaluate quantum noise effects on quantum states by quantum computational operations and perform QEM which are represented as sums of expectation values yielded by ensembles of quantum circuits including original circuits and quantum circuits containing additional operations.
Comparison between our QEM method in kth-order perturbation theory and the other three methods.Here b is a positive constant and ϵerr is a gate error rate.The number of ancilla bits Na is equal to (n − 1)Nq for VD without ancilla bits.

A. Comparison with PEC
First, we make a comparison between our method and PEC [26,27,49,52,57,60,65,69,72,74,75,77].These two methods have a theoretical similarity such that they use additional quantum computational operations to mitigate quantum noise effects, however, their treatments are technically different.In PEC, quantum noise effects are mitigated by constructing inverse processes of quantum noises comprised of the sixteen basis operations and quasiprobabilities called recovery operations.In other words, both original circuits and additional quantum circuits are probabilistically generated owing to quasiprobabilities.Suppose that the recovery operation for the quantum noise under consideration, which acts on a single-qubit state, is composed of N quasp nonzero quasiprobabilities.The elementary operations of the recovery operation (the terms composing the recovery operation) is inserted after each operation of U l and the maximum number of circuits for performing PEC is N Nqd quasp .Once all these circuits are run obeying the quasibrobabilities the quantum noise effects are completely mitigated (non-perturbavtive method with respect to the noise strength).We mathematically describe the recovery operation as follows.By writing the nonzero quasiprobabilities and the associated basis operations as η α and B α (α = 1, . . ., N quasp ), respectively, the recovery operation executed after the operation of U l is described in terms of these quantities as Nquasp αj l =1 η αj l B αj l , where N PEC q,l is the number of qubits on which the recovery operations in the lth layer act, B αj l is the basis operation acting on the jth qubit Q j , and η αj l is the associated quasiprobability of B αj l : note that η αj l include both positive and negative values.Furthermore, let us say that we compute an expectation value ⟨O⟩ with a repetition number N QC and accuracy ϵ and write an associated quantum mechanical variance by When we perform PEC under the same repetition number N QC and the accuracy ϵ the variance of the expectation value with PEC being performed, which we denote by where |η α | = exp(2bϵ err d) with b a positive constant and ϵ err is a gate error rate which is assumed to be gate independent (or a typical gate error rate) [27,49,52,60].Let us say that we generate M quantum circuits (we call M as PEC sampling number) obeying the quasiprobability and an elementary circuit of the M quantum circuits can be either the original circuit or one of the additional circuits.By writing the set of quantum circuits for PEC which is composed of N M [49].How many times a quantum circuit, say C a , appears in the M sampling depends on its quasiprobability.The computational resource, however, is finite in practice and in order to implement PEC in a real circumstance we take the sampling number M such that ⟨O⟩ PEC,M − ⟨O⟩ ideal ≤ ϵ.The sampling number M scales in ϵ as O ϵ −2 [49,52,60,74,75,77].On the other hand, our method QEM is conducted by the estimation of the kthorder quantum noise effect ∆ AD k ρ d•••1 and we subtract it from the noisy expectation value.In contrast to the quasiprobabilities, the coefficients which compose ∆ AD k ρ d•••1 such as ±1 and ± 1  4 are not probabilistic but deterministic (see Eq. ( 9) and the argument in Sec.I in the Supplementary Material, and furthermore, to evaluate ∆ AD k ρ d•••1 we deterministically prepare and execute the quantum-noise-effect circuit group whose size (the number of elementary circuits composing it) is O (dN q ) k .The quality of our QEM method gets better as we increase k and this corresponds to the increasement of M in PEC.Let us make numerical comparisons between these two methods for QAOA discussed in Sec.IV C and show them in Fig. 23.We perform the quantum simulations by taking the AD strength ϑ τ = 0.1 × a with a = 1, 2, . . ., 5. We take k = 1 (first-order perturbation theory) for our QEM method and for PEC we take M = 3dN q + 1, which is the number of quantum circuits we need to perform our first-order QEM scheme, and in this case d = 15, N q = 4 and M = 181.The reason why we take M = 3dN q + 1 is the following.Let us say that we conduct each QEM method with the same amount of quantum computational resource (under the same condition) and this can be regarded as the number of quantum circuits to be used for QEM.This is because, as mentioned previously, both methods are described as sums of expectation values generated by ensembles of quantum circuits consist of original circuits and quantum circuits including additional operations.By taking account this, next let us introduce an error defined by , which represents the absolute of the difference between the ideal expectation value and the expectation value with our QEM method being performed.On the other hand, we introduce an error δ PEC = ⟨ Ô⟩ ρ d•••1 − ⟨O⟩ PEC,M with setting M = 3dN q + 1.The numerical comparison between our method and PEC can be redescribed such that we perform each method using the common quantum computational resource which is the 3dN q + 1 quantum circuits and examine whether δ QEM is smaller than δ PEC or not: the method having smaller error can be represented as the better QEM method.Hereinafter, let us rewrite the expectation value with QEM and that with PEC as → ⟨ Ô⟩ QEM (ϑ τ ) and ⟨O⟩ PEC,M → ⟨O⟩ PEC (ϑ τ , M ), respectively, to express the noise-strength dependencies.Furthermore, we rewrite the ideal expectation value and the noisy expectation value as , respectively.Note that the operator O is chosen to be H C in Eq. (17).The recovery operation of AD acting on single-qubit states after the operation of U l is described by [74] where 1 j l is the identity operator acting on Q j , Z j l [ρ] = Z j ρZ j , (P |0⟩ ) j l is the reset operator acting on Q j , and ϵ AD (ϑ τ ) is the error rate given by ϵ AD (ϑ τ ) = sin 2 ϑτ 2 = 1− e −τ .Let us explain from the result in Fig. 23(a).Here we plot four types of quantities, the ideal cost function ⟨ Ô⟩ ideal (black dotted line), the noisy (no QEM) cost function ⟨ Ô⟩ noisy (ϑ τ ) (blue curve), the cost function with QEM ⟨ Ô⟩ QEM (ϑ τ ) (or-ange curve), and the cost function with PEC ⟨ Ô⟩ PEC (M, ϑ τ , l) with l = 1, . . ., N PEC samp , where N PEC samp is the repetition number of the PEC simulation under the same condition in terms of M and ϑ τ and ⟨ Ô⟩ PEC (M, ϑ τ , l) is the cost function acquired in the lth round of the PEC simulation.We note that N PEC samp is different from N samp introduced in Sec.IV.In contrast to our QEM scheme, PEC is a probabilistic theory and when we perform a PEC simulation for N PEC samp times and obtain N PEC samp data of the cost function ⟨ Ô⟩ PEC (M, ϑ τ , l), in general ⟨ Ô⟩ PEC (M, ϑ τ , l 1 ) ̸ = ⟨ Ô⟩ PEC (M, ϑ τ , l 2 ) for l 1 ̸ = l 2 .This is why we introduce the repetition number N PEC samp for the PEC simulation.On the other hand, ⟨ Ô⟩ QEM (ϑ τ , l 1 ) = ⟨ Ô⟩ QEM (ϑ τ , l 2 ) = ⟨ Ô⟩ QEM (ϑ τ ) for l 1 ̸ = l 2 and the repetition number N PEC samp is unnecessary for our QEM scheme.We compute ⟨ Ô⟩ ideal , ⟨ Ô⟩ noisy (ϑ τ ), ⟨ Ô⟩ QEM (ϑ τ ) with our original code while we compute ⟨ Ô⟩ PEC (M, ϑ τ , l) with the software package Mitiq [75] and Cirq [99], and the PEC simulations are done under the condition M = 181 and N PEC samp = 100 and the data points of ⟨ Ô⟩ PEC (M, ϑ τ , l) are plotted with green circles for each ϑ τ .To quantify which of the two cost functions, ⟨ Ô⟩ QEM (ϑ τ ) and ⟨ Ô⟩ PEC (M, ϑ τ , l), is closer to ⟨ Ô⟩ ideal , we introduce a measure defined by Like RT QEM in Eq. ( 13), the ratio RT PEC/QEM (M, ϑ τ , l) becoming greater than one implies that our QEM scheme is working better than PEC.By denoting the number of ⟨ Ô⟩ PEC (M, ϑ τ , l) data points satisfying RT PEC/QEM (M, ϑ τ , l) > 1 as N PEC/QEM , in Fig. 23(b) we show the ratio N PEC/QEM /N PEC samp .We see that the ratio RT PEC/QEM (M, ϑ τ , l) is greater than 50 percent for every ϑ τ .We also obtain the ratio N PEC/QEM /N samp = 0.99 for the realistic noise strength ϑ τ = 0.045 in the current superconducting qubit devices, as discussed in Sec.??.
In addition to the ratio RT QEM (M, ϑ τ , l) in Eq. ( 22), we compute mean squared errors defined by Note that σ 2 QEM (ϑ τ , N PEC samp ) = σ 2 QEM (ϑ τ ) (no N PEC samp dependency).By using the variances in Eq. ( 23) we can re-describe the comparison between the quality of our QEM and that of PEC as the comparison of the magnitudes of the two variances, i.e., the method exhibiting smaller variance has the better QEM quality.We show this numerical result in Fig. 23(c Note that we have used N PEC samp N QC data of expectation values ⟨ Ô⟩ PEC (M, ϑ τ , N QC , l) to compute σ 2 PEC (M, ϑ τ , N PEC samp , N QC ).We do not observe any significant difference in the variance of PEC between the finite and infinite N QC cases, as can be seen by comparing Figs.23(c) and 23(d).However, due to the finite N QC = 2 10 , we observe a slight increase in the variance of our QEM in the range 0.1 ≤ ϑ τ ≤ 0.3.Nonetheless, our QEM still exhibits As a result, the quality of our QEM method outperforms that of PEC under such conditions.Let us examine for realistic cases when M is sufficiently larger than 181 and we consider that there are two types of effects.The first one is that, owing to the concept of PEC, ⟨O⟩ PEC (ϑ τ , M ) gets closer to ⟨O⟩ ideal owing to an increasement of M , which is a positive effect.The second one, which is a negative effect, is that the total amount of errors associated with the additional operations gets bigger by increasing M which makes ⟨O⟩ PEC (ϑ τ , M ) to become farther from ⟨O⟩ ideal .At some point, say M = M c , we consider that the second (negative) effect gets larger than the first (positive) effect because the second effect is not to be mitigated and gets bigger.As a result, ⟨O⟩ PEC (ϑ τ , M ) becomes farther from ⟨O⟩ ideal for M > M c .It has been shown in [60] that the errors associated with the additional operations can be mitigated by combining with ZNE.The necessity of ZNE, however, implies that we need an additional computational resource to mitigate such errors.In contrast, our QEM method is conducted self-consistently such that the quantum noise effects on both the original circuit and the quantumnoise-effect circuit groups are mitigated.In other words, we do not need additional resources to mitigate the quantum noise effects on the quantum-noise-effect circuit groups.We consider that such a self-consistency is the advantage of our method compared with PEC.Furthermore, the size of the quantum-noise-effect circuit group (the computational cost) is practically polynomial in N q d while the number of quantum circuits for performing PEC is N Nqd quasp (exponential in N q d), and thus the practical computational cost of our QEM method is lower than that of PEC.In total, we consider that our QEM method is superior to PEC even for sufficiently large M.

B. Comparison with ZNE
Next, let us make a comparison between ZNE [26, 27, 49-52, 73, 75, 77] and our method.The similarity between ZNE and our method (as well as PEC) is that both of these methods require additional quantum circuits.In ZNE, first a quantum mechanical expectation of an physical observable, say O phys , is measured which includes a quantum computational error given by a noise strength (or an error rate) γ 0 .Then, we create extra quantum circuits which ideally yield the same expectation value of O phys as the original quantum circuit does but with noise strengths larger than γ 0 .Such boostings of the noise strengths can be done, for instance, by insertions of identity gate operations [27,50,73].Let us denote the original circuit by C org , which is subject to the quantum noise with the strength γ 0 , while we denote the extra circuits subject to noises with strengths γ j by C j with j = 1, 2, . . ., N ZNE , i.e., for ZNE N AQC = N ZNE .Here we labeled the extra circuits C j so that γ 0 < γ 1 < • • • γ NZNE .The highest order of the quantum noise effects which can be mitigated is N ZNE and this is one of the powerful advantage of ZNE.The difference between our method and ZNE, on the other hand, is that in our method (as well as in PEC and ESD and VD) QEM tasks are performed by quantum computational operations (gate operations and measurements) whereas in ZNE quantum computers are used for calculating expectation values while the QEM tasks are done by classical computers.The advantages of ZNE are that it does not require ancilla bits and QEM can be performed without knowledge of quantum noises.There is a drawback, however, that its quality becomes poor when γ 0 is too big.Furthermore, in experiments we need to obtain the values γ 0 , γ 1 , • • • , γ ZNE with high precision.On the other hand, our method requires both the ancilla bits and knowledge of quantum noises but can be applied to quantum noises for arbitrary noise strengths although higher computational cost is demanded for doing higher-order perturbation and both the values of decay rates and gate operation times are needed with high precision.Next let us compare the two methods with respect to the depth of a quantum algorithm and a coherence time with which we identify as a T 1 time.Here we assume that all single-qubit and two-qubit gate times are identical and denote the single-qubit time and the two-qubit gate time by t g 1 and t g 2 , respectively.Moreover, we assume that all T 1 times are identical with respect to qubits.Let us first consider from the conduction of Richardson-extrapolation-based ZNE and for simplicity we make an approximation t g 1 ≪ t g 2 so that the depth of a quantum algorithm d can be identified with the depth of (the number of the layers of) two-qubit gates to be implemented.Let us say that we enhance the noise strength γ i by the factor c i = r i with taking r i = 1 + 0.1 × i with i = 1, 2, . . ... Such a situation can be created by inserting an identity operation for the time 0.1×i×t g 2 after the operation of each U k .In this way, we can effectively make the operation time of U k to be r i t g 2 or we can effectively make the decay rate to be r i γ i while the operation time of U k to be t g 2 .The operation time to execute the quantum algorithm 2 and it must satisfy the condition r i dt g 2 < T 1 .For superconducting qubits t g 2 ≈ 100 nsec and T 1 ≈ 100 µsec [10,11,86] and under such conditions we obtain r i d < 1000.For d = 800, we have r 1 d = 880, r 2 d = 960, r 3 d > 1000, which means that we are able to perform QEM up to second order: For d = 900 we have r 1 d = 990, r 2 d > 1000, which means that we are able to perform QEM up to linear order, and for d > 910 we become unable perform QEM.For larger c i the upper limit of d such that QEM is valid gets smaller.From these considerations, we see that we can perform highorder QEM for small d while for large d we can only perform low-order QEM and for sufficiently large d we cannot apply QEM.Such a characteristics comes from the d dependence of the total additional identity operation time.Furthermore, the exponential extrapolation also becomes invalid for large d since it is only effective for small error rate ϵ err such that ϵ err N q d = O(1) [27]: with the same reason for the order estimation of C 2 PEC we consider that the condition of the exponential extrapolation to become invalid is not ϵ err N g = O(1) but ϵ err N q d = O(1).Next, let us consider the conduction of our QEM.The time for performing our kth-order QEM method is at most dt g 2 + (2t g 2 + t g 1 )k ≃ dt g 2 + 2.1kt g 2 (we assume that all additional operations are P 1 ), where we have set t g 1 ≃ 0.1t g 2 [50,86].Thus, the operational time for implementing the k additional operations, which is 2.1kt g 2 , does not depend on d.For d = 800, 900, and 910 the highest orders of quantum noise effects which we can mitigate are 95, 47, and 42, respectively.Consequently, we consider that for small d satisfying r ZNE dt g 2 < T 1 the quality of ZNE is better while for large d our method has a better quality and such a tendency does not change even T 1 times are extended and the gate times t g 2 get shorter.At the end, we mention that ZNE is not applicable to time-dependent noise [27,49,50] whereas our method is by reformulating τ as a function of time.

C. Comparison with ESD and VD
Finally, we make a comparison between our method and ESD [67,76,77] and VD [66,77].Since these two are similar approaches we bring them together and abbreviate it as ESD/VD.Like PEC and our method, ESD/VD is composed of gate operations and/or quantum measurement on an ancilla bit.ESD/VD has two advantages, it can be applied to various types of quantum noise and additional quantum circuits are unneeded.The procedure of ESD/VD is comprised of three parts, a preparation of n copies of an original circuit which generate (ρ real d•••1 ) ⊗n , a performance of derangement operation, which is a generalization of SWAP operation (in [66] it is called cyclic shift operator), and an operation of controlled-O operator with O denoting the physical operator of which we take an expectation value [67,76]: In [66], the scheme which does not use an ancilla-assisted measurement has been presented and in this case N a = (n − 1)N q (see Table I λ a |ψ a ⟩⟨ψ a |, where the eigenvalues satisfy the descending order λ dom > λ 1 • • • > λ 2 Nq −1 .Due to such a process, we obtain the expectation value ⟨ψ dom |O|ψ dom ⟩.In ESD/VD, there are two problems which need to be handled with, the mismatch between the ideal state ρ d•••1 and the dominant eigenvector state ρ dom = |ψ dom ⟩⟨ψ dom | which is called coherent mismatch or noise floor, and the mitigation of the quantum noise effects associated with the operations of the derangement operator and the controlled-O operator.It was argued in [67] that the quantum computational error coming from the coherent mis-match and that due to the quantum noise effects associated with the operations constructing the quantum circuit for ES-D/VD, which we call ESD/VD circuit, can be mitigated by combining with ZNE.As argued in the previous discussion for ZNE, we consider, however, that such a hybrid scheme works provided that the composite circuit of the original circuit and the ESD/VD circuit can be executed within a coherence time of a qubit.Furthermore, we need large numbers of qubits to perform ESD/VD and this limits sizes of overhead.On the other hand, our method needs the additional quantum circuits (quantum-noise-effect circuit group) but both the errors associated with the operations of U k and those associated with the insertions of the additional operations implemented on the quantum-noise-effect circuit group are self-consistently mitigated and saves an additional computational resource for mitigating the errors associated with the additional operations since our scheme does not need to be combined with other methods for doing this.Moreover, the number of ancilla bits which we need to perform k-th order QEM is k, which does not enlarge the overhead that much.Such properties, in principle, enable us to perform QEM for both short-term (NISQ algorithms) and long-term algorithms.

VII. CONCLUSION AND OUTLOOK
In this paper, we have established our QEM scheme for reducing the quantum noise (decoherence) effects on the singlequbit states which occur during the gate operations.We have formulated it as the perturbation theory with respect to the noise strength (in the case of AD effect it is τ ), which are evaluated by the gate time and decay rate (T 1 time and/or T 2 time), and is represented by the ensemble of quantum circuits, namely the quantum-noise-effect circuit groups.The numbers of quantum circuits composing the quantum-noise-effect circuit groups are polynomial with respect to the product of the depth of the quantum algorithm under consideration and the number of register bits, which can be considered that the conduction of our QEM scheme is not so high-cost computational performance.To demonstrate the validity of our QEM scheme, we have performed the noisy quantum simulations of the qubits under the AD effects for four types of quantum algorithms based on the linear-order perturbation theory.It is to be noted that before we conduct our QEM scheme, we need to be careful with if the expectation values on which we are aiming to perform QEM do not include the linear-order term in τ or if they are equivalent to the ideal simulation results like the computation of the expectation value of the identity operator.As long as these two are not the cases then our linear-order QEM scheme works as we have demonstrated in Sec.IV, and it is valid in a broad region of τ , which implies its effectiveness and powerfulness.Our QEM scheme can be generalized to error mitigation for other types of quantum noises including the generalized amplitude damping, the phase damping, the composition of these two, and the stochastic Pauli noises like the depolarizing channel.Furthermore, it can be extended to cases of error mitigation for higher-order quantum noise effects and once this is established we expect that we become able to perform quantum computations in high accuracy even for long-depth quantum algorithms.In Sec.V, we have discussed the experimental results of the implementation of our QEM scheme.Consequently, we have observed RT QEM > 1 for every quantum computation and our QEM scheme has worked for the IBM Q Experience processors.In this work, we have focused on the decoherence effects acting on the single-qubit states and established the QEM scheme for them.In real quantum devices, however, there exist many types of errors and complex quantum noise channels.We expect that by conducting an elaborate noise (device) characterization we become able to improve the quality of QEM.For instance, by combining our QEM scheme with other errormitigation techniques such as those for state preparation (initialization) and measurement and imperfections of gate operations or those for other quantum noise channels like crosstalk we anticipate to realize QEM with higher quality.
Our QEM scheme is solely conducted by gate operations and measurements on ancilla bits and can be applied to any type of quantum algorithm.Such a characteristic enables us to programmably perform high-accurate quantum computing solely by the quantum-computational operations (software manipulations).Furthermore, our QEM scheme can be performed with any type of quantum hardware such as solid-state systems and atomic-molecule and optical systems and with quantum devices of any generation, and both the computational cost whose order is polynomial in N q d and its accuracy can be coherently controlled.These three characteristics are the big advantages of our scheme.
One of the important outlooks of this work is quantum computations by large-scale quantum devices or future (nextgeneration) quantum devices using various types of quantum hardware such as superconducting circuits and ion-trap qubit systems.In such a case, we consider that we also need to take into account quantum noises acting on many-body quantum states such as collective quantum noises [88][89][90][100][101][102][103] and correlated noises [53,59,[104][105][106][107][108].Another important problem is the establishment of QEM schemes for timedependent quantum noises including non-Markovian quantum noises [109,110].
Our QEM scheme can be extended to mitigation of these complex quantum noise effects provided that they are formulated as groups of quantum circuits.When such formalisms are being constructed, we expect that we become able to realize QEM scheme which mitigates various types of quantum noise effects.We expect that this leads to conduction of quantum computing for big-size problems with high-quality results being obtained.We believe that this paves the way to realize high-quality quantum computing for application to problems in many branches of science and engineering including material science, quantum chemistry, combinatorial optimization problems, and machine learning using large-scale quantum computers.
IV. MEASUREMENT DATA OF T1 AND T2 TIMES In this section, we explain how to extract the experimental values of T 1 and T 2 times with showing our data plots, Fig. 2 (T 1 data points) and Fig. 3 (T 2 data points).The experiment is done by ibmq_belem (five-qubit machine).We note that for both cases the relaxation processes given by the relaxation times t relax a were generated by using a function (command) called "delay" [1].
Let us explain from the measurements of the T 1 times.In Fig. 2, we plot the experimental data of the expectation values ⟨P 1 ⟩(t tot,T1 a ) for the qubits Q j (j = 0, . . ., 4).Here we have taken t relax a = 2 × a (µs) with a = 1, . . ., 50 and N QC = 2 10 .We fit these data points by an exponential function ⟨P 1 ⟩(t tot,T1 and extract the values of T 1 .In Table I we list these experimental values.In addition, we show the data which are open to the public in the parentheses [1].As a result, our experimental data agrees well with them. Next, let us discuss the measurement data of the T 2 times presented in Fig. 3.Here we plot the experimental data points of ⟨P 1 ⟩(t tot,T2 a ) for the qubits Q j and we have taken t relax a = 16  45 a (µs) with a = 1, . . ., 300 and N QC = 2 10 .In this case, instead of a simple exponential decay curve, we fit the data points by a function where a 1,2 , ϕ 1,2 , and b are real numbers.The reasons we fit ⟨P 1 ⟩(t tot,T2 a ) with f T2 (t tot,T2 a ) in Eq. ( 16) are the following.First, a detuning, which is a difference between the qubit frequency and a driving frequency (frequency of a pulse which generates a gate operation), is not actually zero, and instead of an exponential decay the plots of an expectation value of P 1 must be fitted by a damped oscillation curve [1].Here we fit the the oscillation part by a cosine curve.When the detuning, however, is small it is not effective to fit the plots by the damped oscillation function since it is hard to distinguish whether the decreasing behaviors of the plots are coming from the exponential decay (decoherence) or from the oscillation.We avoid such an issue by enhancing the magnitude of the detuning explicitly.Thus, by writing the unitary operation for the measurement of the T 2 times by U meas,T2 , the actual series of the gate operations which have been exploited for the T 2 measurements are not the two Hadamard gates but U meas,T2 = H • I • R z (θ det ) • H with θ det = 16a 225 π (a = 1, . . ., 300).The rotation operation R z (θ det ) plays the role of the enhancement of the detuning whereas I represents the relaxation part and is generated by the delay command as mentioned above.Second, we have examined that the fitting by a damped oscillation function with an oscillation part given by a single type of amplitude, frequency, and phase is not a good fitting, and instead, using the function f T2 (t tot,T2 a ) in Eq. ( 16), which is given by two types of amplitudes, frequencies, and phases, shows much better fitting.We consider the physical reasoning for this as follows.In [2], two different two-level systems having opposite charge parities induced by quasiparticle tunneling have been observed in a transmon qubit system, and in such a case there exist two different frequencies of the two-level systems.Since the ibmq_belem is a transmon qubit system [3,4] we consider that the similar phenomenon has been observed in our experiment and we fitted the data points with the function f T2 (t tot,T2 a ) in Eq. ( 16).As a result, we obtain the values of T 2,j .In Table I, we display the experimental data of T 2 as well as the data open to the public [1] written in the parenthesis.The T 2 times for Q 1 and Q 2 show good agreement whereas those for the other qubits have moderate differences although the orders are the same.

V. QUANTUM DEVICE PROPERTIES
In this section, we present the data of the physical properties of ibmq_belem and ibm_perth (seven-qubit machine) which is open to the public [1]: In Fig. 4, we display illustrations of spatial configurations for qubits in these machines.First, we list the single-qubit properties of ibmq_belem, qubit frequencies and gate operation times in Table II  denote the gate times of identity gate, virtual Z gate, √ X gate, and X gate, respectively.In addition, we present the two-qubit gate properties, CNOT operation times denoted by t CX .The symbol Q i , Q j (i, j = 0, 1, 2, 3, 4 with i ̸ = j) represent nearestneighboring qubits Q i and Q j .The gate operation times in Table II

FIG. 1 .
FIG. 1. Schematic of our proposed QEM method.The original circuit represented by the blue rectangle (left side) describes the quantum circuit for the quantum algorithm to be ran and is composed of the unitary operations U k with k = 1, . . ., d and d denotes the depth of the quantum algorithm.It yields the expectation value ⟨ Ô⟩ ρ real d•••1

FIG. 2 .
FIG. 2. Illustration of quantum computation under AD effects represented as a quantum circuit.Here we show it for d = 3 and Nq = 2.The symbol E AD Q rj expresses the occurrence of AD effect on the register bit Qrj (the j-th register bit).

FIG. 3 .
FIG. 3. (a) Schematic of AD-effect circuit A. When we set ϑ = π and post-select the measurement result of the quantum state of Qa to be |1⟩Q a , we realize the operation of σ− on Qr0.(b) Schematic of AD-effect circuit B. By setting ϑ = π and post-selecting the measurement result of the quantum state of Qa to be |0⟩Q a , we the operation of P 1 on Qr0 is created.

FIG. 4 .
FIG. 4. (a) Schematic of an elementary quantum circuit in the ADeffect circuit group given by the AD-effect circuit A which generates d l=k+1 U l • σ− ρ k•••1 σ+ •
6(a) and 6(b) denote the expectation values of Ô while those in Figs.6(c) and 6(d) describe the ratio RT QEM given in Eq. (13): Figs.6(a) and 6(c) are the results for Ô = X while Figs.6(b) and 6(d) are those for Ô = Z.The dotted lines in Figs.6(a) and 6(b) describe the expectation values without QEM being performed whereas the solid lines represent the expectation values with QEM being performed.For both the dotted and solid lines the red, blue, and green plots in Figs.6(a) and 6(b) are the computational results of ⟨X⟩ and ⟨Z⟩ for d = 2 3 + 1, 2 4 + 1, and 2 5 + 1, respectively.The black dotted lines are the results of the ideal simulations.We have computed

FIG. 6 .
FIG. 6.Quantum simulations for the quantum algorithm U QC pre1 = H ⊗d−1 • X. Plots in (a) and (b) are the results of QEM for the expectation value of X and Z, respectively.The dotted lines are the expectation values without QEM while the solid lines represent the expectation values with QEM.The black dotted lines are the results of the ideal simulations.In (c) and (d), we show the ratio RTQEM for the expectation values of X and Z respectively.The red, blue, and green plots are the results for d = 1 + 2 3 , d = 1 + 2 4 , and d = 1 + 2 5 , respectively.

FIG. 7 .
FIG. 7. Quantum simulations for the quantum algorithm U QC pre2 = (UCH [Qr0; Qr1]) ⊗d−1 ) • XQ r0 ⊗ XQ r1 .Plots in (a) and (b) are the results of QEM for the expectation value of ZX and ZZ, respectively.The dotted lines are the expectation values without QEM while the solid lines represent the expectation values with QEM.The black dotted lines are the ideal simulation results.We plot the ratio RTQEM for expectation value of ZX and ZZ in (c) and (d), respectively.The red, blue, and green plots are the results for d = 1 + 2 3 , d = 1 + 2 4 , and d = 1 + 2 5 , respectively.

FIG. 9 .
FIG. 9.Quantum simulation results for the quantum algorithm U QAA in Eq. (14).Plots in (a) and (b) are the results of the probabilities P110 (probability of |110⟩) and P111 (probability of |111⟩), respectively.The dotted black lines are the ideal simulation results whereas the blue and orange curves are the noisy simulation results without QEM and the ones with QEM, respectively.All of them are obtained by our original code.The blue and orange circles are the noisy simulation results without QEM and the ones with QEM, respectively, and they are obtained by our Qiskit code.Plots in (c) and (d) are the results of the ratio RTQEM for P110 and P111, respectively.The black curves are obtained by our original code while the red circles by our Qiskit code.For each ϑτ we have plotted 100 circles in (a) -(d), i.e., Nsamp = 100.
(a) the ideal minimum value of the cost function is −4.0, and we obtain the maximum cut number four by multiplying minus one.Consequently, our QEM scheme works for the noisy QAOA simulation.

FIG. 20 .
FIG. 20.Plots in (a) and (c) are quantum computation results obtained by ibm_perth and (b) and (d) are quantum simulation results.The quantum algorithm which has been run is U QC pre1 .Both (a) and (b) are the results of expectation values of Z1 and (c) and (d) are the ones of RTQEM.All the quantum circuits have been executed under NQC = 2 14 , Nsamp = 5.The horizontal axes represent d = n rep + 1, where n rep = 16, 32, 64.The quantum computation has been conducted at 07:17, 09/26/2023.

U QC imp, 2 :
Figs. 22(a) and 22(c) are the quantum computation results and Figs.22(b) and 22(d) are the quantum simulation results.As similar to the labeling in Fig. 21, Figs.22(a) and 22(b) (Figs.22(c) and 22(d)) are the results of ⟨Z 1 Z 3 ⟩ (RT QEM ), and all the vertical axes represent d = 3n rep + 1 with n rep = 5, 10, 15.Overall, the characteristics of these results are similar to those for ibm_perth.Here we plot two types of expectation values with QEM which are indicated by the orange open circles and the green squares and two types of RT QEM by the red open circles and the green squares.The quantities plotted by the orange and red open circles have been calculated by using the data of the T 1 and T 2 times obtained at 05:08 whereas the others obtained by the data taken at 12:33: see Table
(Table IV-Table VII) list the physical properties of ibmq_belem (ibm_perth).Here let us briefly explain the gate properties.The native gates of these machines are the following, identity gate, virtual Z gate, √ X gate, X gate, CNOT, and reset operation (reset into |0⟩).The Hadamard gate is decomposed into a series of the native gates as . Here t id , t Rz , t √ X , and t X

TABLE II .
and Table III are used for the implementation of our QEM for the quantum algorithm U QC imp,2 : see the caption of Fig. 22 in the main text.Similarly, in Table IV-VII we list the physical Qubit ω01 (GHz) t id (ns) tRz (ns) t √ X (ns) tX (ns) List of the single-qubit properties of ibmq_belem which are open to the public [1].The above data was taken at 05:08, 09/11/2023.

TABLE III .
[1]t of the two-qubit properties of ibmq_belem which are open to the public[1].The above data was taken at 05:08,09/11/2023.properties of ibm_perth: TableIVand Table V present the data of the single-qubit properties while Table VI and Table VII show those of the two-qubit properties.We have used the data of the gate operation times taken in 09/26/2023 for the implementation of our QEM for the quantum algorithm U QC pre1 (Table V and Table VII) while we have used the ones taken in 09/21/2023 for U QC imp,2 (Table IV and Table VI): see also the captions of Fig. 20 and Fig. 21 in Sec.VB in the main text.Qubit ω01 (GHz) t id (ns) tRz (ns) t √ X (ns) tX (ns)

TABLE IV .
[1]t of the single-qubit properties of ibm_perth which are open to the public[1].The above data was taken at 07:32, 09/21/2023.

TABLE VII .
[1]t of the two-qubit properties of ibm_perth which are open to the public[1].The above data was taken at 04:30, 09/26/2023.