A novel approach to noisy gates for simulating quantum computers

We present a novel method for simulating the noisy behaviour of quantum computers, which allows to efficiently incorporate environmental effects in the driven evolution implementing the gates acting on the qubits. We show how to modify the noiseless gate executed by the computer to include any Markovian noise, hence resulting in what we will call a noisy gate. We compare our method with the IBM Qiskit simulator, and show that it follows more closely both the analytical solution of the Lindblad equation as well as the behaviour of a real quantum computer, where we ran algorithms involving up to 18 qubits; as such, our protocol offers a more accurate simulator for NISQ devices. The method is flexible enough to potentially describe any noise, including non-Markovian ones. The noise simulator based on this work is available as a python package at this link: https://pypi.org/project/quantum-gates.


I. INTRODUCTION
Quantum computers are on the way; currently they manage between dozens and hundreds of qubits [1][2][3][4][5], which does not sound as an impressive number, yet it is already good enough to perform interesting tasks [6,7].As powerful as they promise to be, quantum computers are far from being ideal: since, as for any quantum system, they can hardly be isolated from the surrounding environment, they are prone to errors, which limit their capabilities.Like in the classical case, error correcting schemes have been developed [8][9][10] and first tests have been performed [11,12], but to be implemented they require to the least thousands qubits, which are not available; for the time being, we have to cope with errors.
This stage of development is referred to as Noisy Intermediate-Scale Quantum (NISQ) [6,13] era; the major aim of the research during this near-term period is to maximize the computational power of current devices in view of the long-term goal of fault-tolerant quantum computation [1].
It is clear that NISQ computers require a good understanding of how noises affect quantum circuits and, in order to do so, a proper modeling of the noises is needed.This requires essentially two major tasks: understanding the major sources of noise affecting the qubits, and writing better algorithms for simulating a given noise model on a classical computer.The present work deals with this second task.
To date, the simulations of noisy digital gate-based quantum computers is implemented by adding appropriate quantum operations before and after each ideal gate [14][15][16]: schematically, and working with the density matrix formalism, if an ideal (unitary) gate G is supposed to be executed, the noises affecting it are modeled by adding appropriate operations E 1 (E 2 ) mimicking the noise, before (after) the gate: Such a modeling completely decouples the action of the controlled operation generating the gate G from that of the environment.This approximation works well if G acts almost instantaneously with respect to the noise, i.e. if the gate time t g required to implement the gate is much smaller than the characteristic time scales of the systemenvironment interaction.For instance, in IBM's superconducting devices [17] t g ∼ 10 −8 s, while typical environmental effects such as relaxation and phase damping have characteristic times of order T 1 , T 2 ∼ 10 −4 s.This justifies why this approach has been implemented by the majority of available noise simulators of NISQ computers (see appendix H).
Yet this approach has some limitations.By separating the action of the gate from that of the noise, it does not represent a faithful description of what happens inside a computer, where the controlled action on the qubit(s) generating the gate and the environment act simultaneously and potentially affect each other.Therefore it is expected not to be fully accurate in describing a NISQ computer, especially when the number of gates and qubits is relatively large, which is actually the regime where simulations are more interesting.
In this article we propose an alternative approach, where the noise is integrated into the logical gates, in the sense that the resulting noisy gate is computed by solving for the dynamics generating it, with additional terms describing the noise added to it: where in general G ̸ = E 2 • G • E 1 , and under standard assumptions (e.g., Markovianity) it gives an analytic expression for the solution of the Lindblad equation obtained with perturbative methods.Now G captures, within the limits of validity of Lindblad's equation, the entire physics occurring during the execution of each gate; not only it offers a more accurate description of the system and therefore a better protocol for circuit simulations, but also it helps to understand the different noises acting on the computer, especially in view of possible mitigation strategies.This new approach does not have any computational disadvantage with respect to (1).
As a note, Markovianity, which is the main physical assumption behind the Lindblad equation, and it is a very convenient working hypothesis, can be released in favour of more general noises [14,[18][19][20][21][22]; we will not touch on this possibility here, although the generalization of the approach here introduced is rather straightforward.
Both approaches (1) and ( 2) have a drawback if they are implemented at the density matrix level: the simulation will be slowed down quadratically as a function of the number of qubits.This drawback can be resolved for (1) by replacing the superoperations E 1,2 acting on the density matrix with suitable stochastic operations acting on the state vector [23,24]; in this way, the noisy algorithm becomes random and each single run of the simulation can be seen as a single run of the algorithm on the noisy quantum computer.The same strategy can be adopted for (2); one writes where G ξ is a stochastic gate, solution of a stochastic Schrödinger equation, incorporating both the controlled action generating the (otherwise ideal) gate G and the noise.Here ξ denotes a set of stochastic gaussian variables, and stresses the fact that G ξ , and hence |ψ ′ ξ ⟩, are random; we will omit to indicate ξ in the rest of the paper.Physical quantities are obtained by averaging over the noise.
The general procedure therefore is the following.Given a noiseless algorithm, the corresponding noisy one is obtained by replacing each ideal gate with a noisy gate.The resulting noisy algorithm, which is stochastic, is repeated for different realizations of the random variables, as if they were different runs on a physical quantum computer.This produces a statistics of outcomes, to be compared with those of a real computer, or to be used to predict the behavior of a future NISQ device.
As such, as already mentioned, a first application of our approach is to predict the behaviour of NISQ devices, their potentialities and limitations.But its use goes beyond the NISQ-era horizon: by offering a more accurate modeling of the noise, it allows to better understand the physics underlying the functioning of a quantum computer and to enforce appropriate error mitigation schemes [25][26][27].
The rest of the paper details this program.We present the noisy gates method by designing it on the IBM superconducting computers [17] although the approach is general and can be used to describe any NISQ quantum platform, once the native gate set and the proper noise model are defined.
The paper is organized as follows.In Sec.II we review the main noises affecting superconducting qubits, and how they are described within the Lindblad's formalism; in Sec.III, IV and V we present the general derivation of the noisy gates, specializing it to the native single and two-qubits noisy gates of IBM devices.In Sec.VI we compare the structure of our algorithm with that of IBM Qiskit.
In Sec.VII we present the results of the simulations, which test our algorithm against that of Qiskit in reproducing the solution of the Lindlbad equation, as well as the outcomes from current IBM quantum computers: the simulations show that the proposed method is more accurate and precise compared to that of Qiskit in reproducing the Lindblad equation, with an average improvement between 50% and 90% and more.
The improvement in simulating the real device fluctuates between 10% and 30%, because the underlying noise model is not accurate enough, and also because the devices are not really stable; for a large number of qubits it becomes even lower because the number of runs of the device, which are necessary to recover a good statistic, is too high.In both cases, this is not a limitation of our algorithm, but of the physical model describing the computer.
We conclude our paper with some general remarks and an outlook.

II. REVIEW OF THE NOISE MODEL
The noises which are more relevant in the functioning of superconducting devices have already been characterized in literature [15,16,28]; in this section we briefly present them.With good approximation they are described by a Lindblad dynamics [29,30]: here, H s is the Hamiltonian of the system which implements the ideal gate, and D(ρ) is a Lindblad term describing the effect of the environment.For convenience, we will describe the evolution with a time schedule s ∈ [0, 1], defined as s = t/t g , where t g is the duration of a gate.Apart from state preparation and measurement (SPAM) errors, which happen at the very beginning and very end, during the execution of an algorithm there are two main sources of noise, namely, depolarization and relaxation [28,31].The first, which can be ascribed to the imperfections of the device, tends to bring the state towards the totally mixed one, 1/ √ N , where N = 2 n and n is the number of qubits; for the single qubit, this can be modeled by the following Lindblad term [15,16], where σ 1 = X, σ 2 = Y, σ 3 = Z are the standard Pauli matrices and γ d ≥ 0 is the rate at which depolarization occurs.
The second type of noise is due to the interaction of the physical qubits with the surrounding environment; in particular, due the thermalization towards an equilibrium with the environment, energy exchanges occur.In the scenario of interest, this induces the decay of a qubit towards the ground state |0⟩, an effect which is also known as amplitude damping [15,16].This damping is characterized by a relaxation time T 1 , which identifies the scales at which the initial state decays towards |0⟩; it causes also a damping of the off-diagonal elements of the density matrix in terms of dephasing, which (if only amplitude damping is acting) has a characteristic time 2T 1 .However, at the same time also a contribution of pure dephasing must be taken in account, resulting in an effective dephasing rate 1/T 2 ≥ 1/2T 1 .When also T 1 ≥ T 2 holds (and this is the case of interest to us), the combined action of these two effects, that from now on we will refer to as relaxation or amplitude and phase damping, can be described by the following Lindblad term, where we use the convention σ ± = (X ± iY)/2 and P (1) = |1⟩ ⟨1| is the projector onto |1⟩; the coefficients are related to the characteristic times as γ 1 = t g T −1 We will consider both sources of noise together, meaning that the Lindblad term is D(ρ) = D d (ρ) + D R (ρ), which can be diagonalized in the canonical Lindblad form by standard procedures.Eventually one obtains the Lindblad term where the non normalized Lindblad operators are here, we set , and we defined the parameter ϵ = √ λ.As mentioned in the Introduction, in the case of IBM's superconducting devices the typical order of magnitude of the decoherence times is ∼ 10 −4 s; by contrast, the typical order of magnitude of the time to execute a gate is t g ∼ 10 −8 s, which is small compared to T 1,2 ; in particular, one has This justifies the perturbative expansion we will implement later.While terms of the form (7) describe the dissipation occurring at the single qubit level, one straightforward generalization to the multi-qubit case (the one we will consider in this work) is obtained via the direct sum where the upper index (k) indicates that the Lindblad term (7) acts on the k−th qubit.Such a generalization is based on the assumption that single qubit noises are dominating, therefore neglecting cross talks and correlated noises [32]; they can straightforwardly be implemented in our noisy framework, and they will be the subject of future research.We stress that through Eq. ( 9) we already account for the fact that (for instance, on IBM's devices) multiple-qubit operations are more faulty than single qubit manipulations: when entangling gates are performed, single qubit noises act together, and errors therefore amplify.Before proceeding, one further comment is in order.Casting the behaviour of a real quantum device in a theoretical model is a hard task, and the more accurate the model, the less general it is.As remarked in the Introduction, the purpose of this work is not that of finding the best noise model for a given quantum computer; rather, given a noise model, we are interested in the best way to simulate the device.The noise model we are considering here is therefore ultimately motivated by the fact that it is accurate enough to already give appreciable results in the simulations, but on the other hand it is also simple enough to efficiently enlighten our main points, and general enough to be readily extended to different platforms.It is understood that better results can be achieved only by specializing more the analysis on physical device to be considered.

