# Efficient Simulation of Parametrized Quantum Circuits under Nonunital Noise through Pauli Backpropagation

Victor Martinez, Armando Angrisani, Ekaterina Pankovets, Omar Fawzi, Daniel Stilck França

► To cite this version:

Victor Martinez, Armando Angrisani, Ekaterina Pankovets, Omar Fawzi, Daniel Stilck França. Efficient Simulation of Parametrized Quantum Circuits under Nonunital Noise through Pauli Backpropagation. Physical Review Letters, 2025, 134 (25), pp.250602. 10.1103/j1gg-s6zb. hal-05218980

HAL Id: hal-05218980

https://hal.science/hal-05218980v1

Submitted on 30 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient simulation of parametrized quantum circuits under non-unital noise through Pauli backpropagation

Victor Martinez*

IBM France, Avenue de l'Europe, 92275 Bois-Colombes, France and
Inria, ENS Lyon, UCBL, LIP, F-69342 Lyon Cedex 07, France

Armando Angrisani† and Ekaterina Pankovets‡

Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

Omar Fawzi§

Inria, ENS Lyon, UCBL, LIP, F-69342 Lyon Cedex 07, France

Daniel Stilck França¶

Inria, ENS Lyon, UCBL, LIP, F-69342 Lyon Cedex 07, France and

Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Denmark

As quantum devices continue to grow in size but remain affected by noise, it is crucial to determine when and how they can outperform classical computers on practical tasks. A central piece in this effort is to develop the most efficient classical simulation algorithms possible. Among the most promising approaches are Pauli backpropagation algorithms, which have already demonstrated their ability to efficiently simulate certain classes of parameterized quantum circuits—a leading contender for near-term quantum advantage—under random circuit assumptions and depolarizing noise. However, their efficiency was not previously established for more realistic non-unital noise models, such as amplitude damping, that better capture noise on existing hardware. Here, we close this gap by adapting Pauli backpropagation to non-unital noise, proving that it remains efficient even under these more challenging conditions. Our proof leverages a refined combinatorial analysis to handle the complexities introduced by non-unital channels, thus strengthening Pauli backpropagation as a powerful tool for simulating near-term quantum devices.

## I. INTRODUCTION

With the recent substantial growth in both the size and quality of quantum hardware, the quantum community has searched for quantum algorithms that can be implemented on near-term hardware and provide an advantage over classical algorithms on practically relevant problems [1–4].

In this search for efficient and practical quantum computing, parameterized quantum circuits (PQC) have gained considerable attention [5, 6]. These circuits play a crucial role in variational quantum algorithms (VQA) and are frequently employed in areas such as optimization, with algorithms like QAOA [7] and VQE [8, 9], as well as in machine learning [4, 10]. These algorithms work by optimising the parameters, often rotation angles, of the quantum circuit to minimise or maximise an objective function, made of expectation values of observables. Initially thought to be somewhat noise resistant, strong limitations on their performance in the pres-

ence of noise were later shown [11], essentially showing that the outputs of these algorithms offer no quantum advantage at constant [12] depth or moderately above that [13] and other works established that noisy quantum computers can be simulated in various regimes. However, until recently [14, 15], the vast majority of such results was restricted to noise that is *unital*, i.e. maps the maximally mixed state to itself. However, the noise in current devices often does not satisfy this assumption [16, 17]. Thus, it is imperative to understand to what extent we can simulate noisy quantum computation under non-unital noise. It is important to note that extending classical algorithms to such noise is far from being straightforward. Indeed, the effect of non-unital noise is known to be much more subtle than unital noise and strong qualitative gaps exist between the two. For instance, as shown in [18], without access to fresh qubits we can only implement log depth circuits under depolarizing noise, whereas we can compute for exponential times with non-unital noise.

In this work, we will show how to efficiently simulate noisy PQC s under non-unital noise using Pauli backpropagation methods. Among the various simulation approaches for VQAs, Pauli backpropagation methods are showing significant promise in practice and offer strong theoretical guarantees [14, 15, 19–28], showcasing great accuracy on systems with up to 127 qubits, even in noise-

<footnote>* Victor Martinez victor.martinez@ibm.com</footnote>

<footnote>† Armando Angrisani armando.angrisani@epfl.ch</footnote>

<footnote>‡ Ekaterina Pankovets ekaterina.pankovets@epfl.ch</footnote>

<footnote>§ Omar Fawzi omar.fawzi@ens-lyon.fr</footnote>

<footnote>¶ Daniel Stilck França dsfranca@math.ku.dk</footnote>
less conditions [20, 22, 23, 29]. This approach focuses on structured circuits composed of alternating layers of parameterized rotations and Clifford operations. These circuits, which form a universal family, closely resemble the ones commonly used in PQCs. Using the circuit structure and the effect of unital noise, it has been shown that these simulation methods can recover the energy landscape of the PQC on average over the parameters in polynomial time [19, 23].

Among the Pauli backpropagation methods, LOWESA [19, 20, 24, 25] has demonstrated several notable advantages. The algorithm is both memory-efficient and parallelizable. Additionally, the most resource-intensive part of the algorithm, Pauli backpropagation, only depends on the structure of the circuit and not on the values of the angles. Thus, the algorithm can be seen as a method to efficiently estimate the energy landscape as a function of the angles. Finally, the runtime of the algorithm is independent of both the quantum chip geometry and the depth of the quantum circuit.

From a technical viewpoint, the proof techniques we use require dealing with combinatorial arguments in a more subtle way than in [19], along with new Monte-Carlo techniques to handle the non-unital component of the noise. The main reason for that is that depolarizing noise uniformly dampens all possible Pauli strings, whereas non-unital noise does not. Nevertheless, our final result has a similar scaling as the ones of [19].

Another extension of Pauli backpropagation to random circuits with arbitrary uniform noise, possibly non-unital or dephasing, was also provided in the related paper [30], which leverages the deterministic low-weight truncation scheme previously employed in Refs. [14, 27, 28].

## II. SETUP

As in [19, 20, 24], the parametrized quantum circuits we consider are made of $m$ alternating layers of single-qubit rotations around the Z-axis and Clifford operations $C_i$. The rotations are parametrized by angles $\theta = (\theta_1, \dots, \theta_m)$ and qubit labels $(q_1, \dots, q_m)$ so that the single-qubit gate rotation at layer $i$ is given by $R_z^{(q_i)}(\theta_i) = e^{-j\theta_i/2Z_{qi}}$. The circuits can be written in the following way, and are depicted in Figure 1.

$$U(\theta) = \left( \prod_{i=1}^{m} C_i R_z^{(q_i)}(\theta_i) \right) C_0. \quad (1)$$

Note that these circuits form a universal family of quantum circuits, and that the complexity of classically simulating such quantum circuits is exponential in the number of $R_z$ rotations $m$ [31]. In the following discussions, we will focus only on the gate noise affecting the

<figure><img src="image_2.png" /><figcaption>FIG. 1: Representation of parametrised quantum circuits simulated by LOWESA, made of alternating layers of Clifford operations, amplitude-damping noise and single-qubit rotations</figcaption></figure>

rotations, as it is often the case that noise on the Clifford operations can be more easily managed since they can be implemented transversally [32]. For simplicity, the noise channel we consider is the amplitude damping channel $N_{AD}$. However, we show in the supplemental material that our results hold for most single-qubit noise models considered in the literature, including the composition and probabilistic combinations of amplitude damping, dephasing and depolarizing. Our noisy quantum circuit is thus represented by the channel

$$U_{\theta} = \left( \bigcirc_{i=1}^{m} C_i \circ \mathcal{R}_z^{(q_i)}(\theta_i) \circ N_{AD} \right) \circ C_0, \quad (2)$$

where we always denote in caligraphic font the channel given by conjugation by the corresponding unitary.

In most VQAs, the algorithm requires applying the quantum circuit to the initial state $|0\rangle^{\otimes n}$ and then measuring an observable [7–9]. Our classical algorithm will estimate the expectation value of this observable. Denote by $\mathbb{P}$ the set of Pauli matrices $\mathbb{P} = \{I, X, Y, Z\}$ and by $\mathcal{P}$ a Pauli string in $\mathbb{P}^{\otimes n}$, which will be the observable we consider. In practice, if the observable is not a Pauli string, one can always still decompose it as a sum of Pauli strings and then run the simulation procedure for each of the Pauli strings, as long as there are polynomially many of them. Our quantity of interest is the expectation value as a function of the parameters $\theta \in [0, 2\pi]^m$:

$$f(\theta) = \mathrm{tr}(U_{\theta}(|0\rangle\langle0|)P). \quad (3)$$

We will approximate this expectation value by $\tilde{f}(\theta)$ that can be computed in polynomial time. The approximation error we consider is the $L^2$-norm over the parameter space $\Theta = [0, 2\pi]^m$:

$$\begin{aligned} \Delta(f, \tilde{f}) &= ||f - \tilde{f}||_2 \\ &= \left( \frac{1}{|\Theta|} \int_{\Theta} |f(\theta) - \tilde{f}(\theta)|^2 d\theta \right)^{1/2}. \end{aligned} \quad (4)$$

Note that convergence in this metric does not imply that our approximation scheme succeeds for a given set of angles.
When backpropagating observables in the Heisenberg picture, it is often useful to use the Pauli Transfer Matrix (PTM) formalism [33], which we briefly introduce here. In the PTM, operators and states are represented in the Pauli basis, such that the PTM representation of a quantum state $\rho$ is the real valued vector $|\rho\rangle \in \mathbb{R}^{2^n}$ where $[|\rho\rangle\rangle]_i = \text{tr}(\rho P_i)$ for $P_i \in \mathbb{P}^{\otimes n}$. Similarly, a quantum channel on $n$-qubits $\mathcal{E}$ is represented by the matrix $\mathbf{E}$ with $\mathbf{E}_{ij} = 2^{-n}\text{tr}(P_i\mathcal{E}(P_j))$. Using this formalism, it is clear that the expectation value of a Pauli string can be written as $\text{tr}(\mathcal{P}\mathcal{E}(\rho)) = \langle\langle\mathbf{P}|\mathbf{E}|\rho\rangle\rangle$. Computing expectation values in the Heisenberg picture, where quantum channels act on measurement operators instead of quantum states, can be done using this formalism.

Indeed, the corresponding adjoint operation simply requires taking the transpose of the expectation value, $\langle\langle\mathbf{P}|\mathbf{E}|\rho\rangle\rangle = \langle\langle\rho|\mathbf{E}^T|\mathbf{P}\rangle\rangle$. If $\mathcal{E}$ is made of multiple 2-qubit Clifford unitaries, this gives an efficient way of computing expectation values. Indeed, in the PTM, each Clifford unitary takes the Pauli operator $\mathbf{P}$ to another Pauli operator $\mathbf{P}'$, up to a phase. For 2-qubit unitaries keeping track of this change takes at most $\mathcal{O}(n^2)$ [34]. As long as there are poly($n$) 2-qubit Clifford unitaries, this gives us an efficient way of computing expectation values of Clifford circuits in the Heisenberg picture. Handling the rotations, which are not Clifford gates, requires more work and the use of the structure of the noise channel.

### III. MAIN RESULT

To deal with the rotations, we start by noting that the adjoint of the amplitude damping channel with damping parameter $\gamma > 0$ can easily be defined in the PTM since,

$$ \mathcal{N}_{AD}^{\dagger}(X) = \sqrt{1-\gamma}X, \quad \mathcal{N}_{AD}^{\dagger}(Y) = \sqrt{1-\gamma}Y, \quad \mathcal{N}_{AD}^{\dagger}(Z) = (1-\gamma)Z + \gamma I, \quad \mathcal{N}_{AD}^{\dagger}(I) = I \qquad (5) $$

Let $\mathbf{R}_z(\theta)$ be the PTM representation of the adjoint of $\mathcal{R}_z(\theta)$ and $\mathbf{N}_{AD}$ be the PTM representation of the adjoint of the noise channel. The noisy rotations can be explicitly computed,

$$ \mathbf{R}_z \cdot \mathbf{N}_{AD} = \begin{bmatrix} 1 & 0 & 0 & \gamma \\ 0 & \sqrt{1-\gamma} \cos \theta & -\sqrt{1-\gamma} \sin \theta & 0 \\ 0 & \sqrt{1-\gamma} \sin \theta & \sqrt{1-\gamma} \cos \theta & 0 \\ 0 & 0 & 0 & 1-\gamma \end{bmatrix} \quad (6) \\ = \mathbf{D}_0 + \mathbf{D}_{0z} + \mathbf{D}_{0i} + \mathbf{D}_1 \cos \theta + \mathbf{D}_{-1} \sin \theta $$