III. GENERAL DERIVATION OF NOISY GATES
Let us consider the situation in which the computer executes a gate U g on a set of n qubits.This is achieved by driving the system with an Hamiltonian H s for s ∈ [0, 1], which will induce some unitary evolution U s , defined by iℏdU s /ds = H s U s , and such that U s=1 = U g .However, if noises and imperfections are taken in account, this coherent evolution is replaced by a partially non coherent one, which under the assumptions of Markovianity (and complete positivity) is described by a master equation of the form (4) discussed in the previous section, with the Lindblad term given, in our case, by ( 9) and (7), which needs to be solved in place of the Schrödinger equation.We recall here that the coefficient ϵ is small, ϵ ≪ 1.
In order to switch from the density matrix formalism to the state vector formalism, we perform a linear stochastic unraveling of the Lindbald equation [19,[33][34][35][36]; specifically we consider the following Itô stochastic differential equation for the state vector: [37] d where dW k,s are differentials of standard independent Wiener processes, i.e. stochastic infinitesimal increments such that E dW k,s = 0 and E dW k,s dW k ′ ,s ′ = δ k,k ′ ds.Eq. ( 10) is an unraveling of the Lindblad equation in the sense that the density matrix obtained by averaging the pure states |ψ s ⟩ ⟨ψ s | over the noise: is a solution of Eq. ( 4).In this sense, Eqs. ( 10) and ( 4) have the same physical content; the advantage of the stochastic unraveling is that it allows to work with Schrödinger-like equations for the state vector.
One key property of Eq. ( 10) is that it is linear, and therefore it allows to write the solution as |ψ s=1 ⟩ = N g |ψ 0 ⟩, where N g can be interpreted as a noisy random gate acting on the system.Since Eq. ( 10) in general does not preserve the norm of the state vector, the associated gate N g is not unitary; this is a consequence of the chosen unraveling: one could have chosen normpreserving unravelings [38,39], which however are not linear and therefore do not allow for a gate-like formulation.The lack of norm preservation is not a problem since at the statistical level, i.e. when the average over the noise is taken as in (11), one recovers the Lindblad equation, which is trace preserving.
In general, Eq. ( 10) cannot be solved in a closed form [37,40] except for few specific cases, for example when all operators commute.In Appendix A we show how an approximate solution to order O(ϵ 2 ) can be derived, which results in the following expression for the noisy version of a noiseless gate U g : where we defined the deterministic operator: and the stochastic one: Note that in Eqs. ( 13) and ( 14), L k,s = U † s L k U s are the Lindblad operators in the interaction picture, therefore the noiseless part of the dynamics U s and the noisy one given by the Lindblad operators L k do not factorize, as it might look from a naive understanding of Eq. (12).
As explained in Appendix A, we omitted the additional term −(ϵ 2 /2) (14), which in principle should contribute to order ϵ 2 ; this is legitimate because it is a nested Itô integral of non anticipating functions [37], and hence its stochastic average is 0. For this reason, it drops from all final averaged quantities, and therefore we can neglect it from the start.
Let us also point out that, in the cases of interest to us, the term (13) can always be exponentiated, so that we will always be able to directly calculate e Λ .
The only stochastic term entering the noisy gate N g is Ξ in Eq. ( 14), which is a function of several random variables ξ arising from the stochastic processes W k,s .Let us call L kij,s = L + kij,s + iL − kij,s the ij-th matrix element of the jump operator L k,s in the computational basis, divided in real (+) and imaginary (−) part, respectively.Then, each entry of the stochastic matrix is of the form where we defined the random variables which, being Itô integrals of deterministic functions, are all normally distributed with zero mean, E ξ ± kij = 0, and variances Moreover, one can easily check that they are correlated with each other as The random variables giving Ξ its stochastic character may be defined in several other ways, and the best choice depends on the specific case of interest.In this section we presented one general strategy for defining them, but in practice this lead to an over estimation of the actual number of random variables needed.By straightforwardly counting, one has at most 2N 2 (N 2 −1) real gaussian random variables for a noisy gate acting on n = log 2 N qubits, each random variable being correlated with at most other 2N 2 − 1 ones.In practice, however, we immediately point out that one shall expect neither the number of random variables, nor the number of correlations between them to really follow this scaling.This is mainly due to the fact that real quantum computers usually perform single and two qubit native gates, and single qubit noises are dominating.For instance, given (9), one can upper bound the number of random variables by ∼ 6N 2 log 2 N .In the following sections, as we go through the construction of the native set of noisy gates for IBM's quantum computers, we shall make this claim more clear.Note that our derivation works for any choice of the starting Lindblad master equation, meaning that any Markovian noise model can be treated.In particular, while in this paper we specialize on the noise model described in the previous section, one can add devicemotivated modifications (such as correlated noises and leakages to upper levels in the case of IBM's platform); modifications of this kind are left to future research-as also the generalization of our derivation to non Markovian situations.
More details on the difference between our perturbative approximation and the one used in the standard approach (1) can be found in appendix B.

IV. SINGLE QUBIT NOISY GATES
IBM's superconducting devices implement single qubits operations with unitaries of the form U(θ, ϕ) = e −iθRxy(ϕ)/2 , where we set R xy (ϕ) = cos(ϕ)X + sin(ϕ)Y; such gates are achieved by driving the system with the Hamiltonian [28,41] applied for a time s = 1 [42].The Hamiltonian is driven by time-dependent pulses [28], so that in Eq. ( 17) one should actually consider θ → ω s , and set In this work we consider constant pulses for simplicity, being the generalization to general functions rather straightforward.It should be noted that the functional form of ω s affects the action of the noises on the system, meaning that different pulse shapes might lead to smaller noise effects, i.e. error mitigation; this is a question left for future research.
The task now is to derive the noisy gates N(θ, ϕ) corresponding to the unitaries above, when depolarization and relaxation errors are both taken in account during the evolution.
We begin by computing the evolution of the jump operators in the interaction picture, obtaining the expressions: and where we defined R(θ, ϕ) = cos(θ/2)Z + sin(θ/2)R xy (ϕ) and for a generic angle α we set ᾱ = α + π/2.Then, based on Eq. ( 12), we compute the deterministic, non unitary term Λ(θ, ϕ).Since in the interaction picture the evolution is unitary, one sees that the term corresponding to k = 3 is always vanishing, and one has Λ(θ, , and after integration we get where R(θ, ϕ) = cos(θ/2)Z + sin(θ/2)R xy (ϕ) and φ = ϕ + π/2, so that one has such an expression can be readily exponentiated, leading to where we defined Next, we turn to investigating the stochastic term, Ξ(θ, ϕ).Here, it is convenient to define the following real stochastic variables: 23) whose variances are: while the correlations are: moreover, we define Summing all terms and re-arranging them conveniently, we arrive at the following expression where we defined the following set of complex stochastic coefficients: Since these quantities are all combinations of gaussian random variables with the correlations previously discussed, they can be efficiently sampled with known algorithms; then, the stochastic matrix ( 27) can be assembled and numerically exponentiated.Multiplication by the deterministic term (22) and then by the noiseless gate U(θ, ϕ) eventually lead to the noisy gate N(θ, ϕ) for the single qubit, which, as shown only depends on 8 correlated gaussian variables.
V. TWO-QUBIT NOISY GATES On IBM's quantum chips, two qubit gates are implemented by a driven cross resonance [28,41,43]; labeling with an upper index the qubit each operator acts on, this consists in the execution of the unitary U (1,2) (θ, ϕ) = e −iθZ (1) ⊗R (2)  xy (ϕ)/2 , which can be realised by driving the composite system with the Hamiltonian for a duration s = 1, where, from now on, the tensor product symbol will be dropped, unless otherwise specified.In the proposed approach we take in consideration only noises acting on single qubits, so that the Lindblad term reads where now ρ is the two-qubit statistical operator.The procedure for calculating the noisy gates is the same as in the single qubit case.
First, we compute the Lindblad operators on the first qubit (i = 1) in the interaction picture: s = Z (1) remains constant as it commutes with the Hamiltonian.For the second qubit (i = 2), one has s = R(2sθ, φ), where we defined for convenience The deterministic term Λ(θ, ϕ), see Eq.( 13), can be calculated straightforwardly, leading to (36) notice that again this term can be exponentiated analytically as all the terms involved commute; in particular, one has where F (θ) is the same function defined in the single qubit case.
In order to efficiently write the stochastic term Ξ(θ, ϕ), it is convenient to define, in analogy with the single qubit case, the gaussian random variables k,s sin(sθ), (38) and whose correlations are straightforward to calculate and mimic those already seen in Sec.IV.Then, we can separate Ξ(θ, ϕ) in two parts as Ξ (1) (θ, ϕ) + Ξ (2) (θ, ϕ); the first is equal to while the second part reads where we defined and Again, as in the single qubit case, the stochastic matrix Ξ(θ, ϕ) can be assembled by combining gaussian random variables, and hence it can be efficiently sampled and numerically exponentiated; this, combined with the term U(θ, ϕ)e Λ(θ,ϕ) , gives the noisy gate for two qubits.

VI. COMPARISON OF THE ALGORITHMS
It is instructive to compare the structure of our approach to noise simulation with that of noise simulators based on the standard approach in Eq. ( 1).As shown in appendix H all relevant quantum computing frameworks implement such standard approach, and we chose IBM's Qiskit as term of comparison since it is the most developed one; in the next section we will compare also their performances in simulating the Lindblad equation as well as a real quantum computer.
Both methods rely on the state vector formulation, with important differences though.According to the Qiskit documentation [24,44] the noises are implemented by Kraus maps, which in the density matrix formalism read: The map can be unraveled as a stochastic map on the state vector by imposing that, at a given time, |ψ⟩ changes randomly as follows: with probability: The associate pseudo code is reported in Alg. 1.

Algorithm 1 Qiskit Simulation
Input: Initial state |ψ0⟩, a noiseless circuit C = {U (1) , ..., U (ng ) } composed by ng gates U (i) and number of samples Ns for 0 ≤ k ≤ Ns do while The time complexity of Alg.1 is primarily determined by the matrix vector multiplication step, exhibiting a complexity of O(2 2n ), where n is the number of qubits.The space complexity is dominated by the storage of the state vector and it scales as O(2 n ).It has to be noted that when the Kraus operators are not unitary, as for relaxation, one needs to store the intermediate state vectors, which are necessary in order to compute the probabilites in Eq.( 46).This operation has the same time and space complexity as those of the the previous step.(This can be avoided for mixed unitary error channels: probabilities are known and independent of the current state.) Our noisy gates simulation instead is based on the algorithm summarized in Alg. 2.