where we have defined the quantum processes $\mathbf{D}_0 = |I\rangle\langle I|$, $\mathbf{D}_{0z} = (1-\gamma)|Z\rangle\langle Z|$, $\mathbf{D}_{0i} = \gamma|I\rangle\langle Z|$, $\mathbf{D}_1 = \sqrt{1-\gamma}(|X\rangle\langle X| + |Y\rangle\langle Y|)$ and $\mathbf{D}_{-1} = \sqrt{1-\gamma}(|Y\rangle\langle X| - |X\rangle\langle Y|)$, where $|I\rangle = [1 \ 0 \ 0 \ 0]^T$, $|X\rangle = [0 \ 1 \ 0 \ 0]^T$ ...

Therefore, when backpropagating through the circuit in the Heisenberg picture, the noisy rotation on qubit $q_i$

<figure><img src="image_3.png" /><figcaption>FIG. 2: Tree corresponding to the backpropagation of a single Pauli observable through the circuit. Each node corresponds to a single Pauli string backpropagated through the quantum circuits through the quantum processes $\mathcal{D}_\omega$.</figcaption></figure>

acts in the following way. If $P_{q_i} = I$, nothing happens, and $P_{q_i}$ remains unchanged through $\mathbf{D}_0$. If $P_{q_i} = X$, then $P_{q_i}$ is mapped to $P'_{q_i} = \sqrt{1-\gamma} \cos \theta X$ through process $\mathbf{D}_1$ or to $P'_{q_i} = \sqrt{1-\gamma} \sin \theta Y$ through process $\mathbf{D}_{-1}$. Similarly, if $P_{q_i} = Y$, then it is mapped to $P'_{q_i} = \sqrt{1-\gamma} \cos \theta Y$ through process $\mathbf{D}_1$ or to $P'_{q_i} = -\sqrt{1-\gamma} \sin \theta X$ through process $\mathbf{D}_{-1}$. Finally, if $P_{q_i} = Z$, then it is either mapped to $P'_{q_i} = \gamma I$ through $\mathbf{D}_{0i}$ or to $P'_{q_i} = (1-\gamma)Z$ through $\mathbf{D}_{0z}$.

Except for $I$, we see that when going through the noisy rotation, every Pauli is split into two downscaled Paulis. This gives us a way of computing the expectation value of our noisy circuit. Start with our observable, which is a Pauli string $\mathbf{P}$, then backpropagate it through each layer backwards, alternating layers of Clifford unitaries and single-qubit rotations. The Clifford unitaries update the Pauli string into another Pauli string - up to a phase factor - which can be computed in time $\mathcal{O}(n^2)$. The single-qubit rotation $R_z^{(q_i)}$ takes the Pauli string $\mathbf{P}$ into potentially 2 paths, corresponding to different Pauli strings differing on qubit $q_i$. We then keep track of the new paths created and reiterate this procedure. This creates a tree structure $\mathbf{T}$, depicted in Figure 2, which may have up to $2^m$ branches, where $m$ is our number of parameterized rotations.

Similarly to the unital case [19], formally decomposing our circuit using this tree structure requires us to keep track at each layer of both the rotation angles and which process was encountered for each path. This is done by storing, for each path (i.e. each branch of the tree), the
processes that were encountered during the backpropa-
gation in a vector $\omega \in \{0, 0_z, 0_I, \pm 1\}^m$ and by storing
separately the rotation angles. These two quantities are
then combined using the convenient trigonometric mono-
mials $\phi_{\omega_i}(\theta_i)$, with $\phi_0(\theta_i) = \phi_{0_z}(\theta_i) = \phi_{0_I}(\theta_i) = 1$,
$\phi_1(\theta_i) = \cos\theta_i$, and $\phi_{-1}(\theta_i) = \sin\theta_i$. By defining the
quantity $\Phi_{\omega}(\theta) = \prod_{i=1}^m \phi_{\omega_i}(\theta_i)$, our quantum circuit can
be decomposed using the quantum processes $D_{\omega}$, which
we have previously defined in the PTM:

$$
\mathcal{U}_{\theta} = \sum_{\omega \in \{0,0_z,0_I,1,-1\}^m} \Phi_{\omega}(\theta) \mathcal{D}_{\omega}. \quad (7)
$$

Where each process $\mathcal{D}_\omega$ is given by $\mathcal{D}_\omega = (\bigcirc_{i=1} \mathcal{C}_i \circ \mathcal{D}_{\omega_i}) \circ \mathcal{C}_0$. In this setting, our parametrized expectation value $f$ can then be written as:

$$
f(\theta) = \mathrm{tr}(\mathcal{U}_{\theta}[|0\rangle\langle0|]P) = \sum_{\omega} d_{\omega} \Phi_{\omega}(\theta) \quad (8)
$$

where the Fourier coefficients $d_\omega$ are given by,

$$
d_{\omega} = \mathrm{tr}(\mathcal{D}_{\omega}[|0\rangle\langle0|]P) = \langle\langle P | \mathbf{D}_{\omega} | 0\rangle\rangle. \quad (9)
$$

Recall that when backpropagating a Pauli string $P$, the rotation at each layer acts only on a single qubit $q_i$. Depending on the Pauli operator at qubit $q_i$, at most two processes are valid: $\mathbf{D}_0$ for $I$, $\mathbf{D}_{0z}$ and $\mathbf{D}_{0I}$ for $Z$, and $\mathbf{D}_{\pm 1}$ for $X$ and $Y$.

To approximate this expectation value, we propose two
different methods, each with its own caveats and advan-
tages. The first method, Algorithm 1, considers only the
paths in the tree T that have split fewer than ℓ times,
where ℓ is a chosen cut-off parameter. Indeed, each split
of the paths dampens the corresponding Pauli string,
such that discarding paths that have split sufficiently
could yield a good approximation. Let #ω denote the
number of times a given path ω has split. Algorithm 1
approximates the function f by the quantity:

$$
\tilde{f}(\theta) = \sum_{\substack{\omega \in T \\ \#\omega \le \ell}} d_{\omega} \Phi_{\omega}(\theta) \quad (10)
$$

which is computed in polynomial time for $\ell = O(1)$.
However, unlike in the case of unital noise, it is uncertain
whether each split effectively dampens the expectation
value of the Pauli observable. Specifically, when a split
occurs due to encountering a Pauli Z at qubit $q_i$, it results
in two new branches through processes $\mathbf{D}_{0I}$ and $\mathbf{D}_{0Z}$. In
the case of amplitude damping alone, this split may leave
the approximation error unchanged. Therefore, it is es-
sential to distinguish between splits occurring through
processes $\mathbf{D}_{\pm 1}$, which always dampen the approximation
error, and those occurring through $\mathbf{D}_{0Z}$ and $\mathbf{D}_{0I}$. This

is done by noting which splits occur on discarded paths
when doing the backpropagation.

The second approach, Algorithm 2, gets rid of this issue entirely by randomly sampling only one of the processes $\mathbf{D}_{0z}/\mathbf{D}_{0I}$ in the tree. Any subtree $T_k$ obtained by this sampling procedure only contains splits that effectively dampen the expectation value of the observable. Each of the resulting trees are then truncated in the same manner after enough splitting has occurred. By repeating the randomized sampling process for $K$ independent trees $T_k$, our quantity of interest is approximated by the empirical average

$$
\tilde{f}(\theta) = \frac{1}{K} \sum_{k=1}^{K} \sum_{\substack{\omega \in T_k \\ \#\omega \le l}} \Phi_{\omega}(\theta) d_{\omega} \quad (11)
$$

We present the following theorem for the performance
of both algorithms,

**Theorem 1.** Let $\mathcal{U}_\theta$ be a noisy quantum circuit as defined in Equation 2 with amplitude damping parameter $\gamma > 0$ and $P$ a Pauli observable with expectation value $f(\theta) = \mathrm{tr}(\mathcal{U}_\theta[|0\rangle\langle0|]P)$. For a given cut-off parameter $\ell$, Algorithm 1 outputs in time $O(n^2m2^\ell)$ an approximation $\tilde{f}$ and a lower bound $r$ on the minimum number of times the processes $\mathbf{D}_{\pm 1}$ have been encountered on every discarded branch during the backpropagation, such that,

$$
\Delta(f, \tilde{f}) \le (1 - \gamma)^{r/2} \quad (12)
$$

Given the same cut-off parameter $\ell$, and a sampling
overhead $K$, Algorithm 2 outputs in time $O(Kn^2m2^\ell)$
an approximation $\tilde{f}$ that satisfies with probability at least
$1 - \delta$,

$$
\Delta(f, \tilde{f}) \le (1 - \gamma)^{(\ell+1)/2} + \sqrt{\frac{2 \log (\delta^{-1}/2)}{K}} \quad (13)
$$

where the probability comes from the random sampling of
the trees.

In particular, this means that any quantum circuit
of our universal family can be efficiently simulated un-
der non-unital noise in the L²-error metric. This result
greatly extends what was previously known for simula-
tion under non-unital noise, eliminating many of the as-
sumptions on random circuit/noise introduced in prior
works [14, 15]. Furthermore, Equation 12 (and B13 in
Appendix B) show that for most quantum circuits, cel-
ebrated simulation methods [19, 20, 24] are expected to
work also for non-unital noise. Indeed, although the for-
mal guarantees we give in the result above are generally
better for Alg. 2, since it may be that $r = 0$ for some
circuits and observables, Alg. 1 is closer to what is done
in practice and is deterministic. We show with numerical
evidence that in practice, one expects $r$ to scale linearly with $\ell$ (see Figure 5 in Appendix A). Equation 12 also gives a certificate on how good Algorithm 1 performs for a given circuit and observable $P$, based on the computed value of $r$. We give here a sketch of the proof and refer to Appendices B and C for more detailed ones. The complexity of the algorithms is straightforward to derive. Computing the effect of a Clifford unitary takes $O(n^2)$, and this needs to be done for all paths at every layer. Since there are at most $2^\ell$ paths and $m$ layers, the total complexity is at most $O(n^2m2^\ell)$. For Algorithm 2, the Monte-Carlo simulation requires doing this process $K$ times, such that the runtime is $O(Kn^2m2^\ell)$. Let's now focus on the approximation error made by the algorithm. The $L^2$-error made by truncating a tree can be rewritten (see Lemma 2):

$$ 
\begin{aligned}
\Delta^2(\tilde{f}, f) &= \frac{1}{|\Theta|} \int_{\Theta} |\tilde{f}(\theta) - f(\theta)|^2 d\theta \\
&= \frac{1}{|\Theta|} \int_{\Theta} \left| \sum_{\substack{\omega \in T, \\ \#\omega > \ell}} \Phi_{\omega}(\theta) d_{\omega} \right|^2 d\theta \\
&= \sum_{\substack{(\omega, \omega') \in T^2, \\ h(\omega) = h(\omega'), \\ \#\omega > \ell, \#\omega' > \ell}} 2^{-|h(\omega)|} d_{\omega} d_{\omega'}.
\end{aligned} 
\quad (14) $$

The function $h$ transforms a vector $\omega$ into a vector $h(\omega)$ by assigning the same value, 0, to 0, $0_I$, and $0_Z$, while leaving the rest unchanged. The condition $h(\omega) = h(\omega')$ ensures that only few terms remain in the error sum, and results from the averaging over the rotation angles. The magnitude $|h(\omega)|$ counts the number of non-zero elements in $\omega$, i.e. the number of times processes $D_{\pm 1}$ have been encountered during the backpropagation on path $\omega$. The error depends on all pairs of paths $(\omega, \omega')$ that satisfy certain constraints. For each of these pairs, it is possible to factor out the noise introduced by the processes from $d_\omega$. By defining $Q(\omega) = \prod_{i=1}^m Q(\omega_i)$ as the noise accumulated through the splitting, where $Q(0_Z) = 1-\gamma$, $Q(0_I) = \gamma$, $Q(+1) = Q(-1) = \sqrt{1-\gamma}$, and $Q(0) = 1$, we can rewrite $d_\omega$ as $d_\omega = Q(\omega)d_\omega^0$. By definition of $d_\omega^0$, it is clear that for all $\omega$, $|d_\omega^0| \le 1$. This allows us to derive a simple bound on the error:

$$ \Delta^2(\tilde{f}, f) \le \sum_{\substack{(\omega, \omega') \in T^2, \\ h(\omega) = h(\omega'), \\ \#\omega > \ell, \#\omega' > \ell}} 2^{-|h(\omega)|} Q(\omega)Q(\omega') \quad (15) $$

The error can then be decomposed in two sums, one where $\omega = \omega'$, and one where $\omega \ne \omega'$ but $h(\omega) = h(\omega')$. We show that in both cases, splits occuring through the $D_{\pm 1}$ processes dampen the approximation error by $(1-\gamma)$ and splits occuring through the processes $D_{0I}$ and $D_{0Z}$ leave the approximation error unchanged. To see this,

<figure><img src="image_4.png" /><figcaption>FIG. 3: Branch of a tree resulting from the backpropagation of a Pauli observable, which splits through the process $D_{\pm 1}$.</figcaption></figure>

take a branch of the tree which has had $j$ splits (denote the branch length at this point by $i$), and go through the tree until one more split occurs on this branch. We show that the two resulting branches' contribution to the error term is less than the contribution of the original branch (by computing $Q(\omega)Q(\omega')$ and $2^{-|h(\omega)|}$) in both the cases introduced before, illustrated in Figure 3 and 4. When a pair of paths $(\omega, \omega')$ has a nonzero contribution to the error term, i.e. when $h(\omega) = h(\omega')$, we say that these paths "interact".

* $\omega_{1,...,i} = \omega'_{1,...,i}$: If the split occurs through processes $D_{\pm 1}$, then the two resulting branches can only interact with themselves in order to satisfy $h(\omega) = h(\omega')$, such that the damping introduced is $2 \times (1-\gamma)2^{-1} = 1-\gamma$. If the split comes from processes $D_{0I}$ and $D_{0Z}$, then there will be four terms in the sum as all pairs of paths will satisfy $h(\omega) = h(\omega')$. The approximation error remains unchanged since $\gamma^2 + 2\gamma(1-\gamma) + (1-\gamma)^2 = 1$.

* $\omega_{1,...,i} \ne \omega'_{1,...,i}$: The cross interactions between different branches are more subtle, as a different split (or no split) may occur on one of the branches. First, note that if one of the branches does not split through process $D_0$, then to satisfy $h(\omega) = h(\omega')$, the only valid split for the other branch is through processes $D_{0Z}$ and $D_{0I}$. In this case, the two new branches will interact with the single branch and the error is unchanged since $1 \times \gamma + 1 \times (1-\gamma) = 1$. If both branches split, then to satisfy $h(\omega) = h(\omega')$, they either both split through processes $D_{\pm 1}$ or both split through processes $D_{0Z}$ and $D_{0I}$. The first case leads to two interactions ($-1$ with $-1$ and $+1$ with $+1$), introducing the same damping as before $1-\gamma$. Similarly, the second case leaves the error unchanged since $\gamma^2 + 2\gamma(1-\gamma) + (1-\gamma)^2 = 1$.

Everytime a split occurs through processes $D_{\pm 1}$, its total contribution is damped by a factor $1-\gamma$. Similarly, when a split occurs through processes $D_{0Z}$ or $D_{0I}$, the approximation error is unchanged. Counting for each discarded branch the number of times that processes $D_{\pm 1}$ have been encountered allows us to retrieve Equation 12. The convergence of the Monte-Carlo approach, Equation
<figure><img src="image_5.png" /><figcaption>FIG. 4: Two different branches of a tree resulting from the backpropagation of a Pauli observable that satisfy $h(\omega_1,...,i) = h(\omega'_1,...,i)$. To satisfy $h(\omega) = h(\omega')$ on the entire branches, only few pairs of processes are possible.</figcaption></figure>

13, can be derived similarly. The error made is decomposed in an error term arising from the Monte-Carlo approximation which scales in $1/\sqrt{K}$, our estimator being unbiased, and a truncation error (see Appendix C). The truncation error for each of the sampled trees is easily bounded by $(1 - \gamma)^{(\ell+1)/2}$ using the same analysis and the fact that all splits in those trees occur through processes $D_{\pm 1}$.

IV. CONCLUSION

In this work, we have shown that Pauli backpropagation methods allow us to efficiently compute the expectation value of Pauli observables under non-unital noise.

This simulation method creates a surrogate of the energy landscape by only backpropagating the observable once in the Heisenberg picture. Expectation values for given angles can then be quickly obtained by just computing the trigonometric monomials. This simulation method is parallelizable and memory efficient, and does not require any assumptions on the locality of the observable, on the geometry of the underlying quantum hardware or on the depth of the quantum circuit at hand. This algorithm is particularly adapted to PQC s, as the structure strongly resemble the one commonly used for such circuits. Thus, our work establishes that Pauli backpropagation is an effective and efficient method to simulate noisy quantum devices under realistic hardware noise assumptions.

Furthermore, even though results like the quantum refrigerator [18] show that the hardware noise can have significant impact on its computational power, our result clearly indicate that this difference only arises for highly structured circuits, like the ones used in fault-tolerant constructions. Finally, it would be interesting to show similar results to simulate noisy quantum simulators.

Acknowledgements: The authors thank Manuel S. Rudolph, Zoë Holmes and Antonio Anna Mele for useful feedback and discussion. This work benefited from a government grant managed by the Agence Nationale de la Recherche under the Plan France 2030 with the reference ANR-22-PETQ-0007. AA acknowledges support from the Sandoz Family Foundation-Monique de Meuron program for Academic Promotion.

[1] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Nature **618**, 500 (2023).

[2] H. Yu, Y. Zhao, and T.-C. Wei, Phys. Rev. Res. **5**, 013183 (2023).

[3] A. Carrera Vazquez, C. Tornow, D. Ristè, S. Woerner, M. Takita, and D. J. Egger, Nature **636**, 75-79 (2024).

[4] J. R. Glick, T. P. Gujarati, A. D. Córcoles, Y. Kim, A. Kandala, J. M. Gambetta, and K. Temme, Nature Physics **20**, 479-483 (2024).

[5] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics **3**, 625-644 (2021).

[6] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C.

Kwek, and A. Aspuru-Guzik, Reviews of Modern Physics **94**, 10.1103/revmodphys.94.015004 (2022).

[7] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014), arXiv:1411.4028 [quant-ph].

[8] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O'Brien, Nature Communications **5**, 10.1038/ncomms5213 (2014).

[9] M. Cerezo, K. Sharma, A. Arrasmith, and P. J. Coles, npj Quantum Information **8**, 10.1038/s41534-022-00611-6 (2022).

[10] J. Bowles, S. Ahmed, and M. Schuld, Better than classical? the subtle art of benchmarking quantum machine learning models (2024), arXiv:2403.07059 [quant-ph].

[11] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Nature Communications **12**, 10.1038/s41467-021-27045-6 (2021).
[12] D. Stilck França and R. García-Patrón, Nature Physics 17, 1221–1227 (2021).

[13] R. Takagi, S. Endo, S. Minagawa, and M. Gu, npj Quantum Information 8, 10.1038/s41534-022-00618-z (2022).

[14] T. Schuster, C. Yin, X. Gao, and N. Y. Yao, A polynomial-time classical algorithm for noisy quantum circuits (2024), arXiv:2407.12768 [quant-ph].

[15] A. A. Mele, A. Angrisani, S. Ghosh, S. Khatri, J. Eisert, D. S. França, and Y. Quek, Noise-induced shallow circuits and absence of barren plateaus (2024), arXiv:2403.13927 [quant-ph].

[16] L. Chirolli and G. Burkard, Advances in Physics 57, 225–285 (2008).

[17] R. Blume-Kohout, J. K. Gamble, E. Nielsen, K. Rudinger, J. Mizrahi, K. Fortier, and P. Maunz, Nature Communications 8, 10.1038/ncomms14485 (2017).

[18] M. Ben-Or, D. Gottesman, and A. Hassidim, Quantum refrigerator (2013), arXiv:1301.1995 [quant-ph].

[19] E. Fontana, M. S. Rudolph, R. Duncan, I. Rungger, and C. Cirstoiu, Classical simulations of noisy variational quantum circuits (2023), arXiv:2306.05400 [quant-ph].

[20] M. S. Rudolph, E. Fontana, Z. Holmes, and L. Cincio, Classical surrogate simulation of quantum systems with lowesa (2023), arXiv:2308.09109 [quant-ph].

[21] T. Begušić, J. Gray, and G. K.-L. Chan, Science Advances 10, eadk4321 (2024), https://www.science.org/doi/pdf/10.1126/sciadv.adk4321.

[22] A. Angrisani, A. Schmidhuber, M. S. Rudolph, M. Cerezo, Z. Holmes, and H.-Y. Huang, Classically estimating observables of noiseless quantum circuits (2024), arXiv:2409.01706 [quant-ph].

[23] Y. Shao, F. Wei, S. Cheng, and Z. Liu, Phys. Rev. Lett. 133, 120603 (2024).

[24] S. Lerch, R. Puig, M. S. Rudolph, A. Angrisani, T. Jones, M. Cerezo, S. Thanasilp, and Z. Holmes, Efficient

quantum-enhanced classical simulation for patches of quantum landscapes (2024), arXiv:2411.19896 [quant-ph].

[25] P. Bermejo, P. Braccia, M. S. Rudolph, Z. Holmes, L. Cincio, and M. Cerezo, Quantum convolutional neural networks are (effectively) classically simulable (2024), arXiv:2408.12739 [quant-ph].

[26] P. Rall, D. Liang, J. Cook, and W. Kretschmer, Phys. Rev. A 99, 062337 (2019).

[27] D. Aharonov, X. Gao, Z. Landau, Y. Liu, and U. Vazirani, in *Proceedings of the 55th Annual ACM Symposium on Theory of Computing*, STOC 2023 (Association for Computing Machinery, New York, NY, USA, 2023) p. 945–957.

[28] G. González-García, J. I. Cirac, and R. Trivedi, Pauli path simulations of noisy quantum circuits beyond average case (2024), arXiv:2407.16068 [quant-ph].

[29] A. Tanggara, M. Gu, and K. Bharti, arXiv preprint arXiv:2405.00789 (2024).

[30] A. Angrisani, A. A. Mele, M. S. Rudolph, M. Cerezo, and Z. Holmes, Simulating quantum circuits with arbitrary local noise using pauli propagation (2025).

[31] S. Bravyi and D. Gosset, Physical Review Letters 116, 10.1103/physrevlett.116.250501 (2016).

[32] S. Bravyi and A. Kitaev, Phys. Rev. A 71, 022316 (2005).

[33] J. M. Chow, J. M. Gambetta, A. D. Córcoles, S. T. Merkel, J. A. Smolin, C. Rigetti, S. Poletto, G. A. Keefe, M. B. Rothwell, J. R. Rozen, M. B. Ketchen, and M. Stefken, Phys. Rev. Lett. 109, 060501 (2012).

[34] S. Aaronson and D. Gottesman, Physical Review A 70, 10.1103/physreva.70.052328 (2004).

[35] M. B. Ruskai, S. Szarek, and E. Werner, An analysis of completely-positive trace-preserving maps on 2x2 matrices (2001), arXiv:quant-ph/0101003 [quant-ph].
Appendix A: Average Pauli backpropagation algorithm

In this section we present the pseudocode for our deterministic Pauli backpropagation algorithm in Alg. 1 and empirically evaluate how fast it converges.

<caption>**Algorithm 1: [LOWESA-AD] Simulation of Pauli string observables under non-unital noise**</caption>

**Input:** Quantum circuit made of $m$ alternating layers of Clifford unitaries and $R_z$ rotations, affected by non-unital noise, and defined by process modes $\{D_\omega\}$; measurement Pauli operator $P$; cut-off parameter $\ell$.
**Output:** $\tilde{f}(\theta)$, an approximation of $f(\theta)$ and $r \le \min_{\omega/\#\omega \ge \ell} |h(\omega)|$.
**Procedure LOWESA-AD:** Given a quantum circuit with rotation parameters $\theta$ and Pauli observable $P$, outputs an approximation of the expectation value of $P$, by first creating a surrogate of the energy landscape using Low-weight coefficients and then adjusting the angles.
1: $\tilde{f}(\theta) \leftarrow 0$
2: Run Low-weight coefficients to calculate $d_\omega = \langle\langle 0 | D_\omega^T | P \rangle\rangle$ for all $\omega$ such that $\#\omega \le \ell$.
3: **for all** $\omega$ such that $\#\omega \le \ell$ with non-zero $d_\omega$ **do**
4:     $\tilde{f}(\theta) \leftarrow \tilde{f}(\theta) + d_\omega \Phi_\omega(\theta)$
5: **end for**
6: Return $\tilde{f}(\theta)$ and $r = \min_\omega r_\omega \le \min_{\omega/\#\omega \ge \ell} |h(\omega)|$
**end Procedure**

**Procedure Low-weight coefficients:** Given a quantum circuit and Pauli observable $P$, outputs the backpropagated tree coefficients $d_\omega$ for all $\omega$ such that $\#\omega \le \ell$ in time $O(n^2m2^\ell)$. The tree is constructed iteratively as follows: at each layer of the quantum circuit, the tree $T_i$ is extended at each branch $\omega$ based on the Clifford layer $C_{m-i+1}$ and the Pauli operator on the rotation qubit $q_{m-i+1}$. Note that the initial tree $T_0$ is made of a single node corresponding to the Pauli observable $P$. The branches of the tree that have split more than $\ell$ times are discarded.
1: **for** $i=1$ **to** $m$ **do**
2:     **for** $\omega$ **in** $T_{i-1}$ **do**
3:         $P_\omega \leftarrow C_{m-i+1}^\dagger P_\omega C_{m-i+1}$, up to a phase $\phi_{\omega_i} \in \{\pm 1\}$ which is stored.
4:         **if** $P_\omega^{q_{m-i+1}} = I$
5:             Add branch $\omega$ in $T_i$ with $\omega_i = 0$
6:             $d_{\omega_i}, P_\omega \leftarrow D_{\omega_i}^\dagger (P_\omega)$
7:             **else if** $P_\omega^{q_{m-i+1}} = Z$ and $\#\omega < \ell$ **then**
8:                 Split into two branches $\omega$ and $\omega'$ in $T_i$ with $\omega_i = 0_I$ and $\omega'_i = 0_Z$
9:                 $d_{\omega_i}, P_\omega \leftarrow D_{\omega_i}^\dagger (P_\omega)$
10:              $d_{\omega'_i}, P_{\omega'} \leftarrow D_{\omega'_i}^\dagger (P_\omega)$
11:             $\#\omega, \#\omega' \leftarrow +1$
12:         **else if** $P_\omega^{q_{m-i+1}} \in \{X, Y\}$ and $\#\omega < \ell$ **then**
13:             Split into two branches $\omega$ and $\omega'$ in $T_i$ with $\omega_i = +1$ and $\omega'_i = -1$
14:             $d_{\omega_i}, P_\omega \leftarrow D_{\omega_i}^\dagger (P_\omega)$
15:             $d_{\omega'_i}, P_{\omega'} \leftarrow D_{\omega'_i}^\dagger (P_\omega)$
16:             $\#\omega, \#\omega' \leftarrow +1$
17:         **else if** $P_\omega^{q_{m-i+1}} \in \{X, Y, Z\}$ and $\#\omega = \ell$ **then**
18:             Discard path $\omega$
19:             $r_\omega = |h(\omega)|$
20:         **end if**
21:     **end for**
22: **end for**
23: **for** $\omega$ **in** $T_m$ **do**
24:     $P_\omega \leftarrow C_0^\dagger P_\omega C_0$ and phase $\phi_{\omega_0} \in \{\pm 1\}$
25: **end for**
26: Return $d_\omega = (\prod_{i=1}^m d_{\omega_i} \phi_{\omega_i}) \phi_{\omega_0} (\langle\langle 0 | P_\omega \rangle\rangle)$ and $r_\omega$ for all $\omega$ in $T_m$
**end Procedure**

As described in Theorem 1, it is sufficient to have access to the quantity $\min_{\omega/\#\omega \ge \ell} |h(\omega)|$ to get a certificate on the $L^2$-error made by our approximation. This quantity, which Algorithm 1 computed a lower bound to, is expected to scale linearly with $\ell$ in most cases. We present numerical evidence of such scaling by implementing our Algorithm on common PQC instances, namely QAOA for random 3-regular graphs. For each of these graph instances and associated QAOA quantum circuit, we backpropagate the observable $Z_i Z_j$, keeping track in the resulting tree of the quantity
$r \le \min_{\omega/\#\omega \ge \ell} |h(\omega)|$ for different values of $\ell$. The results are plotted in Figure 5.

<figure><img src="image_6.png" /><figcaption>FIG. 5: Lower bound $r$ on the expected number of splits through processes $\mathbf{D}_{\pm 1}$ for a given cut-off parameter $\ell$. The quantum circuits considered are QAOA circuits on 3-regular graphs, with Pauli observables $Z_i Z_j$.</figcaption></figure>

## Appendix B: Proof of convergence

Recall that the noisy circuits considered are made of alternating layers of Clifford unitary and single qubit $R_z$ rotations, such that,

$$ U_{\theta} = \left( \bigcirc_{i=1}^{m} C_i \circ \mathcal{R}_{z}^{(q_i)}(\theta_i) \circ \mathcal{N}_{AD} \right) \circ C_0 \qquad (B1) $$

The noisy rotations can be decomposed in the PTM using the structure of the noise,

$$ \mathbf{R}_z \cdot \mathbf{N}_{AD} = \mathbf{D}_0 + \mathbf{D}_{0z} + \mathbf{D}_{0i} + \mathbf{D}_1 \cos \theta + \mathbf{D}_{-1} \sin \theta \qquad (B2) $$

where we have defined the quantum processes $\mathbf{D}_0 = |I\rangle\langle I|$, $\mathbf{D}_{0z} = (1-\gamma)|Z\rangle\langle Z|$, $\mathbf{D}_{0i} = \gamma|I\rangle\langle Z|$, $\mathbf{D}_1 = \sqrt{1-\gamma}(|X\rangle\langle X| + |Y\rangle\langle Y|)$ and $\mathbf{D}_{-1} = \sqrt{1-\gamma}(|Y\rangle\langle X| - |X\rangle\langle Y|)$, with $|I\rangle = [1 \ 0 \ 0 \ 0]^T$, $|X\rangle = [0 \ 1 \ 0 \ 0]^T$... Our quantum circuit can then be further decomposed by summing all possible processes and keeping track of the rotation angles with the trigonometric function $\Phi_{\omega}(\theta) = \prod_{i=1}^{m} \phi_{\omega_i}(\theta_i)$, with $\phi_0(\theta_0) = \phi_{0z}(\theta_0) = \phi_{0i}(\theta_0) = 1$, $\phi_1(\theta_1) = \cos \theta_1$, and $\phi_{-1}(\theta_1) = \sin \theta_1$.

$$ U_{\theta} = \sum_{\omega \in \{0, 0_z, 0_i, 1, -1\}^m} \Phi_{\omega}(\theta) D_{\omega} \qquad (B3) $$

Each process $D_{\omega}$ is defined in the Schrödinger picture by $D_{\omega} = (\bigcirc_{i=1} C_i \circ D_{\omega_i}) \circ C_0$. When backpropagating a Pauli string $P$, the rotation acts only on a single qubit $q_i$. Depending on the Pauli operator at qubit $q_i$, at most two processes will be valid: $\mathbf{D}_0$ for $I$, $\mathbf{D}_{0z}$ and $\mathbf{D}_{0i}$ for $Z$, and $\mathbf{D}_{\pm 1}$ for $X$ and $Y$. This suggests the structure of a rooted tree $T = (V, E)$. The set of nodes $V$ represents the backpropagated Pauli strings at each layer, and the edges represent the valid processes $\mathbf{D}_{\omega_i}$ for the given Paulis. Our expectation value can then be rewritten as:

$$ f(\theta) = \mathrm{tr}(U_{\theta}[|0\rangle\langle 0|]P) = \sum_{\omega \in T} d_{\omega} \Phi_{\omega}(\theta) \qquad (B4) $$

where the Fourier coefficients for each branch of the tree are defined by $d_{\omega} = \mathrm{tr}(D_{\omega}[|0\rangle\langle 0|]P) = \langle\langle P|\mathbf{D}_{\omega}|0\rangle\rangle$. By backpropagating the Pauli string in this manner, we retrieve its expectation value, though at the cost of exponential
computation time. The simulation requires updating the Clifford unitary for all paths of the tree at each layer. As there are m layers and at most 2^m branches, the runtime of this exact simulation is at most O(n^2m2^m). Furthermore, as the Pauli observable is backpropagated through the quantum circuit, noise accumulates on each of the branches. Denote Q(ω) the noise damping accumulated throughout the path ω, with Q(ω) = Π_i^m Q(ω_i) where Q(+1) = Q(-1) = √1 - γ, Q(0) = 1, Q(0_I) = γ and Q(0_Z) = 1 - γ. For some d_ω^0 ∈ {-1, 0, 1}, we can rewrite d_ω = Q(ω)d_ω^0, such that,

$$f(\theta) = \mathrm{tr}(\mathcal{U}_{\theta}[0]\langle 0|P) = \sum_{\omega \in T} Q(\omega) d_{\omega}^{0} \Phi_{\omega}(\theta) \qquad (B5)$$

Each time a branch in the tree splits, the resulting Paulis are damped by the noise. Once sufficient damping has occurred, the contribution of this branch to the expectation value can be neglected. We show that keeping only paths in the tree that have split fewer than a constant number of times ℓ yields a good approximation. More precisely, let #ω denote the number of times a given path ω has split. We approximate the function f by the quantity:

$$\tilde{f}(\theta) = \sum_{\substack{\omega \in T \\ \#\omega \le \ell}} d_{\omega} \Phi_{\omega}(\theta) = \sum_{\substack{\omega \in T \\ \#\omega \le \ell}} Q(\omega) d_{\omega}^{0} \Phi_{\omega}(\theta) \qquad (B6)$$

The metric used to judge the quality of our approximation error will be the L²-norm between f and $\tilde{f}$ over the range of angles $\Theta = [0, 2\pi]^m$.

$$\begin{aligned} \Delta^2(\tilde{f}, f) &= \frac{1}{|\Theta|} \int_{\Theta} |\tilde{f}(\theta) - f(\theta)|^2 d\theta \\ &= \frac{1}{|\Theta|} \int_{\Theta} \left| \sum_{\substack{\omega \in T \\ \#\omega > \ell}} \Phi_{\omega}(\theta) d_{\omega} \right|^2 d\theta \end{aligned} \qquad (B7)$$

This error metric is convenient for several reasons. First, note that there is no hope of having a polynomial-time simulation algorithm under non-unital noise for all circuits, as celebrated results [18] have shown that it is possible to compute under such noise for an exponential amount of time and perform error correction. The L²-norm error gives us convergence on average over the rotation angles, which does not contradict previous results since, for a given set of angles, the simulation might not converge. Furthermore, this simulation method creates a surrogate of the parametrized expectation value, which can then be used by plugging in different angles. The L²-norm error therefore tells us how good our approximation is over the entire “landscape”. Finally, when computing this metric, most of the paths will be orthogonal to each other, resulting in only a few terms remaining in the error which can then be analytically bounded.

**Lemma 2.** Denote h a function that transforms a vector ω to a vector h(ω) where 0, 0_Z and 0_I are assigned the same value 0, the rest of the components being unchanged. Denote #ω the number of times the path ω has split and |h(ω)| the number of non-zero elements of h(ω). The error made by the algorithm can be expressed as

$$\Delta^2(\tilde{f}, f) = \sum_{\substack{(\omega, \omega') \in T^2 \\ h(\omega) = h(\omega') \\ \#\omega > \ell, \#\omega' > \ell}} Q(\omega)Q(\omega')2^{-|h(\omega)|}d_{\omega}^{0}d_{\omega'}^{0} \qquad (B8)$$

*Proof.* Using known properties on the integral of product of trigonometric functions, one can compute the useful quantity:

$$\frac{1}{|\Theta|} \int_{\Theta} \Phi_{\omega}(\theta) \Phi_{\omega'}(\theta) d\theta = \frac{1}{|\Theta|} \int_{\Theta} \prod_{i=1}^{m} \phi_{\omega_i}(\theta_i) \prod_{i=1}^{m} \phi_{\omega'_i}(\theta'_i) d\theta = 2^{-|h(\omega)|} \delta_{h(\omega)=h(\omega')} \qquad (B9)$$

The approximation made by our algorithm can therefore be rewritten as:

$$\Delta^2(\tilde{f}, f) = \frac{1}{|\Theta|} \int_{\Theta} \left| \sum_{\substack{\omega \in T \\ \#\omega > \ell}} \Phi_{\omega}(\theta) d_{\omega} \right|^2 d\theta = \sum_{\substack{(\omega, \omega') \in T^2 \\ h(\omega) = h(\omega') \\ \#\omega > \ell, \#\omega' > \ell}} 2^{-|h(\omega)|} d_{\omega} d_{\omega'} \qquad (B10)$$

Finally, factoring out the noise Q(ω) and Q(ω') from the paths ω and ω' gives us the aforementioned equality. □

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6: (a) Tree corresponding to the backpropagation of a single observable through the circuit. (b) Example of a pair of paths $(\omega, \omega')$ such that $v(\omega, \omega') = 2$.</figcaption></figure>

By definition of $d_{\omega}^{0}$, it is clear that for all $\omega$, $|d_{\omega}^{0}| \le 1$. This allows us to derive a simple bound on the error made:

$$ \Delta^2(\tilde{f}, f) \le \sum_{\substack{(\omega, \omega') \in T^2 \\ h(\omega) = h(\omega') \\ \#\omega > \ell, \#\omega' > \ell}} Q(\omega)Q(\omega')2^{-|h(\omega)|} \quad (B11) $$

Our approximation error depends on the sum of pairs of paths $(\omega, \omega')$ of weights greater than $\ell$ satisfying some constraints. To quantify the contribution of each pair of paths, we introduce the quantity $v(\omega, \omega')$, that counts the number of positions where either $\omega$ or $\omega'$ or both have split, see Figure 6. The error can clearly be bounded by:

$$ \Delta^2(\tilde{f}, f) \le \sum_{\substack{(\omega, \omega') \in T^2 \\ h(\omega) = h(\omega') \\ v(\omega, \omega') > \ell}} Q(\omega)Q(\omega')2^{-|h(\omega)|} \quad (B12) $$

We show that splits that occur through processes $D_{\pm 1}$ dampen the approximation error by $1 - \gamma$ and splits that occur through the pair of processes $D_{0z}/D_{0i}$ leave the approximation error unchanged. Since random Clifford maps non-identity Paulis at random, adding a single-qubit random Clifford gate before the rotation allows us to derive the following bound:

**Proposition 3.** Consider quantum circuits as defined previously, i.e. made of alternating layers of Clifford unitaries and 1-qubit rotations followed by a single random Clifford gate. The approximation error made by Algorithm 1 with cut-off $\ell$ is bounded by:

$$ \mathbb{E}[\Delta^2(\tilde{f}, f)] \le [1 - 2\gamma/3]^{\ell+1}, \quad (B13) $$

where the expectation is over the choice of the single-qubit random Clifford gates.

*Proof.* For $i \in \{0, \dots, m\}$ and $\omega_{1, \dots, i}$ and $\omega'_{1, \dots, i}$ such that $h(\omega_{1, \dots, i}) = h(\omega'_{1, \dots, i})$, we define

$$ \alpha_j(\omega_{1, \dots, i}, \omega'_{1, \dots, i}) = \sum_{\substack{\omega_{i+1, \dots, m}, \omega'_{i+1, \dots, m} \\ h(\omega) = h(\omega') \\ v(\omega_{i+1, \dots, m}, \omega'_{i+1, \dots, m}) \ge j}} Q(\omega_{i+1, \dots, m})Q(\omega'_{i+1, \dots, m})2^{-|h(\omega_{i+1, \dots, m})|} $$

Note that the quantity we wish to bound is $\alpha_{\ell+1}$ with the empty word ($i=0$). We show by induction on $j$ that for any $i$ and $\omega_{1, \dots, i}, \omega'_{1, \dots, i}$, we have $\alpha_j(\omega_{1, \dots, i}, \omega'_{1, \dots, i}) \le [1 - 2\gamma/3]^j$.

The initialisation $j=0$ is straightforward as the expectation value of the Pauli string is bounded by 1. Assume now that the property is satisfied for $j$ and notice that both $Q$ and $2^{-|h|}$ can be expressed as products. We introduce
$i'$ which corresponds to the index of the next split after $i$. Note that $\omega_{1,...,i'-1}$ is the (unique) path extending $\omega_{1,...,i}$ up to $i' - 1$ similarly for $\omega'$. In addition, we have $Q(\omega_{i+1,...,i'-1}) = Q(\omega_{i+1,...,i'-1}) = 1$ and $2^{-|h(\omega_{i+1,...,i'-1})|} = 1$ which allows us to write,

$$
\begin{align*}
\alpha_{j+1}(\omega_{1,...,i}, \omega'_{1,...,i}) &= \sum_{\substack{\omega_{i+1,...,m}, \omega'_{i+1,...,m} \\ h(\omega)=h(\omega') \\ v(\omega_{i+1,...,m}, \omega'_{i+1,...,m}) \ge j+1}} Q(\omega_{i+1,...,m}) Q(\omega'_{i+1,...,m}) 2^{-|h(\omega_{i+1,...,m})|} \\
&= \sum_{\substack{\omega_{i+1,...,i'}, \omega'_{i+1,...,i'} \\ h(\omega)=h(\omega')}} Q(\omega_{i+1,...,i'}) Q(\omega'_{i+1,...,i'}) 2^{-|h(\omega_{i+1,...,i'})|} \sum_{\substack{\omega'_{i'+1,...,m}, \omega'_{i'+1,...,m} \\ h(\omega)=h(\omega') \\ v(\omega'_{i'+1,...,m}, \omega'_{i'+1,...,m}) \ge j}} Q(\omega_{i'+1,...,m}) Q(\omega'_{i'+1,...,m}) 2^{-|h(\omega'_{i'+1,...,m})|} \\
&= \sum_{\substack{\omega_{i+1,...,i'}, \omega'_{i+1,...,i'} \\ h(\omega)=h(\omega')}} Q(\omega_{i+1,...,i'}) Q(\omega'_{i+1,...,i'}) 2^{-|h(\omega_{i+1,...,i'})|} \alpha_j(\omega_{1,...,i'}, \omega'_{1,...,i'})
\end{align*}
$$

Using the recursive hypothesis on $\alpha_j(\omega_{1,...,i'}, \omega'_{1,...,i'})$,

$$
\alpha_{j+1}(\omega_{1,...,i}, \omega'_{1,...,i}) \le [1 - 2\gamma/3]^j \sum_{\substack{\omega_{i+1,...,i'}, \omega'_{i+1,...,i'} \\ h(\omega)=h(\omega')}} Q(\omega_{i+1,...,i'}) Q(\omega'_{i+1,...,i'}) 2^{-|h(\omega_{i+1,...,i'})|} \quad (B14)
$$

It remains to bound the sum for pairs of path $(\omega_{i+1,...,i'}, \omega'_{i+1,...,i'})$. Note that by construction, there is at most one split occuring on either or both of these paths. We consider two distinct cases, which we will show lead to the same bound. The first one is the case where $\omega_{i+1,...,i'-1} = \omega'_{i+1,...,i'-1}$ and the second the one where $\omega_{i+1,...,i'-1} \ne \omega'_{i+1,...,i'-1}$ but $h(\omega_{i+1,...,i'-1}) = h(\omega'_{i+1,...,i'-1})$.

1. $\omega_{i+1,...,i'-1} = \omega'_{i+1,...,i'-1}$: In this case, illustrated in Figure 3, it is clear that the same split occurs on both paths. If the split occurs through processes $D_{\pm 1}$ at $i'$, then the two resulting branches can only interact with themselves in order to satisfy $h(\omega) = h(\omega')$, and have contribution $Q(+1)Q(+1)2^{-|h(+1)|} + Q(-1)Q(-1)2^{-|h(-1)|} = (1-\gamma)$ in the sum. If the split occurs through $D_{0Z}/D_{0i}$, then there will be four terms in the sum as all pairs of paths will satisfy $h(\omega) = h(\omega')$. The result of the sum will be $Q(0Z)Q(0Z)2^{-|h(0Z)|} + Q(0I)Q(0I)2^{-|h(0I)|} + Q(0Z)Q(0I)2^{-|h(0Z)|} + Q(0I)Q(0Z)2^{-|h(0I)|} = \gamma^2 + 2\gamma(1-\gamma) + (1-\gamma)^2 = 1$. Since random Clifford gates map the non-identity Paulis uniformly at random between $\{X, Y, Z\}$, the probability of a split occuring through $D_{\pm 1}$ is $2/3$ and $1/3$ through $D_{0Z}/D_{0i}$. This allows us to compute the expected contribution of all the paths after the split:

$$
\mathbb{E} \left[ \sum_{\substack{\omega_{i+1}, \dots, i', \omega'_{i+1}, \dots, i' \\ h(\omega) = h(\omega')}} Q(\omega_{i+1}, \dots, i') Q(\omega'_{i+1}, \dots, i') 2^{-|h(\omega_{i+1}, \dots, i')|} \right] \le [1 - 2\gamma/3] \quad (B15)
$$

2. $\omega_{i+1,...,i'-1} \neq \omega'_{i+1,...,i'-1}$: In this case, illustrated in Figure 4, we need to be more careful as the split occuring on each branch might be different, and two cases are possible. Either both the branches split (the splits might not be the same here), or only one splits, the other encountering process $D_0$. Let us consider the latter case first. If one of the branches do not split, then maintaining $h(\omega) = h(\omega')$ requires the other branch splitting through $D_{0Z}/D_{0i}$. Once again, because of the random Clifford, this occurs with probability $1/3$, and leads to a sum of $Q(0)Q(0Z)2^{-|h(0)|} + Q(0)Q(0I)2^{-|h(0)|} = 1 \times \gamma + 1 \times (1-\gamma) = 1$.

If both of the branches split, we need to differentiate based on the splits occuring. If one of the branches split through $D_{\pm 1}$ and the other through $D_{0Z}/D_{0i}$, then $h(\omega) \neq h(\omega')$. If both split through $D_{\pm 1}$, then the 4 new branches will lead to two terms satisfying $h(\omega) = h(\omega')$, carrying a damping of $Q(+1)Q(+1)2^{-|h(+1)|} + Q(-1)Q(-1)2^{-|h(-1)|} = (1-\gamma)$. Finally, if both split through processes $D_{0Z}/D_{0i}$, then the four terms arising are once again $Q(0Z)Q(0Z)2^{-|h(0Z)|} + Q(0I)Q(0I)2^{-|h(0I)|} + Q(0Z)Q(0I)2^{-|h(0Z)|} + Q(0I)Q(0Z)2^{-|h(0I)|} = \gamma^2 + 2\gamma(1-\gamma) + (1-\gamma)^2 = 1$. Because of the random Clifford, the probability of both split occuring through $D_{0Z}/D_{0i}$ is less than the probability of one of the paths splitting through $D_{\pm 1}$, i.e. less than $1/3$. Once again,
the expected damping is upper bounded by $[1 - 2\gamma/3]$. This bound also holds for the case where only one of the branches split such that,

$$ \mathbb{E}[\sum_{\substack{\omega_{i+1}, \dots, i', \omega'_{i+1}, \dots, i' \\ h(\omega) = h(\omega')}} Q(\omega_{i+1}, \dots, i') Q(\omega'_{i+1}, \dots, i') 2^{-|h(\omega_{i+1}, \dots, i')|}] \le [1 - 2\gamma/3] \quad (B16) $$

Everytime a split occurs on a branch, its expected total contribution is damped by a factor $[1 - 2\gamma/3]$ which concludes the recursive argument and the proof.

$$ \mathbb{E}[\Delta^2(f, \tilde{f})] \le [1 - 2\gamma/3]^{\ell+1} \quad (B17) $$

Obtaining Equation 12 of Theorem 1 can be done using the same proof technique. We've seen that splits occurring through processes $D_{\pm 1}$ dampen the approximation error by $1 - \gamma$ and that splits occurring through processes $D_{0z}/D_{0i}$ leave it unchanged. By counting the minimum number of times all discarded paths have split through processes $D_{\pm 1}$, which we denote $r$, we get that,

$$ \Delta^2(f, \tilde{f}) \le (1 - \gamma)^r \quad (B18) $$

Note that we expect $r$ to scale linearly in $\ell$ for most circuits, as shown in Figure 5 for QAOA [7] instances on 3-regular graphs. $\square$

## Appendix C: Monte-Carlo Pauli backpropagation algorithm

For most quantum circuits, we have established that an extension of the standard LOWESA [19] algorithm allows to approximate any observable in polynomial time in the $L^2$-norm. Provided one is willing to pay an additional sampling toll, we show through Algorithm 2 that a Monte-Carlo approach allows us to obtain classical simulability for any circuit. Using the tools previously introduced, this approach is straightforward to explain.

Consider the complete tree $T$ created by the backpropagation of the Pauli string $P$ through the quantum circuit. As seen previously, splits that occur in this tree due to processes $D_{\pm 1}$ dampen the expectation value of $P$, whereas splits occurring through processes $D_{0z}$ and $D_{0i}$ might leave the norm of the expectation value unchanged, in the case of amplitude damping noise alone. However, these splits still create additional branches, increasing the computational cost of the simulation.

To alleviate this problem, we randomly sample only one of these processes at each split. The process $D_{0z}$ is kept with probability $1 - \gamma$ and $D_{0i}$ with probability $\gamma$. To maintain an unbiased estimator, we make sure to remove the associated noise factors $\gamma$ or $1 - \gamma$ from $d_\omega$. As a result, the nonphysical trees generated by this sampling approach only contain splits that effectively dampen the expectation value of $P$ by $\sqrt{1 - \gamma}$. We can then approximate the trees by truncating paths that have split excessively.

We show that by repeatedly sampling nonphysical subtrees and truncating them at a given depth $\ell$, we obtain an efficient simulation algorithm.

**Theorem 4.** Let $U_\theta$ be a noisy quantum circuit as defined in Equation 2 and $P$ a Pauli observable with expectation value $f(\theta) = \text{tr}(U_\theta[|0\rangle\langle 0|]P)$. It is possible to compute an approximation $\tilde{f}$ of $f$ in time $O(Kn^2m2^\ell)$, where $\ell$ and $K$ are a chosen cut-off parameter and sampling overhead, such that with probability at least $1 - \delta$,

$$ \Delta(f, \tilde{f}) \le (1 - \gamma)^{(\ell+1)/2} + \sqrt{\frac{2 \log (\delta^{-1}/2)}{K}} \quad (C1) $$

where the probability comes from the random sampling of the trees.

*Proof.* Denote by $T_k$ a nonphysical tree created by sampling processes $D_{0z}$ with probability $1 - \gamma$ or $D_{0i}$ with probability $\gamma$ each time a Pauli $Z$ is encountered by the noisy rotation when backpropagating the observable. It is straightforward to see that the following estimator is unbiased:

$$ \hat{f}(\cdot) = \frac{1}{K} \sum_{k=1}^{K} \sum_{\omega \in T_k} \Phi_{\omega}(\cdot) d_{\omega}. \quad (C2) $$
**Algorithm 2:** [MC-LOWESA-AD] Monte-Carlo simulation of Pauli string observables under non-unital noise

**Input:** Quantum circuit affected by non-unital noise, and defined by process modes {$D_\omega$} and noise coefficient $\gamma$; measurement Pauli operator $P$; cut-off parameter $\ell$; sampling overhead $K$

**Output:** $\tilde{f}(\theta)$, an approximation of $f(\theta)$

**Procedure** MC-LOWESA-AD: Given a quantum circuit with rotation parameters $\theta$ and Pauli observable $P$, outputs an approximation of the expectation value of $P$, by randomly sampling $K$ trees and creating a surrogate of the energy landscape for all those trees. The angles are then adjusted.

1: $\tilde{f}(\theta) \leftarrow 0$

2: **for** $k = 1$ **to** $K$ **do**

3: Run *MC Low-weight coefficients* to calculate $d_\omega$ for all $\omega$ such that $\#\omega \le \ell$.

4: $\tilde{f}_k(\theta) \leftarrow 0$

5: **for all** $\omega$ such that $\#\omega \le \ell$ with non-zero $d_\omega$ **do**

6: $\tilde{f}_k(\theta) \leftarrow \tilde{f}_k(\theta) + d_\omega \Phi_\omega(\theta)$

7: **end for**

8: $\tilde{f}(\theta) \leftarrow \frac{1}{K} \tilde{f}_k(\theta)$

9: **end for**

10: Return $\tilde{f}(\theta)$

**end Procedure**

**Procedure** MC-LOWESA-AD: Given a quantum circuit and Pauli observable $P$, outputs a randomly sampled backpropagated tree and its coefficients $d_\omega$ for all $\omega$ such that $\#\omega \le \ell$ in time $O(n^2m2^\ell)$. The sampled tree is constructed as follows: at each layer of the quantum circuit, the tree $T_i$ is extended at each branch $\omega$ based on the Clifford layer $C_{m-i+1}$ and the Pauli operator on the rotation qubit $q_{m-i+1}$. If $P_\omega^{q_{m-i+1}} = Z$, only one of the two branches is kept, with probability $\gamma$ or $1-\gamma$, which does not increase the total number of branches in the tree. Note that the initial tree $T_0$ is made of a single node corresponding to the Pauli observable $P$. The branches of the tree that have split more than $\ell$ times are discarded.

1: **for** $i = 1$ **to** $m$ **do**

2: **for** $\omega$ in $T_{i-1}$ **do**

3: $P_\omega \leftarrow C_{m-i+1}^\dagger P_\omega C_{m-i+1}$, up to a phase $\phi_{\omega_i} \in \{\pm 1\}$ which is stored.

4: **if** $P_\omega^{q_{m-i+1}} = I$ **then**

5: Add branch $\omega$ in $T_i$ with $\omega_i = 0$

6: $d_{\omega_i}, P_\omega \leftarrow D_{\omega_i}^\dagger (P_\omega)$

7: **else if** $P_\omega^{q_{m-i+1}} = Z$ **then**

8: **with probability** $\gamma$ **do**

9: Add branch $\omega$ in $T_i$ with $\omega_i = 0_I$

10: $d_{\omega_i}, P_\omega \leftarrow D_{\omega_i}^\dagger (P_\omega)$

11: **or with probability** $1-\gamma$ **do**

12: Add branch $\omega'$ in $T_i$ with $\omega'_i = 0_Z$

13: $d_{\omega'_i}, P'_\omega \leftarrow D_{\omega'_i}^\dagger (P_\omega)$

14: **else if** $P_\omega^{q_{m-i+1}} \in \{X, Y\}$ and $\#\omega < \ell$ **then**

15: Split into two branches $\omega$ and $\omega'$ in $T_i$ with $\omega_i = +1$ and $\omega'_i = -1$

16: $d_{\omega_i}, P_\omega \leftarrow D_{\omega_i}^\dagger (P_\omega)$

17: $d_{\omega'_i}, P_{\omega'} \leftarrow D_{\omega'_i}^\dagger (P_\omega)$

18: $\#\omega, \#\omega' \leftarrow +1$

19: **else if** $P_\omega^{q_{m-i+1}} \in \{X, Y\}$ and $\#\omega = \ell$ **then**

20: Discard path $\omega$

21: **end if**

22: **end for**

23: **end for**

24: **for** $\omega$ in $T_m$ **do**

25: $P_\omega \leftarrow C_0^\dagger P_\omega C_0$ and phase $\phi_{\omega_0} \in \{\pm 1\}$

26: **end for**

27: Return $d_\omega = (\prod_{i=1}^m d_{\omega_i} \phi_{\omega_i}) \phi_{\omega_0} (\langle 0 | P_\omega \rangle)$ for all $\omega$ in $T_m$

**end Procedure**
Indeed, consider a split occurring through processes $D_{0z}$ and $D_{0i}$. This gives rise to two different branches, whose contributions to $f$ we denote by $(1-\gamma) \times S_1$ and $\gamma \times S_2$, where the damping prefactor comes from the processes. By randomly sampling only one of the branches with the probabilities mentioned, it is clear that the expectation value of this estimator is $(1-\gamma) \times S_1 + \gamma \times S_2$, making it unbiased.

Instead of considering the entire trees created by this Monte Carlo approach, which may still have an exponential number of branches, we propose only taking into account the branches of each tree that have split fewer than a constant number of times $\ell$. This approximation of our unbiased estimator is written as:

$$ \tilde{f}(\cdot) = \frac{1}{K} \sum_{k=1}^{K} \sum_{\substack{\omega \in T_k \\ \#\omega \le \ell}} \Phi_{\omega}(\cdot) d_{\omega} \quad (C3) $$

Note that for the trees $T_k$ considered, all splits come from processes $D_{\pm 1}$. The error we wish to bound is the expected $L^2$-error between our original parametrized expectation value $f$, and the one obtained by the truncated Monte Carlo approach $\tilde{f}$.

$$ \mathbb{E}[\Delta(f, \tilde{f})] = \mathbb{E}[\|f - \tilde{f}\|_2] = \mathbb{E}\left[\left(\frac{1}{|\Theta|} \int_{\Theta} |f(\theta) - \tilde{f}(\theta)|^2 d\theta\right)^{1/2}\right] \quad (C4) $$

This error can be upper bounded by two distinct errors using the triangle inequality: one corresponding to the truncation error and the other to the Monte Carlo approximation, $\|f - \tilde{f}\|_2 = \|(f - \hat{f}) + (\hat{f} - \tilde{f})\|_2 \le \|f - \hat{f}\|_2 + \|\hat{f} - \tilde{f}\|_2$.

$$ \mathbb{E}[\Delta(f, \tilde{f})] \le \mathbb{E}[\|f - \hat{f}\|_2] + \mathbb{E}[\|\hat{f} - \tilde{f}\|_2] \quad (C5) $$

The first term comes from the mean square error of our Monte Carlo simulation, and is bounded by $1/K$, as our estimator is unbiased and the sampled trees yield expectation values in the range of $[-1, 1]$. Computing the second term, the error made by truncating each of our simulated trees, requires additional effort.

$$ \|\hat{f} - \tilde{f}\|_2^2 = \frac{1}{|\Theta|} \int_{\Theta} |\hat{f}(\theta) - \tilde{f}(\theta)|^2 d\theta = \frac{1}{K^2} \sum_{k=1}^{K} \sum_{k'=1}^{K} \frac{1}{|\Theta|} \int_{\Theta} \sum_{\substack{\omega \in T_k \\ \#\omega > \ell}} \Phi_{\omega}(\theta) d_{\omega} \sum_{\substack{\omega' \in T_{k'} \\ \#\omega' > \ell}} \Phi_{\omega'}(\theta) d_{\omega'} d\theta \quad (C6) $$

Using the Cauchy-Schwarz inequality, Lemma 2 (since the only splits in the trees are from $D_{\pm 1}$, the condition $h(\omega) = h(\omega')$ can be reduced to $\omega = \omega'$) and the fact that all splits in a tree $T_k$ come from processes $D_{\pm 1}$, the error term can be rewritten:

$$ \begin{aligned} \|\hat{f} - \tilde{f}\|_2^2 &\le \frac{1}{K^2} \sum_{k=1}^K \sum_{k'=1}^K \left\| \sum_{\substack{\omega \in T_k \\ \#\omega > \ell}} \Phi_{\omega} d_{\omega} \right\|_2 \left\| \sum_{\substack{\omega \in T_{k'} \\ \#\omega > \ell}} \Phi_{\omega} d_{\omega} \right\|_2 \\ &\le \frac{1}{K^2} \left( \sum_{k=1}^K \left( \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} d_{\omega} d_{\omega'} \right)^{1/2} \right)^2 \end{aligned} \quad (C7) $$

Once again, it is possible to rewrite $d_{\omega} = Q(\omega)d_{\omega}^0$ where $Q(\omega)$ represents the noise accumulated through path $\omega$, with $Q(\omega) = \prod_i Q(\omega_i)$ with $Q(+1) = Q(-1) = \sqrt{1-\gamma}$ and $Q(0) = 1$. Since $d_{\omega}^0$ is still the expectation value of a single Pauli string, we have $|d_{\omega}^0| \le 1$. The paths having split through the $D_{\pm 1}$ processes at least $\ell + 1$ times,

$$ \|\hat{f} - \tilde{f}\|_2^2 \le \frac{1}{K^2} \left( \sum_{k=1}^{K} \left( \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} Q(\omega)^2 |d_{\omega}^0|^2 \right)^{1/2} \right)^2 \le \frac{(1-\gamma)^{\ell+1}}{K^2} \left( \sum_{k=1}^{K} \left( \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} \right)^{1/2} \right)^2 \quad (C8) $$

Finally, it remains to notice that for any given tree $T_k$,
$$ \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} \le \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \omega = \omega'}} 2^{-\#\omega} = 1 \qquad (C9) $$

This allows us to bound the truncation error and shows,

$$ \mathbb{E}[\Delta(f, \tilde{f})] \le (1 - \gamma)^{(\ell+1)/2} + 1/\sqrt{K} \qquad (C10) $$

By applying Hoeffding's inequality and noting that each of the sampled trees yields an expectation value in the range of $[-1, 1]$, we can conclude that with probability at least $1 - \delta$,

$$ \Delta(f, \tilde{f}) \le (1 - \gamma)^{(\ell+1)/2} + \sqrt{\frac{2 \log (\delta^{-1}/2)}{K}} \qquad (C11) $$

□

## Appendix D: Extension to normal form non-unital noise channel

The results previously obtained for the amplitude damping channel can be extended to any non-unital noise that can be brought into the normal form representation of the noise channel [15, 18, 35] by acting on the left or right with single qubit unitaries. More precisely, recall that for any single-qubit noise channel $\mathcal{N}$, there exist unitaries $U, V$ with the property that $\mathcal{U}\mathcal{N}\mathcal{V}$ can be represented by two vectors $\mathbf{t} = (t_X, t_Y, t_Z) \in \mathbb{R}^3$ and $\mathbf{D} = (D_X, D_Y, D_Z) \in \mathbb{R}^3$ such that its action on a vector of the Bloch sphere can be written,

$$ \mathcal{N} \left( \frac{\mathbf{I} + \mathbf{w} \cdot \boldsymbol{\sigma}}{2} \right) = \frac{\mathbf{I}}{2} + \frac{1}{2} (\mathbf{t} + D\mathbf{w}) \cdot \boldsymbol{\sigma} \qquad (D1) $$

where we have defined $D := \mathrm{diag}(\mathbf{D})$, $\sigma := \{X, Y, Z\}$ the vector of single-qubit Pauli matrices and $\mathbf{w} \in \mathbb{R}^3$ such that $||\mathbf{w}||_2 \le 1$. The adjoint channel $\mathcal{N}^\dagger$ acts on the Pauli basis in the following way,

$$ \forall P \in \{X, Y, Z\}, \quad \mathcal{N}^\dagger(P) = t_P I + D_P P \qquad (D2) $$

We will now consider channels that can be brought into this form where the unitaries $U, V$ are Clifford. Note that the special cases of depolarizing, dephasing and amplitude damping channels are members of this class with $U, V = I$. Furthermore, note that any composition or convex combination of such channels will remain in the class. As such, this class of channels cover probabilistic and composite combinations of the majority of single qubit channels considered in the literature. We consider the case where the noise channel is non-unital, which implies $\mathbf{t} \ne (0, 0, 0)$ [15]. Note that unital noise is already known to be simulable by similar frequency truncation techniques [19]. We will be making extensive use of the following properties of the channel in its normal form.

**Lemma 5** (Constraints on the parameters of the normal form, adapted from Ref. [30]). Let $\mathbf{D}, \mathbf{t}$ be the parameters of a non-unital channel in normal form. Then, for any $\mathbf{b} = (b_X, b_Y, b_Z) \in \mathbb{R}^3$ such that $||\mathbf{b}||_2 = 1$,

$$ \sum_{P \in \{X, Y, Z\}} b_P^2 |D_P| + \sum_{P \in \{X, Y, Z\}} |b_P t_P| \le 1 \qquad (D3) $$

In particular we have the two following properties:

1. $\forall P \in \{X, Y, Z\}, \quad |D_P| + |t_P| \le 1,$

2. and there exists at most one $P \in \{X, Y, Z\}$ such that $|D_P| + |t_P| = 1$.

*Proof.* Let $O$ be the observable and $\rho$ the state defined by,

$$ O = \sum_{P \in \{X, Y, Z\}} |b_P| \cdot \mathrm{sign}(t_P D_P) P \quad \text{and} \quad \rho = \frac{I + O}{2} \qquad (D4) $$

The operator norm of $O$ can be bounded as follows,
$$ ||O||_{\infty} = \max_{\sigma} |\text{tr}(O\sigma)| = \max_{\substack{\mathbf{r} \in \mathbb{R}^3 \\ ||\mathbf{r}||_2 = 1}} \left| \sum_{P \in \{X,Y,Z\}} |b_P| \cdot r_P \text{sign}(t_P) \right| \le \max_{\substack{\mathbf{r} \in \mathbb{R}^3 \\ ||\mathbf{r}||_2 = 1}} ||\mathbf{b}||_2 ||\mathbf{r}||_2 = 1 \quad (D5) $$

Where the last inequality is obtained through Cauchy-Schwarz. Equality is attained by choosing $\sigma = \rho$, such that $||O||_{\infty} = 1$. Furthermore, Hölder's inequality allows us to bound,

$$ 1 = ||O||_{\infty} ||\rho||_1 \ge |\text{tr}[\mathcal{O}\mathcal{N}(\rho)]| \quad (D6) $$

$$ = \frac{1}{2} \left| \text{tr} \left[ \left( \sum_{P \in \{X,Y,Z\}} |b_P| \cdot \text{sign}(t_P D_P) P \right) \mathcal{N} \left( I + \sum_{P \in \{X,Y,Z\}} |b_P| \cdot \text{sign}(t_P D_P) P \right) \right] \right| \quad (D7) $$

$$ = \frac{1}{2} \left| \text{tr} \left[ \left( \sum_{P \in \{X,Y,Z\}} |b_P| \cdot \text{sign}(t_P D_P) P \right) \left( I + \sum_{P \in \{X,Y,Z\}} (|b_P| D_P \cdot \text{sign}(t_P D_P) + t_P) P \right) \right] \right| \quad (D8) $$

$$ = \left| \sum_{P \in \{X,Y,Z\}} b_P^2 D_P + |b_P t_P| \text{sign}(D_P) \right| = \sum_{P \in \{X,Y,Z\}} b_P^2 |D_P| + \sum_{P \in \{X,Y,Z\}} |b_P t_P|, \quad (D9) $$

where in the last step we used the fact that the coefficients $D_P$ have all the same sign. Obtaining the first property mentioned is then done by setting $\mathbf{b}$ to either $(1,0,0)$, $(0,1,0)$ or $(0,0,1)$, which shows that $\forall P \in \{X,Y,Z\}$, $|D_P| + |t_P| \le 1$. Suppose now that for one of the Paulis $P$, $|D_P| + |t_P| = 1$. Then there cannot be another $P' \ne P$ such that $|D_{P'}| + |t_{P'}| = 1$. Indeed, Equation D3 would imply for any $0 < b_P < 1$ and $0 < b_{P'} < 1$ with $b_P^2 + b_{P'}^2 = 1$ that,

$$ b_{P'}^2 D_{P'} + b_P^2 D_P + b_{P'}|t_{P'}| + b_P|t_P| \le 1 \quad (D10) $$

However, this is not the case if both $|D_{P'}| + |t_{P'}| = 1$ and $|D_P| + |t_P| = 1$ and either $t_{P'} \ne 0$ or $t_P \ne 0$,

$$ \begin{aligned} b_{P'}^2 D_{P'} + b_P^2 D_P + b_{P'}|t_{P'}| + b_P|t_P| &= b_{P'}^2(1 - |t_{P'}|) + b_P^2(1 - |t_P|) + b_Z|t_Z| + b_P|t_P| \\ &= 1 + b_{P'}|t_{P'}|(1 - b_{P'}) + b_P|t_P|(1 - b_P) \\ &> 1, \end{aligned} \quad (D11) $$

where in the second step we used the fact that $b_{P'}^2 + b_P^2 = 1$. $\square$

Note that if there is no $P$ such that $|D_P| + |t_P| = 1$, then it is possible to simulate circuits with rotations along any axis, making our simulation method even more powerful. In the following, if the noise channel is such that there exists a $P \in \{X,Y,Z\}$ satisfying $|D_P| + |t_P| = 1$, then we will consider $P = Z$, similarly to the amplitude damping case. This implies that $|D_X| + |t_X| < 1$ and $|D_Y| + |t_Y| < 1$. Going back to the circuit architecture considered, our noisy circuit is written,

$$ U_{\theta} = \left( \bigcirc_{i=1}^{m} C_i \circ \mathcal{R}_z^{(q_i)}(\theta_i) \circ \mathcal{N} \right) \circ C_0 \quad (D12) $$

We will now consider that the Cliffords putting in diagonal form are identity without loss of generality, as the additional diagonalizing Cliffords can be absorbed into the circuit. With the noise in its normal form, the noisy rotation can be written as

$$ \begin{aligned} \mathbf{R}_z \cdot \mathbf{N} &= \begin{bmatrix} 1 & t_X & t_Y & t_Z \\ 0 & D_X \cos \theta & -D_Y \sin \theta & 0 \\ 0 & D_X \sin \theta & D_Y \cos \theta & 0 \\ 0 & 0 & 0 & D_Z \end{bmatrix} \\ &= \mathbf{D}_0 + \sum_{P \in \{X,Y,Z\}} \mathbf{D}_{0_P} + \sum_{P \in \{X,Y\}} (\mathbf{D}_{1_P} \cos \theta + \mathbf{D}_{-1_P} \sin \theta) + \mathbf{D}_{0_I} \end{aligned} \quad (D13) $$

where we have defined the quantum processes $\mathbf{D}_0 = |I\rangle\langle I|$, $\mathbf{D}_{0_P} = t_P|I\rangle\langle P|$, $\mathbf{D}_{0_I} = D_Z|Z\rangle\langle Z|$, $\mathbf{D}_{1_P} = D_P|P\rangle\langle P|$, $\mathbf{D}_{-1_X} = D_X|Y\rangle\langle X|$, and $\mathbf{D}_{-1_Y} = -D_Y|X\rangle\langle Y|$, with $|I\rangle = [1 \ 0 \ 0 \ 0]^T$, $|X\rangle = [0 \ 1 \ 0 \ 0]^T \dots$.
When backpropagating our Pauli string, the noisy rotation acts only on a single qubit $q_i$. Depending on the Pauli operator at qubit $q_i$, at most 3 processes are valid: $D_0$ for $I$, $D_{0z}$ and $D_{0i}$ for $Z$, $D_{0x}$, $D_{1x}$ and $D_{-1x}$ for $X$ and $D_{0y}$, $D_{1y}$ and $D_{-1y}$ for $Y$. This once again suggests the structure of a rooted tree $T = (V, E)$, where the nodes $V$ represent the backpropagated Pauli strings at each layer, and the edges represent the valid processes $D_{\omega_i}$ for the given Paulis. Similarly as in the amplitude damping case, keeping track of the rotation monomials and the splitting of the paths allows us to write our expectation value as

$$f(\theta) = \mathrm{tr}(U_{\theta}[|0\rangle\langle0|]P) = \sum_{\omega \in T} d_{\omega} \Phi_{\omega}(\theta) \qquad (\mathrm{D}14)$$

However, each branch can split into up to three different branches, which means that recovering $f$ exactly would require time $O(n^2m3^m)$ in the worst case. To handle the additional branches coming from the non-unital noise, we once again propose to randomly sample only some of them. More particularly, our Monte Carlo approach creates nonphysical trees as follows:

* If $P_{q_i} = X$ or $Y$, keep process $D_{0P}$ with probability $|t_P|/(|D_P| + |t_P|)$ and both processes $D_{\pm 1P}$ with probability $|D_P|/(|D_P| + |t_P|)$.

* If $P_{q_i} = Z$, keep process $D_{0z}$ with probability $|t_Z|/(|D_Z| + |t_Z|)$ and process $D_{0i}$ with probability $|D_Z|/(|D_Z| + |t_Z|)$.

Note that keeping the estimator unbiased requires us to multiply the branches kept by $(|D_P| + |t_P|)$ at each split. For $P = X$ or $Y$, this effectively contracts the expectation value since for $P \in \{X, Y\}$, $(|D_P| + |t_P|) < 1$. Therefore, truncating the sampled trees after a constant amount of splits $\ell$ through $D_{\pm 1P}$ yields a good approximation.

**Theorem 6.** Let $U_{\theta}$ be a noisy quantum circuit as defined in Equation D12 and $P$ a Pauli observable with expectation value $f(\theta) = \mathrm{tr}(U_{\theta}[|0\rangle\langle0|]P)$. Denote by $|D| + |t| = \max_{P \in \{X, Y\}} |D_P| + |t_P| < 1$. It is possible to compute an approximation $\tilde{f}$ of $f$ in time $O(Kn^2m2^\ell)$, where $\ell$ and $K$ are a chosen cut-off parameter and sampling overhead, such that with probability at least $1 - \delta$,

$$\Delta(f, \hat{f}) \le (|D| + |t|)^{(\ell+1)} + \sqrt{\frac{2 \log (\delta^{-1}/2)}{K}} \qquad (\mathrm{D}15)$$

where the probability comes from the random sampling of the trees.

*Proof.* Denote by $T_k$ a nonphysical tree created by sampling processes $D_{0P}$ with probability $|t_P|/(|D_P| + |t_P|)$ or keeping the other branch(s) (both $D_{\pm 1P}$ if $P \in \{X, Y\}$ or $D_{0i}$ if $P = Z$) with probability $|D_P|/(|D_P| + |t_P|)$ each time a Pauli $P$ is encountered by the noisy rotation when backpropagating the observable. It is straightforward to see that the following estimator is unbiased:

$$\hat{f}(\cdot) = \frac{1}{K} \sum_{k=1}^{K} \sum_{\omega \in T_k} \Phi_{\omega}(\cdot) d_{\omega}. \qquad (\mathrm{D}16)$$

Indeed, consider a branch that splits by encountering a Pauli $P \in \{X, Y\}$ during the backpropagation. This split leads to 3 different branches through processes $D_{0P}$, $D_{1P}$ and $D_{-1P}$, which contribution to the expectation value we can write $D_P \cos \theta S_1 + D_P \sin \theta S_2 + t_P S_3$, where the noise prefactors come from the processes. By randomly sampling the branches with the probabilities mentioned before and multiplying their contributions by the factor $|D_P| + |t_P|$, the expectation value of our estimator is $|t_P|/(|D_P| + |t_P|) \times S_3 \times (|D_P| + |t_P|) + |D_P|/(|D_P| + |t_P|) \times (\cos \theta S_1 + \sin \theta S_2) \times (|D_P| + |t_P|)$. To make this truly unbiased, one also needs to adjust for the sign of $t_P$ or $D_P$ when multiplying by $(|D_P| + |t_P|)$. Similarly, if the split comes from processes $D_{0z}/D_{0i}$, the two resulting branches have contribution $D_Z \times S_1 + t_Z \times S_2$. Sampling either $D_{0z}$ with probability $|t_Z|/(|D_Z| + |t_Z|)$ or $D_{0z}$ with probability $|D_Z|/(|D_Z| + |t_Z|)$ and multiplying the sampled branch by $|D_Z| + |t_Z|$ yields an expectation value of $|D_Z|/(|D_Z| + |t_Z|) \times S_1 \times (|D_Z| + |t_Z|) + |t_Z|/(|D_Z| + |t_Z|) \times S_2 \times (|D_Z| + |t_Z|)$. Adjusting once again for the sign of $t_Z$ and $D_Z$ suffices to make our estimator unbiased.

The trees $T_k$ created now only contain splits that occur through processes $D_{\pm 1x}$ or $D_{\pm 1y}$. There may still be an exponential amount of branches in the tree though, so we complement this Monte-Carlo approach with the same truncation rule as before, and only keep branches that have split less than a fixed number of times $\ell$. At each of these splits we show that the error made by our approximation is dampen by $|D_X| + |t_X| < 1$ or $|D_Y| + |t_Y| < 1$.
This approximation of our unbiased estimator is written as:

$$ \tilde{f}(\cdot) = \frac{1}{K} \sum_{k=1}^{K} \sum_{\substack{\omega \in T_k \\ \#\omega \le \ell}} \Phi_{\omega}(\cdot) d_{\omega} \quad (D17) $$

The error we wish to bound is once again the $L^2$-error between our original parametrized expectation value $f$, and the one obtained by the truncated Monte Carlo approach $\tilde{f}$. This error is still upper bounded by two distinct errors using the triangle inequality: one corresponding to the truncation error and the other to the Monte Carlo approximation, $||f - \tilde{f}||_2 = ||(f - \hat{f}) + (\hat{f} - \tilde{f})||_2 \le ||f - \hat{f}||_2 + ||\hat{f} - \tilde{f}||_2$.

The first term comes from the mean square error of our Monte Carlo simulation, and can be bounded using Hoeffding's inequality, as our estimator is unbiased. Computing the second term, the error made by truncating can be done in a similar fashion as for amplitude damping.

$$ ||\hat{f} - \tilde{f}||_2^2 = \frac{1}{|\Theta|} \int_{\Theta} |\hat{f}(\theta) - \tilde{f}(\theta)|^2 d\theta = \frac{1}{K^2} \sum_{k=1}^{K} \sum_{k'=1}^{K} \frac{1}{|\Theta|} \int_{\Theta} \sum_{\substack{\omega \in T_k \\ \#\omega > \ell}} \Phi_{\omega}(\theta) d_{\omega} \sum_{\substack{\omega' \in T_{k'} \\ \#\omega' > \ell}} \Phi_{\omega'}(\theta) d_{\omega'} \quad (D18) $$

Using the Cauchy-Schwarz inequality, Lemma 2 (since the only splits in the trees are from $D_{\pm 1_x}$ or $D_{\pm 1_y}$, the condition $h(\omega) = h(\omega')$ can be reduced to $\omega = \omega'$) and the fact that all splits in a tree $T_k$ come from processes $D_{\pm 1_x}$ or $D_{\pm 1_y}$, the error term can be rewritten:

$$ \begin{aligned} ||\hat{f} - \tilde{f}||_2^2 &\le \frac{1}{K^2} \sum_{k=1}^{K} \sum_{k'=1}^{K} \left\| \sum_{\substack{\omega \in T_k \\ \#\omega > \ell}} \Phi_{\omega} d_{\omega} \right\|_2 \left\| \sum_{\substack{\omega' \in T_{k'} \\ \#\omega' > \ell}} \Phi_{\omega'} d_{\omega'} \right\|_2 \\ &\le \frac{1}{K^2} \left( \sum_{k=1}^{K} \left( \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} d_{\omega} d_{\omega'} \right)^{1/2} \right)^2 \end{aligned} \quad (D19) $$

Since we have sampled the branches with probability $|t_P|/(|D_P|+|t_P|)$ and $|D_P|/(|D_P|+|t_P|)$, keeping the estimator unbiased required us to multiply each branch by $(|D_P| + |t_P|)$, and adjust for the sign of $t_P$ or $D_P$. It is therefore possible to rewrite for all $\omega \in T_k$, $d_\omega = Q(\omega)d_\omega^0$ where $Q(\pm 1_X) = Q(0_X) = |D_X| + |t_X| < 1$, $Q(\pm 1_Y) = Q(0_Y) = |D_Y| + |t_Y| < 1$, $Q(0_I) = Q(0_Z) = |D_Z| + |t_Z| \le 1$ and $d_\omega^0 \le 1$. Denote by $|D| + |t| = \max_{P \in \{X, Y\}} |D_P| + |t_P|$. Since each path in the sum has split at least $\ell + 1$ times through processes $D_{\pm 1_x}$ or $D_{\pm 1_y}$, we have that,

$$ ||\hat{f} - \tilde{f}||_2^2 \le \frac{1}{K^2} \left( \sum_{k=1}^{K} \left( \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} Q(\omega)^2 |d_\omega^0|^2 \right)^{1/2} \right)^2 \le \frac{(|D| + |t|)^{2(\ell+1)}}{K^2} \left( \sum_{k=1}^{K} \left( \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} \right)^{1/2} \right)^2 \quad (D20) $$

Since the sampled trees contain only splits in up to two branches, and branches only split when encountering processes $D_{\pm 1_x}$ or $D_{\pm 1_y}$, we still have that for any given tree $T_k$,

$$ \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \#\omega > \ell, \#\omega' > \ell \\ \omega = \omega'}} 2^{-\#\omega} \le \sum_{\substack{(\omega, \omega') \in T_k^2 \\ \omega = \omega'}} 2^{-\#\omega} = 1 \quad (D21) $$

This allows us to bound the truncation error,

$$ ||\tilde{f} - \hat{f}||_2 \le (|D| + |t|)^{\ell+1} \quad (D22) $$

Using Hoeffding's inequality and that each of the sampled trees yields an expectation value in the range of $[-1, 1]$, we can conclude that with probability at least $1 - \delta$,

$$ \Delta(f, \hat{f}) \le (|D| + |t|)^{\ell+1} + \sqrt{\frac{2 \log (\delta^{-1}/2)}{K}} \quad (D23) $$

□
## Appendix E: Comparison to other simulation methods

**Comparison to other results:** Pauli backpropagation has proven to be a powerful technique for simulating quantum circuits, both in noiseless [20, 22, 24] and noisy [14, 15, 19, 23, 28] settings. These methods are motivated by theoretical guarantees that lay the foundation for their practical use, which are summarized in Table I. We see that our guarantees have different improvements over existing results in the literature: either they hold for larger classes of circuits, require less randomness or have better runtime guarantees.

For unital noise, such as depolarizing and dephasing noise, [19] and [23] have demonstrated that it is possible to approximate expectation values in average over rotation angles in polynomial time. These algorithms are both space and time efficient, enabling fast computations of expectation values. Similar results were reported in [14], where the approximation holds on average over the input states. Additionally, assuming anti-concentration, accurate sampling from the quantum circuit was shown to be feasible.

However, extending these methods beyond unital noise remained challenging. For instance, [14] investigated only a randomized version of amplitude damping, while [23] had to include depolarizing noise on top of the amplitude damping channel to ensure their approximation held. For random circuits, a quasi-polynomial time Pauli backpropagation algorithm was derived (assuming geometrically local circuits) under non-unital noise [15]. However, under slightly stronger randomness assumptions, it has been shown that even noiseless circuits can be simulated using Pauli backpropagation techniques [22], with slightly worse run-time. In the noiseless regime, recent work has also established approximation bounds on certain regions of the energy landscape [24].

<table><caption>TABLE I: Comparison of the different Pauli backpropagation methods and their guarantees. The runtime may depend on the number of qubits *n*, the depth of the quantum circuit *m*, the circuit geometry *D*, the additive precision *ε* and the noise strength *p*.</caption><thead><tr><th>Reference</th><th>Ref. [19, 23]</th><th>Ref. [14]</th><th>Ref. [15]</th><th>Ref. [22]</th><th>Ref. [24]</th><th>Alg. 1</th><th>Alg. 2</th></tr></thead><tbody><tr><td>Noise model</td><td>Unital noise</td><td>Depolarising +<br>random amplitude damping</td><td>Non-unital noise</td><td>Noiseless</td><td>Noiseless</td><td>Non-unital noise</td><td>Non-unital noise</td></tr><tr><td>Circuit</td><td>Any circuit,<br>random angles</td><td>Any circuit,<br>random input states</td><td>Random circuits</td><td>Random circuits</td><td>Any circuit,<br>random small angles</td><td>Almost any circuit,<br>random angles</td><td>Any circuit,<br>random angles</td></tr><tr><td>Runtime<br>(ε, p, n, m, D)</td><td>n<sup>2</sup>m(1/ε)<sup>1/p</sup></td><td>mn<sup>log(&radic;m+1/ε)/p</sup></td><td>e<sup>log<sup>D</sup>(ε<sup>-1</sup>)</sup></td><td>mn<sup>log 1/ε</sup></td><td>m<sup>log 1/ε</sup></td><td>n<sup>2</sup>m(1/ε<sup>2</sup>)<sup>3/2p</sup></td><td>n<sup>2</sup>m(1/ε<sup>2</sup>)<sup>1+1/p</sup></td></tr></tbody></table>

**Future work:** A significant open question in the field concerns the extension of Pauli backpropagation methods to accommodate any noise model. In this work, we have made a great progress in this direction by handling a broad class of non-unital noise, but there still exist noise channels that fall beyond the scope of this work. Adapting to even more general noise models demands additional efforts and poses challenges due to the complexity of the backpropagated tree splitting. Furthermore, we have shown that for many noise models, it is possible to handle rotations along *any* axis. It still remains to show whether this is true for all non-unital noise models. Another important problem is the simulation of noisy continuous time dynamics with noise, which could offer valuable insights into the performance of noisy quantum simulators. Further research in these areas could help understand the conditions under which quantum computations could offer real advantages.