Algorithm 2 Noisy Gates Simulation
Input: Initial state |ψ0⟩, a noiseless circuit C = {U (1) , ..., U (ng ) } composed by ng gates U (i) and number of samples Ns for 0 ≤ k ≤ Ns do map a noisy circuit C = {N (1) , ..., N (ng ) } on C sample stochastic processes ξ inside noisy gates N (i) The time complexity of Alg. 2 is again O(2 2n ), determined by the matrix vector multiplication step.Analogously, the space complexity is O(2 n ).We notice that in Alg. 2 there is no need to perform the scalar product in Eq. (46).Moreover, all optimization to reduce the time complexity that are possible for the first step of Alg. 1 are also possible for Alg. 2. Finally, both algorithms perform samples of random numbers, but this operation has a constant scaling.

VII. SIMULATIONS
We now study the performances of our noisy gates method, and compare them with those of Qiskit's simulator [44].First, in subsection VII A we test the two approaches against the solution of Lindblad equation ( 4), by studying a repeated application of IBM's native gate set.Then, in subsection VII B we compare the predictions of both methods with the behaviour of an actual quantum computer, by running the inverse QFT algorithm on the IBM's quantum processors ibmq kolkata and ibmq oslo.In appendix G we perform the same analysis by running the GHZ algorithm on ibm oslo.All simulations are performed by using the noise model described in section II (see also appendices C and D).The implementation of the work proposed in this paper is open source and available as a python package at this link.It allows the user to run noisy simulations.

A. Comparison with the numerical solution of Lindblad equations
First, let us compare our method with the one implemented in the Qiskit simulator for the task of simulating the Lindblad equation.To this purpose, we simulate the same Lindblad equation with both methods, obtaining the density matrix ρ ng , from the noisy gates simulation, and the density matrix ρ ibm from the Qiskit simulation.We then benchmark the results with the density matrix σ obtained by directly solving numerically the Lindblad equation with Mathematica [45].We compare these density matrices by computing the Hellinger distances where the Hellinger distance is defined by with ρ kk (σ kk ) the diagonal elements of ρ (σ).Note that the Hellinger distance is a classical measure of the distance between the readout probability distributions: while it cannot be interpreted as a distance between quantum states (it does not take in account the coherences), it directly compares the concrete outputs of the real device, which are classical (the oucomes of Z measurements).In appendix F we also compute the fidelities We run the simulations on both single and two qubit gates.Considering the native gate set of IBM's quantum computers, {R z (ϕ), X, SX, CNOT}, we remind that R z (ϕ) are implemented as virtual gates [28,41], i.e. they are noiseless, and the CNOT gates are implemented by combining single qubit gates in Eq. ( 17) and CR gates in Eq. ( 31) [28,41,46].Moreover, X and SX gates are both rotations around the X-axis for different values of θ, see Eq. (17).Thus for our purposes, it is sufficient to simulate the X, CR and CNOT gates affected by noises.The fact that noises drive the system towards the maximally mixed state is the reason why the improvement decreases in time.The noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times, as the noise is dominant over the unitary evolution.In the interesting regime [0, 2000 • tg] before decoherence dominates, our improvement is always above 60%.
Single qubit simulations.We first simulate a repetition of X gates, each of which can be obtained by setting θ = π and ϕ = 0 in Eq. ( 17); we initialize the qubit in |0⟩ and we use the qubit noise parameters of ibmq manila (more details on the device can be found in appendix E).We evolve the state of the qubit for a time T = N t g , with N = 15000.In the upper panels of Fig. 1 we plot the time evolution of the population of the ground state, ρ 00 = ⟨0|ρ|0⟩, as obtained with the three methods.In the noiseless case, ρ 00 should oscillate between 0 and 1 with period 2t g , as at each step of t g a complete X rotation is performed; in the presence of noises, the oscillations are damped due to the relaxation of the qubit, while the depolarization drives probabilities towards the asymptotic value ρ 00 → 0.5.
Both our simulation and that obtained using Qiskit's simlulator qualitatively reproduce this behaviour.In Fig. 1 we have also highlighted with dashed vertical lines the characteristic times of relaxation and depolarization (see the caption); for times approaching these values the state is not a reliable quantum state anymore, as the density matrix becomes completely mixed.Given this consideration, in the lower plots we stop at N = 2000.
In order to inspect which of the two models reproduces more accurately and precisely the Lindblad evolution, we have run 100 independent simulations with both the noisy gates simulator and the Qiskit simulator, computing for each run the Hellinger distances H ng σ , H ibm σ .We computed the means over the 100 independent simulations, Hng σ , Hibm σ and the standard deviations ∆H ng σ , ∆H ibm σ .These quantities are shown in the lower panels of Fig. 1.During the relevant time interval [0, T ] the Hellinger distance of the noisy gates simulator is closer to zero, than that obtained with the Qiskit simulator.Both results are compatible within the error bars, however the standard deviations associated to the noisy gates simulations are significantly smaller than those associated to the Qiskit simulations, as also highlighted in the inset of Fig. 1 (e).We notice that the difference between Hng The fact that noises drive the system towards the maximally mixed state is the reason why the improvement decreases in time.The noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times, as the noise is dominant over the unitary evolution.In the interesting regime [0, 100 • tg] our improvement is always above 88%.
increases.The relative improvement is shown in Fig. 1 (f).The fact that noises drive the system towards the maximally mixed state is the reason why the improvement decreases over time.The noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times.Indeed after such times the strength of the noise is dominant over the unitary evolution, or the Hamiltonian contribution is negligible with respect to the Lindblad term (see Eq.( 4)), which is the same in the two approaches.In the interesting regime [0, T ] our improvement is always above 60%.In appendix F we repeat a similar analysis for the fidelities.
Two qubits simulations.Next, we simulate a repetition of Cross-resonance gates as defined in Eq. ( 31), where we choose ϕ = 0 and θ = π.We initialize the system in the state |10⟩ and we use the qubit noise parameters of ibmq manila.In the three upper panels of Fig. 2 we show the time evolution of the entry ρ 22 = ⟨10|ρ|10⟩; the x-axis is normalized in terms of the two-qubit gate time t g .The two-qubit state goes asymptotically towards the completely mixed state as ρ 22 reaches the asymptotic value 0.25.The probability ρ 22 , which in the ideal case should flip between one and zero, is again damped over time by relaxation effects.Again, we have highlighted with vertical dashed lines the characteristic time scales of the noises, showing only the T 1 and T 2 values of the target qubit as representative values.The depolarizing error is the dominant one, spoiling the quantum state already after ∼ 100 CR gates; for this reason, the the lower panels we consider a total duration N ∼ 100.As before, we report the Hellinger distances, showing the different results of 100 independent simulations together with their mean and standard deviation in the three lower panels of Fig. 2.
As in the single qubit case, within the relevant time interval [0, T ] the Hellinger distances obtained with the noisy gates simulations are closer to zero than those obtained with the Qiskit simulator.However now, differently from the single qubit case, the two results are not compatible within error bars.Moreover the difference between Hng σ and Hibm σ is now of the order ∼ 10 −1 .This corresponds to a relative improvement in the range from 90% to 88% as time increases, shown in Fig. 2 (f).In the interesting regime [0, T ] our improvement is always above 88%.We notice that in Fig. 2 (e) the value of Hibm σ approaches that of Hng σ for times close to 100 CR gate times.The reason why this happens is the same explained above for the single qubit case.In appendix F shows the relative improvement.The fact that noises drive the system towards the maximally mixed state is again the reason why the improvement decreases in time: the noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times, as the noise is dominant over the unitary evolution.In the interesting regime [0, 100 • tg] our improvement is always above 55%.The upper subplots of panel (f) show the time evolution of the ρ11 entry of the density matrix for the same sequence of CR gates in Fig. 2 and the lower subplots show the time evolution of the ρ11 entry of the density matrix for the sequence of CNOT gates.Colors have the same meaning as for Fig. 1.For the CR gates, the Qiskit simulation of ρ11 is visibly different from the Lindblad evolution, thus explaining the higher improvement of the noisy gates simulation in the Hellinger distance in Fig. 2.
we repeat a similar analysis for the fidelities.
We then perform the analysis for a repetition of CNOT gates, for an initial state given by |10⟩ and qubit noise parameters of ibmq quito (see appendix E).We notice that in this simulation we implement each CNOT gate directly without expressing it as a combination of single qubit gates and CR gates, as it is done in IBM devices.We make this choice because in this way it is easier to solve numerically the target Lindblad equation.At each time step of the evolution we simulate a circuit with an increasing number of CNOT gates and measurements at the end.Thus we add SPAM channels (see appendices C and D) to model measurements errors.This allows to extend the analysis to runs on real hardware, that involve measurements, as we will show later in subsection VII B. In the three upper panels of Fig. 3 we show the time evolution of the ρ 22 = ⟨10|ρ|10⟩ entry of the density matrix.The relevant time interval is again given by a total duration of N ∼ 100 gates.Indeed the depolarizing error in this case spoils the quantum state after ∼ 120 CNOT gates.Fig. 3 (d) shows the mean of the Hellinger distances Hng σ (in blue) and Hibm σ (in red) and their standard deviations ∆H ng σ and ∆H ibm σ , also shown in the inset.Once more, within the relevant time interval [0, T ] the Hellinger distances obtained with the noisy gates simulations are closer to zero than those obtained with the Qiskit simulator and the difference between Hng σ and Hibm σ is of the order ∼ 10 −2 .This corresponds to a relative improvement in the range from 80% to 55% as time increases.This is shown in Fig. 3 (e).In the interesting regime [0, T ] the relative improvement is always above 55%.
By looking at Figs. 2 (e) and 3 (d), we notice that the improvement in the Hellinger distance gained by using the noisy gates approach is much higher for CR gates with respect to CNOT gates.The reason why this happens is clarified in Fig. 3  and red curves are obtained with Qiskit simulations.The noisy gates simulations make good predictions for both gate sequences, as the blue curves follow closely the orange curves.On the other hand, the Qiskit simulation for the CR gates is visibly different from the numerical solution of the Lindblad equation.This might be due to the fact that the CR gate is a block diagonal matrix with X(θ) in the upper block and X(−θ) in the lower block while the CNOT gate is block diagonal with an identity in the upper block and X(θ) in the lower block.The identity in the CNOT might lead to a lower influence of noises on the ρ 00 and ρ 11 entries of the density matrix.These observations explain why the Hellinger distances obtained with the noisy gates in different simulations are very good and similar to each other, while the Hellinger distance obtained with Qiskit is better for the CNOT with respect to the CR.Nevertheless, the noisy gates approach always outperforms the standard one by a significant amount, as shown by the relative improvements.

B. Comparison with the behaviour of a real quantum computer
Now, we inspect the performances of the noisy gates approach when trying to reproduce the behaviour of a real quantum computation.To this purpose, we first ex-tend the analysis of the CNOT gates sequence in subsection VII A, and then we focus on the inverse Quantum Fourier Transform (QFT † ).When dealing with a real hardware, we must take into account that the noise model we implement in this analysis (see section II) might not be accurate enough in describing the device, and that different quantum devices might behave very differently from one another.As we will show, despite the choice of a simple noise model and the instability of ibmq devices, our approach is still able to outperform the standard one also when compared with the real hardware.
CNOT simulations.We run the sequence of CNOT gates of subsection VII A on ibmq quito, available on the cloud and comprising 7 superconducting transmon qubits [47] (see appendix E) and we reconstruct the density matrix χ obtained from the physical device, to be compared with the density matrices ρ ng , ρ ibm and σ obtained for the CNOT simulations discussed in subsection VII A. We remark again that for the CNOT simulations of subsection VII A, we implemented each CNOT gate directly without expressing it as a combination of single qubit gates and CR gates, as it is done in IBM devices, because in this way it is easier to solve numerically the target Lindblad equation.We create a list of circuits, each consisting of an increasing number of CNOT gates, and measure each circuit 1000 times to obtain the output probability distributions, thus deriving the evolution of the outcome probabilities as the number of gates increases.As noted above, since each circuit involves measurements we added a SPAM error to model measurement errors.
The Hellinger distance H χ σ = H(χ, σ) between the Lindblad evolution and the evolution obtained with ibmq quito is shown in Fig. 4 (a).This distance is three to tens time larger with respect to Hng σ and Hibm σ that are shown in Fig. 3 (d).While the standard approach and the noisy gates approach have a certain level of agreement with the Lindblad equation, the latter is deviating from the quantum hardware by a significantly higher level.This is also the reason why it is not possible to appreciate the difference between the mean Hellinger distance Hng χ = H(ρ ng , χ) of the noisy gates with ibmq quito and the mean Hellinger distance Hibm χ = H(ρ ibm , χ) of Qiskit with ibmq quito, as shown in Fig. 4 (b) and Fig. 4 (c).Fig. 4 (d) shows the relative improvement with respect to the device, which is calculated as | Hibm χ − Hng χ |/ Hibm χ .The relative improvement is around 10%.The smaller relative improvement with respect to those shown in the previous subsection is only to a small extend due to the fact that we do not decompose CNOT gates.The main reason, as we explain when discussing the simulations of the QFT (see below), is that additional noises are present in ibmq devices, i.e. crosstalks, correlated noises and coherent errors [48,49].The simple noise model that we consider in this work does not take such noises into account.
QFT simulations The (QFT † ) is a subroutine of many important quantum algorithms, as for example the Shor's algorithm [50,51].An important feature of QFT † is that the circuit for n qubits is readily extendable to n + 1 qubits; thus we can efficiently test the robustness of the method as the circuit's width and depth increase.We run QFT † for n = 2, . . ., 5 on ibmq oslo and for n = 2, . . ., 18 on ibmq kolkata.These devices are available on the cloud, comprising respectively 7 and 27 superconducting transmon qubits [47], see appendix E for further details.We set as input of QFT † the state |+⟩ ⊗n , obtained by applying a layer of Hadamard gates on each qubit initialized in |0⟩.In this way the ideal output of QFT † should be |0⟩ ⊗n .Runs on real quantum computers are performed by taking 1000 shots, i.e. measurements.We also run the corresponding noisy gates and Qiskit simulations.(In appendix G we perform a similar analysis for the GHZ algorithm [52].) Implementing QFT † circuit on ibmq devices requires to transpile the circuit into their native gate set.We have defined a custom noise model in Qiskit, by adding after each gate of the transpiled circuit the depolarizing and relaxation channels, and the SPAM channel before measurements, see appendix C. Similarly, in the noisy gates simulation each gate is replaced with its noisy version according to the noise model in section II.During idle-times of qubits we put the relaxation noise gates (see appendix D) in order to take into account the stand-by times of the physical qubits; before measurements, we apply SPAM noise gates (see appendix D) which accounts for read-out errors.In these simulations the CNOT gates inside the circuits are decomposed in terms of single-qubit and CR gates, as in ibmq devices.
In order to measure the performance of different approaches in simulating the behaviour of the quantum computer, we look at their distance with the outcomes of the real device; this is achieved by computing the Hellinger distance between the probability distributions, and it can be done without performing full tomography on the quantum states, which scales exponentially with the number of qubits and becomes unfeasible for the cur-rent simulations.
In Fig. 5 (a) we plot the average values of H ng χ = H(ρ ng , χ), H ibm χ = H(ρ ibm , χ) as the number of qubits n increases from 2 to 5, where now the diagonal elements of χ are the outcome probabilities of ibmq oslo.As in the previous section, we have run 100 independent simulations, each including 1000 samples, for both methods and for each n, in order to compute the standard deviations ∆H ng χ , ∆H ibm χ shown in the insets of Fig. 5.Only for Fig. 5 (c) we have run a single simulation of 1000 samples, thus standard deviations are not present.We notice that for every n we get Hng χ < Hibm χ and ∆H ng χ < ∆H ibm χ .The relative improvement, shown in green in the insets of of Fig. 5, changes significantly between different devices and also for the same device but in different moments, namely with different noise parameters, meaning that the performances of such devices are not very stable.For example at n = 3, in the left panel the relative improvement is ∼ 25%, in the central panel it is ∼ 5% and in the right panel it is ∼ 25%.The highest relative improvement obtained with the run on ibmq oslo is ∼ 30% and for runs on ibmq kolkata is ∼ 35%.
The results show that our method is more accurate than existing ones.Actually, it reproduces the Lindblad dynamics better (Figures 1 and 2) than the dynamics of the quantum devices.The reason, mentioned before, is that quantum devices are affected by additional and more complicated noises, which are not taken into account by the noise model we are using; we stress again that to find a better noise model is not the scope of this work, and will be subject of future research.
The simulations on ibmq kolkata have been extended to 18 qubits to test the computational scalability of the from n = 9 to n = 18 qubits for the QFT † executed on ibmq kolkata.Since for n ≥ 9 the depth of the circuit is such that noises make the resulting probability distribution very flat, 1000 runs of the circuit on the quantum device, which returns a single computational basis state in each run, are not sufficient to reconstruct faithfully the probability distribution over the 2 n basis states; by increasing the number of qubits, a 0 probability is associated to an increasing number of basis states.The simulator instead does does not suffer from this limitation: in each run, it returns a non-zero value for each possible output.Then the respective probability distributions differ more and more, and this is the reason why the Hellinger distance rapidly increases.Despite the fact that the number of runs of the device is not sufficient to derive clear conclusions, we notice that the noisy gates approach still performs better than the standard one.noisy gates simulator.In Fig. 6 we show the mean Hellinger distances Hng χ and Hibm χ and their standard deviations ∆H ng χ and ∆H ibm χ from n = 9 to n = 18 qubits: simulations apparently become rapidly bad, since the Hellinger distance approaches 1, its maximum value.There is a clear reason behind that, which does not represent a limitation of our simulator.First of all, for such an high number of qubits, the depth of the transpiled circuit is so large that noises dominate [53] and the resulting probability distributions are very flat.Then, to recover a faithful probability distribution over the 2 n possible outcomes by the quantum device, which returns a single outcome in each run, the number of circuit runs must be significantly larger than 2 n .Therefore, for n ≥ 9 a number of runs equal to 1000 is not sufficient (and increasing this number becomes soon impractical): the output distribution from the device is increasingly dominated by 0's, while our simulator returns (in general) a non-zero probability for each output state: this makes the Hellinger distances of Fig. 6 approach 1.Nevertheless, also in this case the noisy gates approach performs better than the standard one, even if the number of runs of the quantum device are not enough to properly recover the full probability distribution.
As a final remark, we stress that we obtain better results with respect to Qiskit, despite the fact that we have chosen the simplest time dependent pulse shape in the Hamiltonians (see Eq. ( 17) and Eq. ( 31)).

VIII. CONCLUSIONS AND OUTLOOK
We have developed a novel approach, called noisy quantum gates, to improve classical simulations of NISQ computers: it is based on integrating the noise into the gates, rather than keeping gates and noise as two separate dynamics.We have shown that our approach is very successful in simulating the Lindblad dynamics, with a relative improvement between 50% and 90% and more, compared with the standard gate-noise separation method.
When compared against real quantum devices, the improvement fluctuates between 10% and 30%; this is largely due to the fact that the underlying noise model is too simple to accurately represent the dynamics of the device, as discussed in connection to the simulation of the CNOT gate.This is not a weakness of the noisy gate approach here presented, but of the underlying noise model, which we used since it is rather standard in the literature.
There is a number of potential improvements that can be straightforwardly implemented; all of them require an update of the noise model, not of the simulation strategy, which is already very good.First of all, there are likely additional single-qubit errors which should be taken into account, for example those induced by the driving pulses.Secondly, in the present work we considered only noncorrelated single-qubit errors, but the method can easily accommodate also correlated two-qubits errors [48,54] by introducing proper correlated noises into the stochastic equations.Another possible extension of the approach is to add in the Hamiltonians small interactions between adjacent qubits in order to mimic cross talk errors [49,55].Last, the current version of the noisy gates approach relies on the Lindblad equation that works in the Markovian limit; this is reflected in the fact that we used stochastic equations based on white noises.The approach can be generalized to non-Markovian dynamics by using colored noises, as already discussed in the literature in different contexts [19][20][21][22].
Furthermore, our approach is also useful for other purposes that go beyond plane error analysis.For example, the shape of the pulse in the driving Hamiltonians, (see Eq. ( 17) and Eq. ( 31)), can affect the noise.In our work we chose for simplicity a rectangular shape, but usually in real devices different shapes can be used, for example Gaussian ones.Consequentially, a natural application of our approach is error mitigation [56,57], by optimizing the parameters of the pulse in order to minimize the effect of the noise [58][59][60]; the optimization can be performed for example by exploiting machine learning techniques, to find the best pulse parameters, which can be tested on real quantum hardware.
In this work we specified our approach to the native gate set and noise model of IBM devices; clearly the approach is general and can be used to describe in principle any NISQ platform.get a system of SDEs: which must be solved with the initial conditions ψ 0 0 = |ψ 0 ⟩.The zeroth order differential equation is the deterministic one given by the Hamiltonian evolution alone, hence its solution is simply ψ 0 s = U s |ψ 0 ⟩.The first order SDE is an example of a time-dependent Ornstein-Uhlenbeck process [37]: the solution is where we defined S s = s 0 dW τ L τ .Finally, the solution to the second order SDE is where L s = U † s LU s .Then, the solution at order ϵ 2 is given by |ψ 1 ⟩ = N |ψ 0 ⟩ + O(ϵ 3 ), where the evolution operator is N = U g N ′ , with In order to evaluate the solution in the form given in the main text, we make use of the following equality: s (A16) obtained by using the Itô rule [37] for each entry of the stochastic matrices.Substituting this expression into Eq.(A15),we get to second order: where Λ, Ξ and C = 1 0 dW s L s , S s are the same quantities defined in the main text.

Appendix B: Comparison of the approximations
We focus on the main differences between the standard approximation made in error analysis against the one considered in the noisy gates approach.
Given the following Lindblad master equation where L [ρ t ] = Lρ t L † − 1 2 {L † L, ρ t }, let's move to the interaction picture by defining χ t = U † t,t0 ρ t U t,t0 and χ t0 = ρ t0 .Then where The formal solution of Eq. (B2) is where T [•] is the time ordering.Thus in the Schrödinger picture we can write the formal solution of Eq. (B1) as -Standard approximation -The main approximation that can be found in the literature is to separate the Hamiltonian dynamics from the noise one [15,16].This choice is based on the observation that in general in quantum devices ω >> γ, where ω is the pulse frequency of the Hamiltonian.Thus, the noise dynamics can be seen as frozen with respect to the faster Hamiltonian one.It means that in Eq. (B4) one assumes We notice that indeed in Eq. (B6) the two dynamics are independent.
-Noisy gates approximation -Also in this case the approximation is based on ω >> γ, but we assume that γ is not small enough to completely separate the dynamics.An example of this can be seen in the devices of IBM where the noise evolution can be influenced in a nonnegligible manner by the pulse of the drive Hamiltonian [46,61].Thus we make a first order approximation over γ in Eq. (B4) T e and we get In Eq. (B8) the noise depends on the Hamiltonian dynamics through L(s).We stress that the perturbative solution of the SDE in the noisy gates model reproduce density matrices of the form of Eq. (B8).
and the stochastic term of the relative Itô equation reads: (D5) With this stochastic term the Itô equation is analytically solvable [40] and we get the following non-unitary noisy gate N relax (t, t 0 ) = e iα W2(t,t0) iS(t, t 0 )e iα W2(t,t0) 0 e − γ 1 2 (t−t0) e −iα W2(t,t0) , (D6) where we defined for simplicity α := γ pd /4, and is a complex stochastic Itô process.In principle, such a term is problematic in view of a simulation, since it is not easy to sample.To understand this, look for instance at the real part, S R (t, t 0 ) = √ γ 1 t t0 e − γ 1 2 (s−t0) cos 2α W2 (s, t 0 ) dW s,1 ; (D8) this is an Itô integral of a stochastic function, and it is not easy to derive its probability distribution; thus, sampling S(t, t 0 ) may be problematic.We can avoid such a difficulty by adequately substituting N relax (t, t 0 ) with some modified noisy gate, which is equivalent to the former once the average is carried out, in the sense that Eq. (C3) still holds even if the new noisy gate is not a solution of the unraveling (D5) anymore.For instance, it is straightforward to verify that this holds for the following choice: Ñ relax (t, t 0 ) = e iα W2(t,t0) i S(t, t 0 )e −iα W2(t,t0) 0 e − γ 1 2 (t−t0) e −iα W2(t,t0) , (D9) with the definition S(t, t 0 ) = √ γ 1 (D11) The difference is that now the process S(t, t 0 ) is just the Itô integral of a deterministic function, hence we know that it must have a Gaussian statistics [37], which makes it more convenient for a simulation.

Appendix E: Device parameters
For the simulations in Sec.VII and in appendix G we used the device parameters provided by IBM.Here we report the average value of such parameters.ibmq manila contains 5 fixed-frequency transmons qubits [47], with median fundamental transition frequency of 4.962 GHz and median anharmonicity of −0.34358 GHz.The median qubit lifetime T 1 of the qubits is 149.11µs, the median coherence time T 2 is 44.43µs and the median readout error is 0.0217.The single qubit gate error varies between 1.975 × 10 −4 and 6.138 × 10 −4 , while the CNOT error varies between 7.072 × 10 −3 and 1.125 × 10 −2 , depending on the specific connection.In the simulations that reproduce the Lindblad equations, parameters of qubits zero and one were used.ibmq kolkata contains 27 fixedfrequency transmons qubits, with median fundamental transition frequency of 5.102 GHz and median anharmonicity of −0.34345 GHz.The median qubit lifetime T 1 of the qubits is 127.39µs, the median coherence time T 2 is 86.41µs and the median readout error is 0.0132.The single qubit gate error varies between 1.443 × 10 −4 and 5.410 × 10 −3 , while the CNOT error varies between 4.214 × 10 −3 and 1 × 10 −2 , depending on the specific connection.The qubits, which are used to run QFT † algorithm, belong to the list [0, 1,4,7,10,12,15,18,21,23,24,25,22,19,16,14,11,8,5,3,2].ibmq quito contains 7 fixed-frequency transmons qubits, with the median fundamental transition frequency of 5.164 GHz and median anharmonicity of −0.3315 GHz.The median qubit lifetime T 1 of the qubits is 105.84µs, the median coherence time T 2 is 84.05µs and the median readout error is 0.044.The single qubit gate error varies between 3.054 × 10 −4 and 6.929 × 10 −4 , while the CNOT error varies between 9.682 × 10 −3 and 1.463 × 10 −2 , depending on the specific connection.The qubits, which are used for CNOT gate sequence are 0 and 1. ibmq oslo contains 7 fixed-frequency transmons qubits, with the median fundamental transition frequency of 5.046 GHz and median anharmonicity of −0.3429 GHz.The median qubit lifetime T 1 of the qubits is 128.12µs, the median coherence time T 2 is 58.57µs and the median readout error is 0.0216.The single qubit gate error varies between 1.648 × 10 −4 and 6.698 × 10 −4 , while the CNOT error varies between 6.471 × 10 −3 and 2.067 × 10 −2 , depending on the specific connection.The qubits, which are used to run GHZ algorithm, belong to the list [0, 1,3,5,4].FIG. 9. On panel (a), probabilities histograms for 4 qubits of a single independent simulation of the GHZ algorithm.In orange the results for ibmq oslo, in blue for the noisy gates and in red for the Qiskit simulator.On panel (b), Hellinger distance for the GHZ algorithm for n = 2, . . ., 5 qubits.Each value is the mean of 100 independent simulations for the noisy gates, in blue, and for the Qiskit simulations, in red.The left inset shows the relative improvement, calculated as | Hibm σ − Hng σ |/ Hσibm, while the right inset shows the standard deviations as functions of the number of qubits.

FIG. 1 .
FIG. 1. Repetition of X gates.The three upper panels (a), (b), (c) show the time evolution of the ρ00 = ⟨0|ρ|0⟩ entry of the density matrix.The numerical solution of the Lindblad equation is displayed in orange (a), that of the noisy gates simulation in blue (b), and that of the Qiskit simulation in red (c).The noisy gates and Qiskit simulations are obtained with 1000 samples, and qualitatively they reproduce the time evolution of the Lindblad equation.Vertical dashed lines in the three top panels represent the time scales of relaxation T1 (green), T2 (yellow) and depolarization T d (grey).Panel (d) shows the Hellinger distances H ng σ , in blue, and H ibm σ , in red, as a function of time.Different curves are obtained from 100 independent runs of the two methods (for better readability only five are shown), where each simulation is obtained by averaging over 1000 samples.Panel (e) shows the mean of the Hellinger distances Hng σ , and Hibm σ , obtained from the 100 independent runs, and vertical error bars show their standard deviations ∆H ng σ , ∆H ibm σ .The inset displays ∆H ng σ and ∆H ibm σ as functions of time.Panel (f) shows the relative improvement of the distance Hng σ with respect to Hibm σ , calculated as | Hibm σ − Hng σ |/ Hibm σ .The fact that noises drive the system towards the maximally mixed state is the reason why the improvement decreases in time.The noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times, as the noise is dominant over the unitary evolution.In the interesting regime [0, 2000 • tg] before decoherence dominates, our improvement is always above 60%.

FIG. 2 .
FIG. 2. Repetition of CR gates.The three upper panels (a), (b), (c) show the time evolution of the ρ22 entry of the density matrix for the CR gate with θ = π and ϕ = 0. Colors have the same meaning as for Fig. 1.Vertical dashed lines represent the time scales of relaxation, T1 (in green) and T2 (in yellow) of the target qubit, and depolarization T d (grey).The noisy gates simulations reproduce qualitatively better the time evolution obtained from the direct numerical solution of the Lindblad equation.Panels (d) and (e) display the Hellinger distances H ng σ , in blue, and H ibm σ , in red, as a function of time, for a repetition of CR gates.The plots have the same meaning as for Fig.1.Panel (f) shows the relative improvement of the distance Hng σ with respect to Hibm σ , calculated as | Hibm σ − Hng σ |/ Hibm σ .The fact that noises drive the system towards the maximally mixed state is the reason why the improvement decreases in time.The noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times, as the noise is dominant over the unitary evolution.In the interesting regime [0, 100 • tg] our improvement is always above 88%.

FIG. 3 .
FIG. 3. Repetition of CNOT gates.The three upper panels (a), (b), (c) show the time evolution of the ρ22 entry of the density matrix for the CNOT gate.Colors have the same meaning as for Fig. 1.The noisy gates simulations reproduce qualitatively better the time evolution obtained from the direct numerical solution of the Lindblad equation.Panel (d) displays mean of the Hellinger distances Hngσ , in blue, and Hibm σ , in red, and their standard deviations as functions of time.Panel (e) shows the relative improvement.The fact that noises drive the system towards the maximally mixed state is again the reason why the improvement decreases in time: the noisy gates and the standard approaches lead to the same predictions when one is close to decoherence times, as the noise is dominant over the unitary evolution.In the interesting regime [0, 100 • tg] our improvement is always above 55%.The upper subplots of panel (f) show the time evolution of the ρ11 entry of the density matrix for the same sequence of CR gates in Fig.2and the lower subplots show the time evolution of the ρ11 entry of the density matrix for the sequence of CNOT gates.Colors have the same meaning as for Fig.1.For the CR gates, the Qiskit simulation of ρ11 is visibly different from the Lindblad evolution, thus explaining the higher improvement of the noisy gates simulation in the Hellinger distance in Fig.2.

FIG. 4 .
FIG. 4. Repetition of CNOT gates.Panel (a) shows the Hellinger distance H χ σ between the Lindblad evolution and ibmq quito for the repetition of CNOT gates.Panel (b) shows the mean Hellinger distance Hng χ and the standard deviations between the noisy gates simulation and ibmq quito.Panel (c) shows the mean Hellinger distance Hibm χ and the standard deviations between the Qiskit simulation and ibmq quito.Panel (d) shows the relative improvement calculated as | Hibm χ − Hng χ |/ Hibm χ .The relative improvement is around 10%.The smaller relative improvement with respect to those shown in the previous figures, is mainly due to additional noises present in ibmq devices, i.e. crosstalks, correlated noises and coherent errors.

FIG. 5 .
FIG. 5.Quantum Fourier Transform.Panel (a)  shows the Hellinger distances between the noisy gate approach and ibmq oslo, and between the Qiskit simulator and ibmq oslo, when executing the QFT † algorithm for n = 2, . . ., 5 qubits.Each value is the mean of 100 independent simulations for the noisy gates, in blue, and for the Qiskit simulations, in red.The left inset shows the relative improvement, calculated as | Hibm χ − Hng χ |/ Hibm χ , while the right inset shows the standard deviations as functions of the number of qubits.Panel (b) is the same as (a), with ibmq oslo replaced by ibmq kolkata and the number of qubit going up to 8. In panel (c) the comparison presented in (b) has been repeated a second time on ibmq kolkata; in this case, only a single simulation of 1000 samples is considered.The inset shows the relative improvement.
Fig. 5 (b) displays again the average values of H ng χ = H(ρ ng , χ), H ibm χ = H(ρ ibm , χ) up to 8 qubits, where now the diagonal elements of χ come from ibmq kolkata.As shown in Fig. 5 (c) we compute again H ng χ , H ibm χ to test the stability of ibmq kolkata in different runs.

FIG. 7 .
FIG. 7. Fidelities F ng σ , in blue, and F ibm σ , in red, as a function of time, for a repetition of X gates.On panel (a), the fidelities obtained from 100 independent runs of the two methods are pictured (for better readability only five are shown), where each simulation is obtained by averaging over 1000 samples.On panel (b), the means Fng σ , Fibm σ of the same simulations and their standard deviations ∆F ng σ , ∆F ibm σ are displayed.The inset shows the standard deviations ∆F ng σ , ∆F ibm σ as functions of time.

FIG. 8 .
FIG. 8. Fidelities F ng σ , in blue, and F ibm σ , in red, as a function of time, for a repetition of CR gates.Panels (a) and (b) have the same meaning as for Fig.7.