# Operational Markovianization in Randomized Benchmarking

P. Figueroa-Romero<EMAIL_ADDRESS>IQM Germany GmbH,
Nymphenburgerstrasse 86, 80636 Munich, Germany  M. Papič  IQM Germany GmbH,
Nymphenburgerstrasse 86, 80636 Munich, Germany  Department of Physics and
Arnold Sommerfeld Center for Theoretical Physics, Ludwig-Maximilians-
Universität München, Theresienstrasse 37, 80333 Munich, Germany  A. Auer  IQM
Germany GmbH, Nymphenburgerstrasse 86, 80636 Munich, Germany  M.-H. Hsieh  Hon
Hai Quantum Computing Research Center, Taipei, Taiwan  K. Modi  School of
Physics and Astronomy, Monash University, VIC 3800, Australia  Centre for
Quantum Technology, Transport for New South Wales, Sydney, NSW 2000, Australia
I. de Vega  IQM Germany GmbH, Nymphenburgerstrasse 86, 80636 Munich, Germany
Department of Physics and Arnold Sommerfeld Center for Theoretical Physics,
Ludwig-Maximilians-Universität München, Theresienstrasse 37, 80333 Munich,
Germany

###### Abstract

A crucial task to obtain optimal and reliable quantum devices is to quantify
their overall performance. The average fidelity of quantum gates is a
particular figure of merit that can be estimated efficiently by Randomized
Benchmarking (RB). However, the concept of gate-fidelity itself relies on the
crucial assumption that noise behaves in a predictable, time-local, or so-
called Markovian manner, whose breakdown can naturally become the leading
source of errors as quantum devices scale in size and depth. We analytically
show that error suppression techniques such as Dynamical Decoupling (DD) and
Randomized Compiling (RC) can operationally _Markovianize_ RB: _i_) fast DD
reduces non-Markovian RB to an exponential decay plus longer-time corrections,
while on the other hand, _ii_) RC generally does not affect the average, but
_iii_) it always suppresses the variance of such RB outputs. We demonstrate
these effects numerically with a qubit noise model. Our results show that
simple and efficient error suppression methods can simultaneously tame non-
Markovian noise and allow for standard and reliable gate quality estimation, a
fundamentally important task in the path toward fully functional quantum
devices.

##  I Introduction

The characterization of noise in quantum information processors will remain a
necessary and unavoidable task to build fault-tolerant devices. Some of the
most common sets of techniques to achieve this can be said to fall within an
interval comprising tomographic techniques on one end, and Randomized
Benchmarking (RB) techniques on the other: while the first can provide a
detailed description of noise with an exponential measurement and sampling
overhead in system size, the latter can estimate coarse average error rates
efficiently [1, 2].

Due to the simplicity of RB-based protocols, they have become ubiquitous both
for small systems [3, 4, 5, 6, 7, 8], as well as a steppingstone for scalable
techniques for more ambitious learning of quantum noise [9, 10, 11, 12, 13].
However, the manageable analytical behavior of the RB data –namely, an
exponential decay capturing average gate-fidelities, with State Preparation
and Measurement (SPAM) errors isolated as multiplicative and offset constants–
is guaranteed only with highly-simplified and unrealistic assumptions about
the noise. More realistic regimes have actively been under investigation [14,
15] and can still benefit from RB’s simplicity.

Arguably, however, one of the most difficult simplifications to relax in any
characterization technique is the one assuming that noise can be associated
with each individual quantum gate separately. In reality, noise is e.g.,
introduced by the nature of the quantum device itself or by inherent limits in
the control that the experimenter has, and errors occurring at a given time
can propagate and affect future errors. These kinds of correlations and their
effects are generally known as non-Markovianity [16, 17, 18] and can be fully
described in the quantum realm through multi-time generalizations of quantum
channels [19, 20, 21]. A signature aspect of non-Markovian noise in RB is the
non-trivial non-exponential decay of the data rendered by it [22, 23, 24, 25,
26], which makes the information about the noise and its correlations
remarkably hard to analyze. Moreover, the key feature of RB of robustness
against SPAM errors ceases to hold if any of these are correlated. This is
depicted schematically in Figure 1.

Figure 1: Randomized Benchmarking (RB) under Markovian and non-Markovian
noise: Cartoons of sample circuits of the RB protocol for a given initial
state $\rho_{\mathsf{S}}$, sequence of randomly sampled quantum gates
$\\{\mathcal{G}_{i}\\}$ with “undo” (compiled inverse of the sequence) gate
$\mathcal{G}_{m+1}$, and final measurement $M$, subject to (a) Markovian
noise, where errors are uncorrelated with each other and can be associated to
each individual gate, and (b) non-Markovian noise, where an external (quantum)
system –i.e., an environment– can serve as a memory, correlating errors in
time and propagating their information to the final measurement. Example
outputs of RB associated with (a) and (b) are shown on the top and bottom
right, correspondingly; detail in notation and meaning will be explained
throughout the manuscript, here we only point out the distinction that non-
Markovian RB decays are generally non-trivially non-exponential, with the
standard notion of an individual average gate-fidelity being ill-defined.

While characterizing and better understanding non-Markovian noise is an active
area of research [27, 28, 29, 30, 31, 32, 33], various control and error
suppression techniques, such as Dynamical Decoupling (DD) [34] and Randomized
Compiling (RC) [35] have been shown to remove non-Markovian noise effects to a
statistically significant extent. This can further be viewed with a resource-
theoretic lens [36] as keeping local-system information in exchange of the
consumption of temporal correlations. Generally, in the Markovian regime, both
DD and RC are well-known to enhance the quality of quantum computations with a
low resource overhead [37, 38, 39, 40], and combining RB with quantum control
has previously been done successfully, e.g., to optimize the quantum control
itself [41], or to demonstrate enhanced average gate-fidelities [42].

Main Results (informal): In this manuscript, we show that for a broad class of
non-Markovian noise models, DD can effectively and efficiently _Markovianize_
RB, i.e., remove non-Markovian non-exponential deviations, allowing for a
straightforward prediction of enhanced average gate-fidelities with low
overhead. In other words, DD converts non-Markovian correlations into more
tractable quantum noise. We exemplify this effect numerically on a qubit with
an $XY4$ sequence. Moreover, we analyze the effect of tailoring non-Markovian
noise into Pauli noise within the accessible subsystem of interest in non-
Markovian RB, finding that, while RC does not Markovianize RB in the same
sense that DD does, the uncertainty in the outputs get suppressed so that
average outputs can be almost as precise as analytical ones. We exemplify
these effects with the same numerical model for a single qubit strongly
interacting with another.

Our results show that coherent noise suppression and decoupling schemes in RB
can both efficiently Markovianize and allow the extraction of enhanced average
gate fidelities that accurately capture –and do not overestimate– all
statistically-relevant error rates.

##  II RB and operational modeling of quantum noise

###  II.1 Setup and notation

Detail on the notation we employ in the main text, as well as in the
derivations for our results can be read in full in Appendix A.1. Here we will
consider a composite Hilbert space labeled $\mathsf{S}\mathsf{E}$, comprised
by a system of interest $\mathsf{S}$, and an environment $\mathsf{E}$, which
we assume to be inaccessible. $\mathsf{E}$ could simply be a subset of idle or
ancillary qubits within a quantum device [26]. We will label their respective
dimensions as $d_{\mathsf{S}}=2^{n_{s}}$, for $n_{s}$ qubits in $\mathsf{S}$,
and $d_{\mathsf{E}}$.

We present our results with the uniformly distributed $n_{s}$-qubit Clifford
group as our gate-set and denote individual gates with the symbol
$\mathcal{G}_{i}$, $i=1,2,\ldots$, although the general treatment with finite
groups can be seen in the Appendices. We denote initial states as
$\rho_{\mathsf{S}}$ and measurements as $M$, both on system $\mathsf{S}$, and
which can be chosen arbitrarily, as long as the noiseless expectation
$\operatorname{tr}(M\rho_{\mathsf{S}})$ is known. While the RB protocol can be
read in detail in Appendix A.3, the three main input components in a RB
experiment are _i_) the gate-set $\\{\mathcal{G}_{i}\\}_{i=1}^{m}$ for a given
integer $m$, _ii_) the initial state $\rho_{\mathsf{S}}$, and _iii_) the
measurement element or observable $M$.

We will denote by $\mathbf{E}_{\mathcal{G}}$ the uniform average over the
gate-set $\\{\mathcal{G}_{i}\\}$, and reserve the notation $\mathcal{F}_{m}$
for a so-called average-sequence fidelity, defined here as
$\mathcal{F}_{m}:=\mathbf{E}_{\mathcal{G}}\operatorname{tr}[M^{\prime}\mathcal{G}_{m+1}^{\prime}\circ\mathcal{G}^{\prime}_{m}\circ\cdots\mathcal{G}^{\prime}_{1}(\rho_{\mathsf{S}}^{\prime})]$,
where the primed terms denote the real (noisy) implementations of the
corresponding initial state, measurement, and gates, and where
$\mathcal{G}_{m+1}:=\mathcal{G}_{m}^{-1}\circ\cdots\mathcal{G}_{1}^{-1}$.

Finally, we refer to a quantum channel being a “$\mathsf{X}$ channel”, or as a
“channel acting on $\mathsf{X}$”, when it maps inputs from space $\mathsf{X}$
to outputs in $\mathsf{X}$.

###  II.2 RB: Markovian vs. non-Markovian

Up to a few more parameters, like the number of gates sampled or the number of
samples to generate, a full-fledged framework for RB exists [15] when the
underlying noise is assumed to be effectively uncorrelated in time, i.e.,
Markovian. This allows us to estimate average gate fidelities, which despite
being a rough figure of merit, are an essential component for the
characterization of noise in a quantum device. In the Markovian, time-
stationary, and gate-independent noise regime, the outputs of a RB experiment
over gate sequences of length $m$ are estimates that can be fitted to a
function of the form

$\mathcal{F}_{m}={A}\,p^{m}+B,$ (1)

so-called an average sequence fidelity, where here $p\lesssim 1$ is a _quality
factor_ capturing the noise solely associated to gates, and both
$0\leq{A,B}\leq 1$ are constants isolating the SPAM errors. The average gate-
fidelity 111 We point out that the average gate-fidelity of a quantum channel
$\tilde{\mathcal{G}}$ with respect to a gate $\mathcal{G}$, is a measure of
their “average orthogonality”, rather than their distinguishability (which
albeit related, are not quite the same thing):
$\mathrm{F}_{\mathsf{avg}}:=\int{d}\psi\operatorname{tr}[\tilde{\mathcal{G}}(\psi)\mathcal{G}(\psi)]$,
where $\psi$ are (uniformly distributed) pure states. of the noisy gates with
respect to the ideal ones, $\mathrm{F}_{\mathsf{avg}}$, is then related to the
quality factor $p$ simply as (see e.g., [4])

$\mathrm{F}_{\mathsf{avg}}=p+\frac{(1-p)}{d_{\mathsf{S}}}.$ (2)

The Markovian assumption is effectively equivalent to an environment,
$\mathsf{E}$, dissipating or _forgetting_ any information of its interaction
with system $\mathsf{S}$ at any given time, thereby just introducing noise
locally at such time and reducing the purity of the noisy outputs, effectively
as if $\mathsf{E}$ was not there at all. But several factors in the
advancement and scaling up of quantum devices, including being able to probe
smaller timescales and having larger systems strongly interacting with each
other, make the Markov approximation implausible. Non-Markovianity, i.e., the
presence of (non-negligible) temporal correlations, implies that anything
other than $\mathsf{S}$, in particular other qubits, can play the role of an
$\mathsf{E}$ _with memory_ propagating undesired noise correlations in time
[44], and that noisy processes cannot be thought of as individual quantum
channels associated independently to the ideal quantum gates.

A formal definition of non-Markovianity can be seen in Appendix A.2; in
particular, we work within the process tensor framework [19], which
generalizes the notion of stochastic processes to quantum theories [45] and
contains both the classical notion and the various criteria for quantum non-
Markovianity known hitherto [20, 21].

It has been shown in [25] that the generalization of Eq. (1) to the non-
Markovian gate-independent 222 Non-Markovianity _can_ be said to introduce a
type of gate-dependence; nevertheless, here by explicit gate-independence we
mean noise associated to $\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G}$, where
$\mathcal{I}_{\mathsf{E}}$ is an identity map on $\mathsf{E}$ and
$\mathcal{G}$ is an ideal gate, is not explicitly dependent on $\mathcal{G}$.
noise case, has the (generally non-exponential) form

$\mathcal{F}_{m}=\operatorname{tr}\left\\{M\operatorname{tr}_{\mathsf{E}}\left[\left(\mathcal{Q}^{\prime}_{m,A}+\mathcal{Q}^{\prime}_{m,B}\right)\rho\right]\right\\},$
(3)

where the $\mathcal{Q}^{\prime}_{m,\bullet}$ are generalizations of quality
factors, so-called _quality maps_ of associated sequence length $m$, and here
$\rho$ is a noisy $\mathsf{S}\mathsf{E}$ initial state depending on the
prepared $\rho_{\mathsf{S}}$ (which can generally get correlated with
$\mathsf{E}$). Quality maps are, generically, Completely Positive (CP) maps
which can be further understood as averaged multi-linear maps (so-called
process-tensors [21]) taking the set of $m$ ideal digital gates
$\\{\mathcal{G}_{i}\\}_{i=1}^{m}$ as input and being uniformly averaged over
all of these. More generally than Markovian quality factors, they encode both
SPAM and fidelity-like information of the noise across a whole RB process,
reducing to $A$, $B$ and $p$ of Eq. (1) in the Markovian limit. In general,
Eq. (3) is difficult to work with even in the time-stationary 333 Here with
noise being time-stationary, albeit non-Markovian, we mean that noise
associated to $\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G}_{i}$, where
$\mathcal{I}_{\mathsf{E}}$ is an identity channel on $\mathsf{E}$ and
$\mathcal{G}_{i}$ is the ideal gate in $\mathsf{S}$ at timestep $i$, is
independent of the timestep $i$. noise regime, as opposed to the Markovian
case, although it has been previously studied in [24, 25].

To derive both the decays in Eq. (1) and Eq. (3), and an explicit mathematical
expression of all the quantities involved, we need a mathematical model of
both Markovian and non-Markovian noise.

###  II.3 Modeling of quantum noise

Markovian noise can be generically modeled by defining noisy gates, now
possibly non-unitary quantum channels, as
$\tilde{\mathcal{G}}:=\Lambda\circ\mathcal{G}$, where $\mathcal{G}$ is the
ideal digital gate, and $\Lambda$ is any CP (trace non-increasing) map on
$\mathsf{S}$ which we refer to as a noise channel. In general, noise channels
could depend on which specific gate is applied,
$\tilde{\mathcal{G}}_{i}=\Lambda_{\mathcal{G}}\circ\mathcal{G}_{i}$, or at
which time-step $i$ in a gate-sequence such gate is applied,
$\tilde{\mathcal{G}}_{i}=\Lambda_{i}\circ\mathcal{G}$ (or, of course, depend
on both, $\tilde{\mathcal{G}}_{i}=\Lambda_{\mathcal{G}_{i}}\circ\mathcal{G}$).
To generalize this to non-Markovian (temporally correlated) noise we have to
take $\mathsf{E}$ into account. That is, we now model the whole
$\mathsf{S}\mathsf{E}$ noisy map
$\tilde{\mathcal{G}}:=\Lambda\circ(\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G})$,
with $\Lambda$ now a CP map on $\mathsf{S}\mathsf{E}$. Sequential applications
$\tilde{\mathcal{G}}_{j}\circ\cdots\circ\tilde{\mathcal{G}}_{k}$ can give rise
to temporal correlations among the corresponding $\Lambda$ noise maps, as
formalized by the process tensor framework [19, 21].

Figure 2: Modeling RB with non-Markovian noise: Sample circuit of a RB
experiment on a quantum system $\mathsf{S}$ with input state
$\rho_{\mathsf{S}}$, gate sequence $\mathcal{G}_{1}$,…,$\mathcal{G}_{m}$,
undo-gate
$\mathcal{G}_{m+1}:=\mathcal{G}_{1}^{-1}\circ\cdots\circ\mathcal{G}_{m}^{-1}$,
and measurement operator $M$; an underlying non-Markovian noise process can be
represented by a fiducial initial state $\rho_{\mathsf{E}}$ on an
environmental quantum system $\mathsf{E}$ together with a set of CP, trace
non-increasing maps $\Lambda_{0}$, …, $\Lambda_{m+1}$ acting jointly on
$\mathsf{S}\mathsf{E}$.

The average sequence fidelity in Eq. (3) corresponds simply to the explicit
evaluation of the uniform average over Clifford gates,
$\mathbf{E}_{\mathcal{G}}$, in

$\mathcal{F}_{m}:=\mathbf{E}_{\mathcal{G}}\\!\left\\{\operatorname{tr}\\!\left[M\operatorname{tr}_{\mathsf{E}}\circ\Lambda\\!\circ\\!\mathcal{G}_{m+1}\\!\circ\\!\cdots\\!\circ\\!\Lambda\\!\circ\\!\mathcal{G}_{1}(\rho)\right]\right\\},$
(4)

where
$\mathcal{G}_{m+1}:=\mathcal{G}_{1}^{-1}\circ\cdots\circ\mathcal{G}_{m}^{-1}$
and $\rho:=\Lambda(\rho_{\mathsf{E}}\otimes\rho_{\mathsf{S}})$ for given
fiducial initial states of $\mathsf{S}$ and $\mathsf{E}$, $\rho_{\mathsf{S}}$
and $\rho_{\mathsf{E}}$, respectively. This is such that when
$\Lambda=\mathcal{I}$, the average sequence fidelity equals
$\operatorname{tr}\\!\left(M\rho_{\mathsf{S}}\right)$. Implicitly, together
the gate-independent and time-stationary noise assumptions mean
$\Lambda_{\mathcal{G}_{i}}=\Lambda_{i}=\Lambda$, for all gates
$\mathcal{G}_{i}$ and time-steps $i$. A circuit representation of a RB
sequence sample with non-Markovian time-non-stationary noise can be seen in
Fig. 2; terms $\Lambda_{0}$ and $\Lambda_{m+1}$ are interpreted jointly as
SPAM noise terms.

As detailed in Appendix A.4, in the Markovian case, Eq. (4) leads to Eq. (1)
with $p$ depending solely on the average gate-fidelity of $\Lambda$, as in Eq
(20), and $A$, $B$ depending solely on SPAM errors, as in Eq. (21). More
generally, in [24], Eq. (4) was shown to be of the form of Eq. (3), with the
quality factors $\mathcal{Q}_{m,A}^{\prime}$ and $\mathcal{Q}_{m,B}^{\prime}$
given explicitly as in Eq. (26) of Appendix B. This implies as well that, in
the Markovian, time-stationary noise approximation regime, Eq. (3) reduces to
the exponential decay of Eq. (1).

We wish to be able to exploit the simplicity of the RB protocol in spite of
noise being non-Markovian. In the following, we show that precisely this
reduction of the non-Markovian case to the Markovian one can be effectively
achieved operationally by interleaving DD sequences within the RB protocol.

##  III Markovianization of RB with DD

The way DD works is by applying a sequence of pulses via a control
Hamiltonian, $H_{\mathsf{ctrl}}$, in a way that effectively averages out
undesired coupling terms in the total free-evolution Hamiltonian, $H$, which
for finite-dimensional systems other than having a finite largest singular
value, can be arbitrary [48]. For infinite-dimensional $\mathsf{E}$, at least
a frequency cut-off is required [48], but a broader classification for which
Hamiltonians are amenable to DD exists [49, 50]. Generally, DD can effectively
decouple, albeit partially, a wide class of $\mathsf{S}\mathsf{E}$
Hamiltonians. Here we will consider Universal Dynamical Decoupling (UDD),
which is universal in the sense that it averages out errors up to the first
order in the Magnus expansion of the evolution [34], independently of $H$.
Concretely, we refer to a unitary group $\mathds{V}$ on $\mathsf{S}$ as a
universally decoupling group whenever
$\sum_{v\in\mathds{V}}vXv^{\dagger}=O_{\mathsf{E}}\otimes\mathds{1}_{\mathsf{S}}$,
for any $\mathsf{S}\mathsf{E}$ operator $X$ and some $\mathsf{E}$ operator
$O_{\mathsf{E}}$ and $\mathds{1}_{\mathsf{S}}$ the identity on $\mathsf{S}$.

In the setting of RB, we can model free evolution at any given time-step as
the underlying noise $\Lambda$ on the whole $\mathsf{S}\mathsf{E}$, so
_Markovianizing_ Eq. (3) can effectively be accomplished by applying UDD
between RB gates. As detailed in Appendix C, we consider ideal pulses
generated by
$H_{\mathsf{ctrl}}(t)=\frac{\pi}{2}\sum_{k}\delta(t\\!-\\!t_{k})v_{k}$ at
times $t_{1}<t_{2}<\ldots<t_{\eta}$, where $\delta$ is a Dirac delta and
$\eta=|\mathds{V}|$ is the number of elements in the decoupling group, i.e.,
infinitely strong instantaneous pulses with decoupling operators $v_{k}$ at
times $t_{k}$.

Figure 3: Interleaved UDD on RB: Circuit depictions of (a) standard RB with
initial state $\rho_{\mathsf{S}}$, measurement $M$, and time-intervals
$\tau_{\mathsf{rb}}$ with noise $\Lambda^{(\tau_{\mathsf{rb}})}$ between
application of gates $\mathcal{G}_{i}$, (b) the same RB-sequence, now with a
single interleaved UDD sequence
$\mathcal{S}_{\mathsf{dd}}^{(\tau_{\mathsf{dd}})}$, leaving environment
correlations only at $\mathcal{O}(\tau_{\mathsf{dd}}^{2})$, and (c) the
sequence $\mathcal{S}_{\mathsf{dd}}^{(\tau_{\mathsf{dd}})}$ with non-identity
pulses $\tilde{\mathcal{V}}_{i}\in\mathds{V}$ applied in time-intervals of
$\tau_{\mathsf{dd}}$, and total elapsed time $\eta\tau_{\mathsf{dd}}$, where
$\eta=|\mathds{V}|$.

Let us label the noise maps by an elapsed-time $t$ between subsequent
applications of any two gates as $\Lambda^{(t)}$. Denoting the channels
associated to the ideal decoupling operators as
$\mathcal{V}(\cdot):=v(\cdot)v^{\dagger}$, these pulses can be applied at
evenly-spaced time-intervals $\tau_{\mathsf{dd}}$ through
$\mathcal{S}_{\mathsf{dd}}^{(\tau_{\mathsf{dd}})}:=\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{\mathcal{V}\in\mathds{V}}\left(\mathcal{V}\circ\Lambda^{(\tau_{\mathsf{dd}})}\circ\mathcal{V}^{\dagger}\right)$.
Notice that we are taking the application of DD pulses to be strictly evenly
spaced. As we will interleave these among the random RB gates, we will
necessarily have free-evolution (noise) on the edge terms between the
application of gates; taking either the first or last pulse as the identity,
we can equivalently write

$\mathcal{S}^{(\tau_{\mathsf{dd}})}_{\mathsf{dd}}=\Lambda^{(\tau_{\mathsf{dd}}/2)}\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{\tilde{\mathcal{V}}}\left(\tilde{\mathcal{V}}\circ\Lambda^{(\tau_{\mathsf{dd}})}\circ\tilde{\mathcal{V}}^{\dagger}\right)\circ\Lambda^{(\tau_{\mathsf{dd}}/2)},$
(5)

with $\tilde{\mathcal{V}}$ being non-identity elements of $\mathds{V}$.

For example, for a single-qubit, $\mathds{V}=\\{\mathds{1},X,Y,Z\\}$, the
Pauli group, so that
$\mathcal{S}^{(\tau_{\mathsf{dd}})}_{\mathsf{dd}}=\mathcal{Z}\circ\Lambda^{(\tau_{\mathsf{dd}})}\circ\mathcal{Z}\circ\mathcal{Y}\circ\Lambda^{(\tau_{\mathsf{dd}})}\circ\mathcal{Y}\circ\mathcal{X}\circ\Lambda^{(\tau_{\mathsf{dd}})}\circ\mathcal{X}\circ\Lambda^{(\tau_{\mathsf{dd}})}$,
where curly letters here are the maps associated to each Pauli operator, and
which can be seen to be equivalent to the so-called $XY4$ and $XZ4$ sequences.

If we label by $\tau_{\mathsf{rb}}$ the original elapsed time between
application of RB Clifford gates, and if $\tau_{\mathsf{rb}}$ further turns
out to be the _minimum_ time between application of two subsequent gates on
the given device, then we would generally need to consider no less than
$\tau_{\mathsf{dd}}=2\tau_{\mathsf{rb}}$, due to the edge terms in Eq. (5).
Other considerations could come into play here in order to optimize
$\tau_{\mathsf{dd}}$, e.g., the fact that Clifford gates are composite or in
general whether DD pulses, in particular, can be applied with
$\tau_{\mathsf{dd}}<\tau_{\mathsf{rb}}$.

We can model the application of ideal UDD within a RB protocol by interleaving
a single sequence $\mathcal{S}^{(\tau_{\mathsf{dd}})}_{\mathsf{dd}}$ among
each ideal gate of a RB circuit, i.e.,
$\mathcal{S}^{(\tau_{\mathsf{dd}})}_{\mathsf{dd}}\circ\mathcal{G}_{m+1}\circ\cdots\circ\mathcal{S}^{(\tau_{\mathsf{dd}})}_{\mathsf{dd}}\circ\mathcal{G}_{1}\circ\mathcal{S}^{(\tau_{\mathsf{dd}})}_{\mathsf{dd}}$,
as depicted in Fig. 3. Notice this will change the circuit depth, and in
particular it will modify the total time between application of the RB gates
to $\eta\,\tau_{\mathsf{dd}}$; nevertheless, we obtain the following:

Result 1.  Let $\Lambda^{(t)}=\mathrm{e}^{t\mathcal{L}}$, where
$\mathcal{L}(\cdot):=-i[H,\cdot]+\mathcal{D}(\cdot)$ for a
$\mathsf{S}\mathsf{E}$ time-independent Hamiltonian $H$ and a dissipator
$\mathcal{D}=\mathcal{D}_{\mathsf{S}}+\mathcal{D}_{\mathsf{E}}$ with only
local $\mathsf{S}$ and/or $\mathsf{E}$ contributions, and let $\gamma$ be the
diagonal matrix of relaxation rates of
$\mathcal{D}_{\mathsf{S}}:=\sum_{k}\gamma_{k}[L_{k}(\cdot)L_{k}^{\dagger}-\frac{1}{2}\\{L_{k}^{\dagger}{L}_{k},\cdot\\}]$,
for some (unit-less) traceless and orthonormal operators $L_{k}$. Then for
$\tau_{\mathsf{dd}}\ll 1/\operatorname{tr}(\gamma)$, the average sequence
fidelity of length $m$ for a RB experiment under time-stationary noise with
single interleaved UDD sequences of the form of Eq. (5) satisfies,

$\mathcal{F}_{m}=A\,p_{\tau_{\mathsf{dd}}}^{m}+B+\mathcal{O}(\tau_{\mathsf{dd}}^{2}),$
(6)

where $A,B$ are SPAM constants for fixed $\tau_{\mathsf{dd}}$, and the quality
factor is a $\mathcal{O}(\tau_{\mathsf{dd}})$ term

$p_{\tau_{\mathsf{dd}}}=1-\eta\,\tau_{\mathsf{dd}}\,\frac{\operatorname{tr}(\gamma)}{d_{\mathsf{S}}-\frac{1}{d_{\mathsf{S}}}},$
(7)

where $\eta=|\mathds{V}|$ is the number of elements of the decoupling group
$\mathds{V}$.

The proof is shown in Appendix C and the value of all constants is shown
explicitly in Eq. (50-52).

The main message in Result 3 is that ideal-pulse UDD Markovianizes RB for a
broad class of non-Markovian noise models, in the sense of reducing the
average sequence fidelity of Eq. (3) to that of Eq. (1), to the first order in
the DD sequence interval time $\tau_{\mathsf{dd}}$. Notice the limiting case
of $\tau_{\mathsf{dd}}\to 0$ would be that when both the ideal DD pulses _and_
the random gates in the RB sequence are implemented infinitely fast, so that,
of course, $\mathcal{F}_{m}=1$, as no added idling time is being considered.

The continuous dynamical model for the noise as a $\mathsf{S}\mathsf{E}$
Lindblad evolution encompasses a broad class of noise models, where the
$\mathsf{S}\mathsf{E}$ dynamics are Markovian while the reduced ones on
$\mathsf{S}$ are non-Markovian. The only restriction is that the dissipator
term, $\mathcal{D}$, cannot have global $\mathsf{S}\mathsf{E}$ contributions
to the first $\tau_{\mathsf{dd}}$ order for our result to hold, which would be
the case if correlated $\mathsf{S}\mathsf{E}$ information is dissipated. While
in Section V.1 we take as an example a two-qubit system, motivated by more
realistic two-level defect noise models, relevant in superconducting systems,
many other realistic models of non-Markovianity would consider an infinite-
dimensional $\mathsf{E}$. We must point out that the only limiting factor to
our theory in such cases would be exotic models having parameter domains where
the dynamics cannot be dynamically decoupled [51].

In the absence of dissipation terms,
$\mathcal{F}_{m}=1+\mathcal{O}(\tau_{\mathsf{dd}}^{2})$, i.e., the average
sequence fidelity is identically one, albeit only to the first
$\tau_{\mathsf{dd}}$ order. Notice that the modeling of noise in Result 3
explicitly takes into account the presence or absence of dissipation terms
between applications of RB gates and DD pulses: e.g., Result 3 will hold with
a dilated Hamiltonian $H^{\prime}$ on a larger system
$\mathsf{S}\mathsf{E}\mathsf{E}^{\prime}$, only if the absence of global
dissipation on $\mathsf{S}\mathsf{E}$ is also implied.

The quality factor $p_{\tau_{\mathsf{dd}}}$ is derived in Appendix C and it is
related to the average gate-fidelity of the channel
$\mathcal{I}+\tau_{\mathsf{dd}}\,\mathcal{D}_{\mathsf{S}}$ with respect to the
identity $\mathcal{I}$; i.e., it quantifies the noise contribution from the
dissipator $\mathcal{D}_{\mathsf{S}}$, which is the sole generator of
Markovian noise once the Hamiltonian has been averaged out in $\mathsf{S}$ and
it is a trace-zero map, not CP in general. The final expression in terms of
the relaxation rates, $\gamma_{k}$, can then be seen to follow.

Higher-order, $\mathcal{O}\left(\tau_{\mathsf{dd}}^{2}\right)$ terms, albeit
still also containing Markovian contributions, will, in general, be non-
Markovian; so while non-exponential deviations may get suppressed and allow to
fit an exponential, the resulting decay is purely Markovian only to first
order. Thus we can interpret the efficacy of digital UDD sequences, in both
Markovianization of the average sequence fidelity and improvement in overall
average gate-fidelity, mainly in terms of how fast they can be applied,
relative to the time between application of the RB gates,
$\tau_{\mathsf{rb}}$. Given that multi-qubit Cliffords are composite gates, it
might be possible to implement DD pulses in time-scales
$\tau_{\mathsf{dd}}<\tau_{\mathsf{rb}}$. As we see in Section V.1, in practice
this might be the main limiting factor, together with considerations such as
non-ideal, finite-width pulses, as well as imperfect applications of these.

Finally, realistic DD cannot be implemented as ideal instantaneous pulses, and
these themselves can introduce control errors [52]. Furthermore, it is known
that finite-width DD does not achieve perfect decoupling even to first time-
order [34]. However, as long as pulses are sufficiently narrow, the noise they
introduce is sufficiently small [52] and local, UDD will Markovianize RB to an
extent close to the one predicted by Result 3. One can furthermore employ more
elaborated DD techniques such as Concatenated DD [48], or devised optimized
decoupling sequences, although we do not pursue that here.

##  IV RB under Pauli-twirled noise

While the Markovianizing effect of DD in RB is somehow expected, a prominent
error-suppression technique that has recently been shown to suppress non-
Markovian noise in a statistically significant way [35] is Randomized
Compiling (RC) [53, 40]. RC can be understood as the operational way of
tailoring arbitrary Markovian noise quantum channels into Pauli channels,
which mathematically corresponds to a mapping known as Pauli-twirling. As
opposed to DD, RC is not a quantum control technique, but instead it relies on
compiling a set of logically-equivalent circuits (i.e., with no increase in
depth and effectively implementing the same quantum operations) where noisy
gates are dressed with uniformly sampled random single-qubit Pauli gates;
averaging over all such circuits, approximately and efficiently implements a
Pauli twirl.

Figure 4: RB under non-Markovian $\mathsf{S}$-Pauli twirled noise: In (a), a
sample circuit for a standard RB experiment where noise
$\Lambda_{i}^{\mathsf{P}_{\mathsf{S}}}$ is time-non-stationary, non-Markovian,
and has been Pauli-twirled on subsystem $\mathsf{S}$, while in (b) the
operational definition for $\Lambda_{i}^{\mathsf{P}_{\mathsf{S}}}$ is shown
explicitly as an average over ideal $n_{s}$-qubit Pauli terms $P$; this is
known as a $\mathcal{G}_{i}$-twisted twirl [12] on $\mathsf{S}$, and RC
accurately approximates it efficiently by randomly sampling single-qubit Pauli
terms and compiling the outputs [53]. Only here we employ the notation
$\mathcal{G}^{\dagger}(\cdot):=G^{\dagger}(\cdot)G$.

Given any $n_{s}$-qubit quantum channel $\Phi$, RC effectively and efficiently
maps its $\chi$-matrix representation,
$\Phi(\cdot):=\sum_{ik}\chi_{ik}P_{i}\,(\cdot)P_{k}$, to
$\Phi^{\mathsf{P}}(\cdot)=\sum_{i}\chi_{ii}P_{i}\,(\cdot)P_{i}$, where here
$P$ are $n_{s}$-qubit Pauli operators. The channel $\Phi^{\mathsf{P}}$ is
generally a non-unitary (stochastic or incoherent) quantum channel known as a
Pauli channel, and the $\alpha_{i}:=\chi_{ii}$ define a probability
distribution called Pauli error rates. It can be seen that average gate-
fidelity (and thus RB) estimates precisely a $\chi_{00}=\alpha_{0}$ term,
corresponding to the identity term, $P_{0}:=\mathds{1}$ (probability of no
error happening). Crucially, RC suppresses all off-diagonal terms of the
$\chi$ matrix, including terms associated with coherent errors, which can give
worst-case error rates orders of magnitude higher than stochastic errors [54],
while leaving average gate-fidelity unchanged.

For the non-Markovian case, a Pauli-twirl solely on system $\mathsf{S}$ on any
given $\mathsf{S}\mathsf{E}$ quantum channel $\Lambda$, in the $\chi$-matrix
representation, takes the form

$\Lambda^{\mathsf{P}_{\mathsf{S}}}(\cdot):=\sum\chi_{\mu\nu,ii}(R_{\mu}\otimes{P}_{i})(\cdot)(R_{\nu}^{\dagger}\otimes{P}_{i}),$
(8)

for some basis operators $R$ on $\mathsf{E}$, and coefficients
$\chi_{\mu\nu,ii}$ of the corresponding $\mathsf{S}\mathsf{E}$ Hermitian
$\chi$-matrix. If we were to trace out $\mathsf{E}$ on this channel, then, of
course, the reduced channel is a Pauli channel,
$\operatorname{tr}_{\mathsf{E}}\Lambda^{{\mathsf{P}_{\mathsf{S}}}}\left(\rho_{\mathsf{E}}\otimes\cdot\right)=\sum\alpha_{i}(\rho_{\mathsf{E}}^{\prime})P_{i}(\cdot)P_{i}$,
with $\alpha_{i}$ being Pauli error-rates dependent on the
$\mathsf{E}$-reduced output state $\rho_{\mathsf{E}}^{\prime}$. This tracing,
however, generally occurs at the end of a computation or any multi-time
process.

Operationally, however, the mapping
$\Lambda\mapsto\Lambda^{\mathsf{P}_{\mathsf{S}}}$ can only be implemented by
modifying the noisy gates
$\tilde{\mathcal{G}}:=\Lambda\circ(\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G})$,
as we don’t have direct access to $\Lambda$. This can be done through a so-
called $\mathcal{G}$-twisted twirl [12] on system $\mathsf{S}$, defined as

$\tilde{\mathcal{G}}(\cdot)\mapsto
4^{-n_{s}}\sum_{P}P\tilde{\mathcal{G}}\left(\mathcal{G}^{\dagger}(P)\,(\cdot)\,\mathcal{G}(P)\right)P,$
(9)

with sum over $P$ all $n_{s}$-qubit ideal Pauli operators, and where (abusing
notation) we denote as $\mathcal{G}^{\dagger}(\cdot):=G^{\dagger}(\cdot)G$ the
adjoint map of the noiseless Clifford gate $\mathcal{G}$. This is depicted
schematically in Fig. 4. RC thus approximates $\mathsf{S}$-Pauli twirls by
randomly sampling $\mathcal{G}_{i}$-twisted twirls on all time steps and
recompiling the Pauli gates, which can be done with a number of samples much
smaller than $4^{n_{s}}$ for larger $n_{s}$ [40].

Similar to DD, where we worked with the ideal pulse limit, here we focus on
the case where RC fully and perfectly tailors $\Lambda$ into
$\Lambda^{\mathsf{P}}$. While for a single-qubit this can be done exactly, in
general, the tailoring of noise into Pauli noise by RC can be quantified to
occur to a large percentage with a small number of samples [53]. We thus point
out the following:

Result 2.  For any RB sequence fidelity given as
$\mathrm{f}_{m}[\\{\Lambda_{i}\\}_{i=0}^{m+1}]:=\operatorname{tr}[{M}\operatorname{tr}_{\mathsf{E}}\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{i=1}^{m+1}(\Lambda_{i}\circ\mathcal{G}_{i})\,\rho]$
of length $m$, where here
$\rho:=\Lambda_{0}(\rho_{\mathsf{E}}\otimes\rho_{\mathsf{S}})$, the
corresponding average sequence fidelity with $\mathsf{S}$-Pauli twirled noise,
$\mathbf{E}_{\mathcal{G}}\mathrm{f}_{m}[\\{\Lambda_{i}^{\mathsf{P}_{\mathsf{S}}}\\}_{i=0}^{m+1}]$,
remains in general non-exponential.

Furthermore, when averaging is over uniformly distributed $n_{s}$-qubit
Clifford gates,

$\mathbf{E}_{\mathcal{G}}\mathrm{f}_{m}\left[\\{\Lambda_{i}\\}_{i=0}^{m+1}\right]\\!=\\!\mathbf{E}_{\mathcal{G}}\mathrm{f}_{m}\left[\\{\Lambda_{0},\Lambda_{m+1}\\}\\!\cup\\!\\{\Lambda_{i}^{\mathsf{P}_{\mathsf{S}}}\\}_{i=1}^{m}\right],$
(10)

where here $\cup$ denotes the union of sets. On the other hand,

$\mathbf{V}_{\mathcal{G}}\mathrm{f}_{m}\left[\\{\Lambda_{i}^{\mathsf{P}_{\mathsf{S}}}\\}_{i=0}^{m+1}\right]\leq\mathbf{V}_{\mathcal{G}}\mathrm{f}_{m}\left[\\{\Lambda_{i}\\}_{i=0}^{m+1}\right],$
(11)

where $\mathbf{V}_{\mathcal{G}}$ is variance over arbitrary gates
$\mathcal{G}$.

The proof can be seen in Appendix D.

The first statement means that Pauli twirling does not Markovianize the
average sequence fidelity as DD does, i.e., it does not turn non-Markovian
non-exponential RB decays into exponential ones. In particular, this statement
holds for any gate set other than the multi-qubit Clifford group, at least as
long as it forms a finite group.

On the other hand, the second statement in Eq. (10) holds only for the
Clifford group but it essentially means that RC, or any Pauli-noise tailoring
technique, would at most have an effect (not necessarily Markovianizing or
increasing the average sequence fidelity) on SPAM contributions when it comes
to non-Markovian RB data. The reason for this is not RB per-se, but rather
average gate-fidelity as a figure of merit, as it only takes into account the
probability of no error happening on $\mathsf{S}$ and not all the other error
terms that Pauli twirling eliminates. That is, average gate-fidelity by
definition is proportional only to the zeroth (identity) element of the
$\chi$-matrix; explicitly, see Eq. (64,65) in Appendix D.1. This being said it
is still possible that there is a class of noise profiles where it is
sufficient to $\mathsf{S}$-Pauli-twirl the SPAM noise in order to Markovianize
the average sequence fidelity.

Despite not Markovianizing the average sequence fidelity, Eq. (11) establishes
that RC will, in general, reduce the variance of sequence fidelities 444 The
statement in Eq. (11) holds more generally for $\mathrm{f}_{m}$ being noisy
expectation values.. The importance of this is twofold, it allows for both a
confident diagnosis of non-exponential deviations and for an accurate
estimation of meaningful error rates. While mathematically this result can be
seen to follow simply because Pauli-twirling removes additive terms to the
variance, it can be argued that, physically, the reason is that it reduces the
coherence of noise, as precisely the expectation of the squared sequence
fidelity is proportional to the so-called average unitarity. The unitarity,
$u$, is a figure of merit quantifying loss of purity due to noise [7], and
satisfies $u=1$ if noise is coherent, i.e., due to a unitary map, and $u\leq
1$ otherwise. In the Markov case, $u\geq{p^{2}}$ [56], which is saturated for
depolarizing noise, i.e., Pauli channels with all $\chi_{ii}$ terms, for
$i\neq 0$, being equal. In the non-Markov case, it is _expected_ as well that
$\mathsf{S}$-Pauli-twirling would only decrease or leave the total unitarity
unchanged 555 In the $\mathsf{S}\mathsf{E}$ case, it is unclear whether
$\mathsf{S}$-Pauli-twirling _indeed will always_ either only decrease or leave
the total unitarity unchanged, due precisely to $\mathsf{E}$; here the role of
the average sequence fidelity would also come into play to always decrease the
variance.. This, in a sense, was already pointed out in the original unitarity
benchmarking proposal of [7], although as also mentioned there, it is less
straightforward to give a concrete bound.

Finally, similar to the case of DD, here we have considered an exact
implementation of Pauli twirling, i.e., a perfect application of RC. While
this is unrealistic, in general, RC incurs in only a small Pauli sampling
overhead to closely approximate an exact Pauli twirl [40]. We now give two
concrete numerical examples of our results for both DD and RC within non-
Markovian RB.

##  V Numerical Examples

###  V.1 XY4 to Markovianize a qubit

As proof-of-principle, we consider both one $\mathsf{S}$ qubit and one
$\mathsf{E}$ qubit with time-stationary noise
$\Lambda^{(t)}=\mathrm{e}^{t\mathcal{L}}$ acting jointly on both qubits, where
$\mathcal{L}(\cdot):=-i[H,\cdot]+\mathcal{D}_{\mathsf{S}}(\cdot)$ for a
Hamiltonian $H=JXX+h_{x}(XI+IX)+h_{y}(YI+IY)$ with some constants
$J,h_{x},h_{y}$ and local $\mathsf{S}$-dissipator
$\mathcal{D}_{\mathsf{S}}(\cdot):=\sum_{k}\gamma_{k}[L_{k}(\cdot)L_{k}^{\dagger}-\frac{1}{2}\\{L_{k}^{\dagger}{L}_{k},\cdot\\}]$
with $L_{0}=X$ and $L_{1}=Z$. While this is an arbitrary dynamical model,
similar qubit-to-qubit noise mechanisms are conceivable to arise e.g., with
two-level system defects in superconducting qubits [58], albeit with distinct
spin interactions, dissipators, and varying strengths.

Figure 5: Interleaved DD in non-Markovian RB: Average sequence fidelities
($\mathcal{F}_{m}$) for RB experiments on increasing sequence lengths ($m$)
with the full, generally non-Markovian noise model
$\Lambda^{(\tau_{\mathsf{rb}})}$ (blue triangles), and interleaved ideal $XY4$
sequences (green closed markers and black open markers) with time-intervals
between pulses of $\tau_{\mathsf{dd}}=\tau/2$, $\tau$ and $2\tau$, for a fixed
$\tau=\tau_{\mathsf{rb}}=0.03$. All numerical averages were taken over 40
samples, with red bands denoting uncertainty of mean; non-Markovian analytical
decay (blue upper triangle) computed according to [24] and
$\mathcal{O}(\tau_{\mathsf{dd}})$ analytical decays (black open markers) up to
first time-order according to Eq. (6). 1$\mathsf{S}$+1$\mathsf{E}$ qubit model
given by $\Lambda^{(t)}=\mathrm{e}^{t\mathcal{L}}$ where
$\mathcal{L}(\cdot):=-i[H,\cdot]+\mathcal{D}_{\mathsf{S}}(\cdot)$ for
$H=JXX+h_{x}(XI+IX)+h_{y}(YI+IY)$ and local $\mathsf{S}$ dissipator
$\mathcal{D}_{\mathsf{S}}(\cdot):=\sum_{k}\gamma_{k}[L_{k}(\cdot)L_{k}^{\dagger}-\frac{1}{2}\\{L_{k}^{\dagger}{L}_{k},\cdot\\}]$;
chosen parameters $\rho_{\mathsf{S}}=M=|0\rangle\\!\langle 0|$, $J=1.7$,
$h_{x}=1.47$, $h_{y}=-1.05$ and Lindblad terms $L_{0}=X$, $\gamma_{0}=0.002$,
$L_{1}=Z$, $\gamma_{1}=0.007$; Lindblad evolution approximated up to
$\mathcal{O}(t^{10})$.

We compare interleaved DD pulses given by the $XY4$ sequence, briefly
described below Eq. (5) and depicted in Fig. 3, for DD sequences of
$\tau_{\mathsf{dd}}=\tau/2$, $\tau$, and $2\tau$ time-intervals, where we
fixed $\tau:=\tau_{\mathsf{rb}}$ to be the evolution time of the noise between
applications of random Clifford gates in the original non-Markovian RB
experiment. Outputs for particular constants in the noise model
$\Lambda^{(t)}$ can be seen in Fig. 5, where the non-Markovian analytical
average sequence fidelity was computed according to [24], and the analytical
decays to first $\tau_{\mathsf{dd}}$ orders were computed according to Eq. (6)
where $\eta=4$.

The main message is that DD effectively removes non-Markovian non-exponential
deviations for a time scale at which DD sequences are applied
$\tau_{\mathsf{dd}}\ll 1/\operatorname{tr}(\gamma)$, and how well it does so
–whether it outputs decays close to the purely exponential
$\mathcal{O}(\tau_{\mathsf{dd}})$ decay–, depends mainly on whether the
$\tau_{\mathsf{dd}}$ is small relative to that of the noise,
$\tau_{\mathsf{rb}}$, in the original non-Markovian RB experiment.

In particular, interleaving DD pulses in time-intervals of
$\tau_{\mathsf{dd}}=2\tau_{\mathsf{rb}}$ does reduce non-exponential
deviations, but still, non-Markovian noise dominates the decay. In general,
this can be assessed by inspecting the separation of the numerical,
all-$\tau_{\mathsf{dd}}$-order average decays with respect to the
corresponding only first-$\tau_{\mathsf{dd}}$-order analytical decays; further
discussion and an analysis in the time-scales for $\tau_{\mathsf{dd}}$
relative to the chosen $\tau_{\mathsf{rb}}$ can be seen in Appendix C.1.4, as
captured by Fig. 9, where it is clear that this separation is suppressed and
decays asymptotically become purely Markovian for
$\tau_{\mathsf{dd}}<\tau_{\mathsf{rb}}$.

Taking into account that pulses will have a finite width and might themselves
be noisy due to implementation errors, the ability to perform good pulses in
time scales shorter than that of the RB gates, would lead to a useful
reduction of non-exponential deviations and enhancement of average sequence
fidelities, as shown by the plots for time-intervals
$\tau_{\mathsf{dd}}=\tau_{\mathsf{rb}}$ and
$\tau_{\mathsf{dd}}=\tau_{\mathsf{rb}}/2$. Comparison with the corresponding
analytical $\tau_{\mathsf{dd}}$ decays, up to the first time-order, can help
to estimate the contribution of non-Markovian and higher-order Markovian terms
in the model.

###  V.2 Subsystem Pauli-twirled noise

We now consider the same 1$\mathsf{S}$+1$\mathsf{E}$ qubit model but where
instead of interleaving $XY4$ pulses, we consider noise that has been Pauli-
twirled on qubit $\mathsf{S}$ (e.g., via RC). We now fix a value for $\tau$
and simply denote $\Lambda$ as the noise channel for all time steps. We
distinguish between SPAM noise terms $\Lambda_{0}=\Lambda_{m+1}$ and the
_bulk_ noise terms $\Lambda_{1}=\Lambda_{2}=\ldots=\Lambda_{m}$ (i.e., all
non-SPAM noise). Since we deal with a single-qubit, and RC does not change the
logical structure of quantum circuits (e.g., does not increase their depth),
we take the limit of perfect Pauli-twirled noise on $\mathsf{S}$.

Figure 6: The effect of RC on non-Markovian RB: For the same model and
parameters (fixed $\tau_{\mathsf{rb}}=0.03$ and
$\Lambda^{(\tau_{\mathsf{rb}})}$ truncated at
$\mathcal{O}(\tau_{\mathsf{rb}}^{10})$) as Fig. 5, in (a), a comparison of the
analytical average sequence fidelity decay (computed according to [24]) for
combinations of bare noise ($\Lambda$) and $\mathsf{S}$-Pauli twirled noise
($\Lambda^{\mathsf{P}_{\mathsf{S}}}$) on bulk terms
$\Lambda_{1}=\Lambda_{2}=\ldots=\Lambda_{m}$ and SPAM terms
$\Lambda_{0}=\Lambda_{m+1}$, showing that decays coincide up to SPAM terms;
while in (b) the numerical variance over 40 samples of RB circuits, either all
under bare noise or $\mathsf{S}$-Pauli twirled noise, showing that the
variance generally gets suppressed for the latter.

In Fig. 6, we numerically demonstrate Result 4 by comparing the analytical
average sequence fidelities (computed according to [24]) with different
choices of either bare noise or $\mathsf{S}$-Pauli twirled noise for the bulk
terms and the SPAM ones; we also show the behavior of the variance of the
sequence fidelity with respect to sequence length for a set of 40 numerical
samples, all either with bare noise or $\mathsf{S}$-Pauli twirled noise.

For the behavior of the average sequence fidelity with respect to SPAM, we
should point out that an overall increase in fidelity due to
$\mathsf{S}$-Pauli twirled noise is incidental in this example; other than
this, it is clear that decays remain non-exponential and that whenever SPAM
terms ($\Lambda_{0}$ and $\Lambda_{m+1}$) coincide, it is irrelevant for the
decay whether all other terms were $\mathsf{S}$-Pauli twirled or not. While we
do not argue about the feasibility of operationally twirling SPAM noise, it is
conceivable that this could be achieved at least partially on either the state
preparation term or the measurement one (which in practice are
indistinguishable).

For the case of the variance, we employed numerical averages on 40 RB samples.
The $\mathsf{S}$-Pauli twirled variance gets clearly suppressed, and is lower
than that of bare noise, for all sequence lengths. It stands out that the
variance for bare noise itself changes considerably with increasing sequence
length. While the variance for twirled noise is not entirely vanishing, it is
considerably suppressed, given its relation to the unitarity, as mentioned
above in Section IV, so deriving such guarantees mathematically would make a
strong point for the employment of RC in obtaining fidelities.

##  VI Conclusions and Discussion

We have shown that noise suppression techniques, such as Dynamical Decoupling
(DD) and Randomized Compiling (RC), can be efficient tools for dealing with
non-Markovian noise sources in Randomized Benchmarking (RB). In particular,
_i_) Universal Dynamical Decoupling (UDD) applied with fast and narrow pulses
reduces a wide class of non-Markovian non-exponential RB decays to an
exponential decay plus perturbative corrections in time, _ii_) RC does not
Markovianize RB in the same sense that DD does, and in fact leaves the
sequence fidelity decay invariant up to State Preparation and Measurement
(SPAM) noise, however, _iii_) RC is ensured to decrease, or at worst leave
unchanged, the variance of RB sequence fidelities.

Our results imply that standard noise suppression techniques can be valuable
tools in taming, benchmarking, and optimizing non-Markovian noise. In
particular, while we have dealt with the _standard_ (or original) RB protocol
as a benchmarking framework, this approach is amenable to be adapted to any
RB-based technique considering e.g., scalability or specific purpose metrics
[15].

In the case of DD, while Markovianization is always achieved at short times,
the full effective RB sequence fidelities can still be dominated by non-
Markovian terms, so the main limitation for DD is the time-scale at which
pulses can be applied, which for an enhancement of fidelities would need to be
shorter than those between application of the RB gates. Our main Result 3
regarding DD is stated for a broad class of continuous noise models with mild
restrictions on the global Hamiltonian and dissipators. More generally,
however, for any noise dynamics, the timescales at which decoupling can be
efficiently achieved are connected to such dynamics [59, 60], but not
necessarily to whether they display non-Markovianity [51]. Furthermore, while
they might be generally non-Markovian, some dynamics can hide non-Markovian
effects [61, 62], in particular for our case in RB, and not display
significant deviations. Finally, there is the fact that realistic DD pulses
have a finite width and themselves can introduce errors; for our results,
however, it is sufficient for them to be narrow and contain only local, time,
and gate-independent small errors.

For RC in our main Result 4, it is somewhat surprising that $\mathsf{S}$-Pauli
twirling does not Markovianize RB data since twirling also effectively
decouples by averaging out all non-Pauli error terms [35]. Concretely, RB with
the Clifford group, while it is clear that in the Markovian case it leaves
average gate-fidelities unchanged, as opposed to DD, it was expected for it to
have a decoupling effect that would impact _directly_ the average of the non-
Markovian RB outputs. Nevertheless, RC can be extremely valuable in
suppressing the uncertainty in the average sequence fidelity, which can ensure
an accurate diagnosis of non-exponential deviations (as opposed to statistical
fluctuations) and reliable estimation of error rates. Furthermore, as argued
before, this reduction in variance due to RC is related to the amount of
coherence of the noise, which in turn is intimately connected and relevant to
average gate-fidelity and fault tolerant-relevant metrics such as diamond norm
[54].

Noise suppression techniques, such as DD and RC, are a vital ingredient
allowing for the possibility to take quantum computing beyond a noise-
intermediate regime. Our results highlight the importance of incorporating
such basic noise suppression techniques to deal with one of the most
complicated sources of errors, namely non-Markovianity, not only for
deployment but for basic error diagnostics and benchmarking. As a perspective,
we expect these ideas to be useful, easily adapted, and enhanced, to other
scalable and holistic benchmarking and characterization techniques, e.g., [9,
10, 11, 12, 13], most of which are the intellectual progeny of RB.

###### Acknowledgements.

We thank Jay Nath for valuable discussions, and PFR thanks Robin Blume-Kohout
for helpful comments during the APS March Meeting 2023. KM acknowledges
support from the Australian Research Council Discovery Projects DP210100597
and DP220101793. PFR, MP, AA, and IdV acknowledge support from the German
Federal Ministry of Education and Research (BMBF) under Q-Exa (grant No.
13N16062), QSolid (grant No. 13N16161), and MUNIQC-SC (grant No. 13N16185).

## References

  * Eisert  _et al._ [2020] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi, Nat. Rev. Phys. 2, 382–390 (2020).
  * Kliesch and Roth [2021] M. Kliesch and I. Roth, PRX Quantum 2, 010201 (2021).
  * Emerson  _et al._ [2005] J. Emerson, R. Alicki, and K. Życzkowski, J. Optics B: Quantum Semiclass. Opt. 7, S347 (2005).
  * Magesan  _et al._ [2012] E. Magesan, J. M. Gambetta, and J. Emerson, Phys. Rev. A 85, 042311 (2012).
  * Knill  _et al._ [2008] E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland, Phys. Rev. A 77, 012307 (2008).
  * Gambetta  _et al._ [2012] J. M. Gambetta, A. D. Córcoles, S. T. Merkel, B. R. Johnson, J. A. Smolin, J. M. Chow, C. A. Ryan, C. Rigetti, S. Poletto, T. A. Ohki, M. B. Ketchen, and M. Steffen, Phys. Rev. Lett. 109, 240504 (2012).
  * Wallman  _et al._ [2015] J. Wallman, C. Granade, R. Harper, and S. T. Flammia, New J. Phys. 17, 113020 (2015).
  * Wood and Gambetta [2018] C. J. Wood and J. M. Gambetta, Phys. Rev. A 97, 032306 (2018).
  * Proctor  _et al._ [2022] T. Proctor, S. Seritan, K. Rudinger, E. Nielsen, R. Blume-Kohout, and K. Young, Phys. Rev. Lett. 129, 150502 (2022).
  * Erhard  _et al._ [2019] A. Erhard, J. J. Wallman, L. Postler, M. Meth, R. Stricker, E. A. Martinez, P. Schindler, T. Monz, J. Emerson, and R. Blatt, Nat. Commun. 10 (2019).
  * Helsen  _et al._ [2021] J. Helsen, M. Ioannou, I. Roth, J. Kitzinger, E. Onorati, A. H. Werner, and J. Eisert, Estimating gate-set properties from random sequences (2021), arXiv:2110.13178 [quant-ph] .
  * Flammia [2021] S. T. Flammia, Averaged circuit eigenvalue sampling (2021), arXiv:2108.05803 [quant-ph] .
  * Harper and Flammia [2023] R. Harper and S. T. Flammia, Learning correlated noise in a 39-qubit quantum processor (2023), arXiv:2303.00780 [quant-ph] .
  * Chen  _et al._ [2022] J. Chen, D. Ding, and C. Huang, PRX Quantum 3, 030320 (2022).
  * Helsen  _et al._ [2022] J. Helsen, I. Roth, E. Onorati, A. Werner, and J. Eisert, PRX Quantum 3, 020357 (2022).
  * Rivas  _et al._ [2014] A. Rivas, S. F. Huelga, and M. B. Plenio, Rep. Prog. Phys. 77, 094001 (2014).
  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, Rev. Mod. Phys. 89, 015001 (2017).
  * Pollock  _et al._ [2018a] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Phys. Rev. A 97, 012127 (2018a).
  * Pollock  _et al._ [2018b] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Phys. Rev. Lett. 120, 040405 (2018b).
  * Milz and Modi [2021] S. Milz and K. Modi, PRX Quantum 2, 030201 (2021).
  * Epstein  _et al._ [2014] J. M. Epstein, A. W. Cross, E. Magesan, and J. M. Gambetta, Phys. Rev. A 89, 062321 (2014).
  * Ball  _et al._ [2016] H. Ball, T. M. Stace, S. T. Flammia, and M. J. Biercuk, Phys. Rev. A 93, 022303 (2016).
  * Figueroa-Romero  _et al._ [2021] P. Figueroa-Romero, K. Modi, R. J. Harris, T. M. Stace, and M.-H. Hsieh, PRX Quantum 2, 040351 (2021).
  * Figueroa-Romero  _et al._ [2022] P. Figueroa-Romero, K. Modi, and M.-H. Hsieh, Quantum 6, 868 (2022).
  * Ceasura  _et al._ [2022] A. Ceasura, P. Iyer, J. J. Wallman, and H. Pashayan, Non-exponential behaviour in logical randomized benchmarking (2022), arXiv:2212.05488 [quant-ph] .
  * Chen  _et al._ [2020] Y.-Q. Chen, K.-L. Ma, Y.-C. Zheng, J. Allcock, S. Zhang, and C.-Y. Hsieh, Phys. Rev. Appl. 13, 034045 (2020).
  * White  _et al._ [2020] G. A. L. White, C. D. Hill, F. A. Pollock, L. C. L. Hollenberg, and K. Modi, Nat. Commun. 11, 6301 (2020).
  * Goswami  _et al._ [2021] K. Goswami, C. Giarmatzi, C. Monterola, S. Shrapnel, J. Romero, and F. Costa, Phys. Rev. A 104, 022432 (2021).
  * White  _et al._ [2022] G. White, F. Pollock, L. Hollenberg, K. Modi, and C. Hill, PRX Quantum 3, 020344 (2022).
  * Papič and de Vega [2022] M. Papič and I. de Vega, Phys. Rev. A 105, 022605 (2022).
  * Guo [2022] C. Guo, SciPost Phys. 13, 028 (2022).
  * White [2022] G. White, Nat. Rev. Phys. 4, 287 (2022).
  * Ezzell  _et al._ [2022] N. Ezzell, B. Pokharel, L. Tewala, G. Quiroz, and D. A. Lidar, Dynamical decoupling for superconducting qubits: a performance survey (2022), arXiv:2207.03670 [quant-ph] .
  * Winick  _et al._ [2022] A. Winick, J. J. Wallman, D. Dahlen, I. Hincks, E. Ospadov, and J. Emerson, Concepts and conditions for error suppression through randomized compiling (2022), arXiv:2212.07500 [quant-ph] .
  * Berk  _et al._ [2021] G. D. Berk, S. Milz, F. A. Pollock, and K. Modi, Extracting quantum dynamical resources: Consumption of non-Markovianity for noise reduction (2021), arXiv:2110.02613 [quant-ph] .
  * Ng  _et al._ [2011] H. K. Ng, D. A. Lidar, and J. Preskill, Phys. Rev. A 84, 012305 (2011).
  * Pokharel  _et al._ [2018] B. Pokharel, N. Anand, B. Fortman, and D. A. Lidar, Phys. Rev. Lett. 121, 220502 (2018).
  * Ware  _et al._ [2021] M. Ware, G. Ribeill, D. Ristè, C. A. Ryan, B. Johnson, and M. P. da Silva, Phys. Rev. A 103, 042604 (2021).
  * Hashim  _et al._ [2021] A. Hashim, R. K. Naik, A. Morvan, J.-L. Ville, B. Mitchell, J. M. Kreikebaum, M. Davis, E. Smith, C. Iancu, K. P. O’Brien, I. Hincks, J. J. Wallman, J. Emerson, and I. Siddiqi, Phys. Rev. X 11, 041039 (2021).
  * Kelly  _et al._ [2014] J. Kelly, R. Barends, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, I.-C. Hoi, E. Jeffrey, A. Megrant, J. Mutus, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White, A. N. Cleland, and J. M. Martinis, Phys. Rev. Lett. 112, 240504 (2014).
  * Souza  _et al._ [2015] A. M. Souza, R. S. Sarthour, I. S. Oliveira, and D. Suter, Phys. Rev. A 92, 062332 (2015).
  * Note [1] We point out that the average gate-fidelity of a quantum channel $\tilde{\mathcal{G}}$ with respect to a gate $\mathcal{G}$, is a measure of their “average orthogonality”, rather than their distinguishability (which albeit related, are not quite the same thing): $\mathrm{F}_{\mathsf{avg}}:=\intop\ilimits@{d}\psi\operatorname{tr}[\tilde{\mathcal{G}}(\psi)\mathcal{G}(\psi)]$, where $\psi$ are (uniformly distributed) pure states.
  * White  _et al._ [2023] G. A. L. White, K. Modi, and C. D. Hill, Phys. Rev. Lett. 130, 160401 (2023).
  * Milz  _et al._ [2020] S. Milz, F. Sakuldee, F. A. Pollock, and K. Modi, Quantum 4, 255 (2020).
  * Note [2] Non-Markovianity _can_ be said to introduce a type of gate-dependence; nevertheless, here by explicit gate-independence we mean noise associated to $\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G}$, where $\mathcal{I}_{\mathsf{E}}$ is an identity map on $\mathsf{E}$ and $\mathcal{G}$ is an ideal gate, is not explicitly dependent on $\mathcal{G}$.
  * Note [3] Here with noise being time-stationary, albeit non-Markovian, we mean that noise associated to $\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G}_{i}$, where $\mathcal{I}_{\mathsf{E}}$ is an identity channel on $\mathsf{E}$ and $\mathcal{G}_{i}$ is the ideal gate in $\mathsf{S}$ at timestep $i$, is independent of the timestep $i$.
  * Khodjasteh and Lidar [2007] K. Khodjasteh and D. A. Lidar, Phys. Rev. A 75, 062310 (2007).
  * Arenz  _et al._ [2017] C. Arenz, D. Burgarth, and R. Hillier, J. Phys. A: Math. Theor. 50, 135303 (2017).
  * Arenz  _et al._ [2018] C. Arenz, D. Burgarth, P. Facchi, and R. Hillier, J. Math. Phys. 59 (2018).
  * Burgarth  _et al._ [2021a] D. Burgarth, P. Facchi, M. Fraas, and R. Hillier, SciPost Phys. 11, 027 (2021a).
  * Qi  _et al._ [2023] J. Qi, X. Xu, D. Poletti, and H. K. Ng, Phys. Rev. A 107, 032615 (2023).
  * Wallman and Emerson [2016] J. J. Wallman and J. Emerson, Phys. Rev. A 94, 052325 (2016).
  * Sanders  _et al._ [2015] Y. R. Sanders, J. J. Wallman, and B. C. Sanders, New J. Phys. 18, 012002 (2015).
  * Note [4] The statement in Eq. (11) holds more generally for $\mathrm{f}_{m}$ being noisy expectation values.
  * Dirkse  _et al._ [2019] B. Dirkse, J. Helsen, and S. Wehner, Phys. Rev. A 99, 012315 (2019).
  * Note [5] In the $\mathsf{S}\mathsf{E}$ case, it is unclear whether $\mathsf{S}$-Pauli-twirling _indeed will always_ either only decrease or leave the total unitarity unchanged, due precisely to $\mathsf{E}$; here the role of the average sequence fidelity would also come into play to always decrease the variance.
  * Abdurakhimov  _et al._ [2022] L. V. Abdurakhimov, I. Mahboob, H. Toida, K. Kakuyanagi, Y. Matsuzaki, and S. Saito, PRX Quantum 3, 040332 (2022).
  * Addis  _et al._ [2015] C. Addis, F. Ciccarello, M. Cascio, G. M. Palma, and S. Maniscalco, New Journal of Physics 17, 123004 (2015).
  * D’Arrigo  _et al._ [2019] A. D’Arrigo, G. Falci, and E. Paladino, The European Physical Journal Special Topics 227, 2189 (2019).
  * Burgarth  _et al._ [2021b] D. Burgarth, P. Facchi, M. Ligabò, and D. Lonigro, Phys. Rev. A 103, 012203 (2021b).
  * Burgarth  _et al._ [2021c] D. Burgarth, P. Facchi, D. Lonigro, and K. Modi, Phys. Rev. A 104, L050404 (2021c).
  * Claes  _et al._ [2021] J. Claes, E. Rieffel, and Z. Wang, PRX Quantum 2, 010351 (2021).
  * Chruściński and Kossakowski [2006] D. Chruściński and A. Kossakowski, Phys. Rev. A 73, 062314 (2006).
  * Chiribella  _et al._ [2008] G. Chiribella, G. M. D’Ariano, and P. Perinotti, Phys. Rev. Lett. 101, 060401 (2008).
  * Chiribella  _et al._ [2009] G. Chiribella, G. M. D’Ariano, and P. Perinotti, Phys. Rev. A 80, 022339 (2009).
  * Portmann  _et al._ [2017] C. Portmann, C. Matt, U. Maurer, R. Renner, and B. Tackmann, IEEE Trans. Inf. Theory , 1–1 (2017).
  * Nurdin and Gough [2021] H. I. Nurdin and J. Gough, 2021 60th IEEE Conference on Decision and Control (CDC) 10.1109/cdc45484.2021.9683765 (2021).
  * Costa and Shrapnel [2016] F. Costa and S. Shrapnel, New J. Phys. 18, 063032 (2016).
  * Kretschmann and Werner [2005] D. Kretschmann and R. F. Werner, Phys. Rev. A 72, 062323 (2005).
  * Gutoski and Watrous [2007] G. Gutoski and J. Watrous,  STOC ’07, 565–574 (2007).
  * Milz  _et al._ [2017] S. Milz, F. A. Pollock, and K. Modi, Open Syst. Inf. Dyn. 24, 1740016 (2017).
  * Taranto [2020] P. Taranto, Int. J. Quantum Inf. 18, 1941002 (2020).
  * Milz  _et al._ [2019] S. Milz, M. S. Kim, F. A. Pollock, and K. Modi, Phys. Rev. Lett. 123, 040401 (2019).
  * Note [6] The Clifford group on $n$-qubits is defined as the set of unitaries normalizing the $n$-qubit Pauli group, $\mathds{P}_{n}^{*}$ modulo the identity, i.e., $\mathds{C}_{n}:=\\{G:\mathds{U}(2^{n})|P\in\pm\mathds{P}_{n}^{*}\,\Rightarrow\,GPG^{\dagger}\in\pm\mathds{P}_{n}^{*}\\}$. The crucial property for standard RB is that the $n$-qubit Clifford group forms a so-called unitary 2-design, while the main limitation to scalability comes from the fact that $n$-qubit Clifford gates are composite gates, with both the number of elements and elementary components scaling non-favorably in $n$.
  * Blume-Kohout  _et al._ [2022] R. Blume-Kohout, M. P. da Silva, E. Nielsen, T. Proctor, K. Rudinger, M. Sarovar, and K. Young, PRX Quantum 3, 020335 (2022).
  * Note [7] A unitary $t$-design is an ensemble of unitary matrices that reproduce up to the $t$th statistical moment of the unitary group with the uniform, so-called Haar measure; see e.g. [4].
  * Milz  _et al._ [2021] S. Milz, C. Spee, Z.-P. Xu, F. Pollock, K. Modi, and O. Gühne, SciPost Phys. 10 (2021).
  * Note [8] The only assumption we don’t explicitly relax is gate-dependence, i.e., the noise maps $\Lambda_{i}$ at any given time-step $i$ do not depend explicitly on which gate $\mathcal{G}_{i}$ was applied.
  * Kong [2021] L. Kong, A framework for randomized benchmarking over compact groups (2021), arXiv:2111.10357 [quant-ph] .
  * Note [9] A representation of a finite group $\mathds{G}$ can be defined as a map from $\mathds{G}$ to a vector space of unitary linear operators, e.g., that of $d$-dimensional complex matrices $\phi:\mathds{G}\to\mathds{U}(\mathds{C}_{d\times{d}})$. We call a representation reducible if it can be expressed as a direct sum of irreducible representations (irreps), and we say this representation is multiplicity-free if it contains no more than a single copy of such irreps: $\phi=\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigoplus$}}}}_{\mu}\phi_{\mu}^{\otimes{n_{\mu}}}$ with $n_{\mu}=1$ for all labels $\mu$. This is a technical constraint that can potentially be relaxed [63, 15] with some added but potentially non-crucial complications..
  * Note [10] We use the notation $\hat{\mathcal{P}}_{\pi}$ in order to not excessively conflate terms in Eq. (22), however, all $\hat{\mathcal{P}}_{\pi}$ here are literally projection operators, not stemming _by definition_ from a quantum channel. Of course one can still also associate a channel $\mathcal{P}_{\pi}$ to them, which however, does not necessarily has the action of a projection superoperator.
  * Helsen  _et al._ [2019] J. Helsen, X. Xue, L. M. K. Vandersypen, and S. Wehner, npj Quantum Inf. 5, 71 (2019).
  * Gough and Nurdin [2017] J. E. Gough and H. I. Nurdin, in  _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ (IEEE, 2017) pp. 6155–6160.
  * Note [11] Equivalently a unitary 2-design; the projectors are $\hat{\mathcal{P}}_{1}=\Psi$ and $\hat{\mathcal{P}}_{2}=\mathds{1}-\Psi$ where $\Psi:=\sum@\slimits@|ii\rangle\\!\langle{jj}|/d_{\mathsf{S}}$. See e.g. [64]. For details on the particular case of the Average Sequence Fidelity (ASF) in the Clifford case, see [25].
  * Sacchi [2005] M. F. Sacchi, Phys. Rev. A 71, 062340 (2005).
  * Wallman and Flammia [2014] J. J. Wallman and S. T. Flammia, New J. Phys. 16, 103032 (2014).
  * Kern  _et al._ [2005] O. Kern, G. Alber, and D. L. Shepelyansky, Eur. Phys. J. D. 32, 153–156 (2005).
  * Note [12] Strictly, the ordering in the vectorized (hat $\hat{\mathcal{X}}$) representation, as we defined it, is $\mathsf{E}\mathsf{S}\otimes\mathsf{E}\mathsf{S}$; we are slightly abusing notation for clarity of presentation.
  * Note [13] Notice that the time-stationary case, $\Lambda_{i}=\Lambda$ for all $i$, would leave the ASF in a _time-non-stationary_ noise form due to environment dependence, as we would drop all single $i$ indices, but the remaining $\mathsf{E}$ indices, $\mu,\sigma$, and thus all quality factors $p^{(\mu\sigma)}$ too, would remain distinct for each time-step.
  * Note [14] (Complete) Positivity of a quantum channel $\Lambda$ does not correspond to positive semi-definiteness of $\hat{\Lambda}$; however, it does correspond to positive semi-definiteness of the so-called Choi state of $\Lambda$, which can be defined as $\Upsilon_{\Lambda}:=(\Lambda\otimes\mathcal{I})\Psi$, where $\Psi=\sum@\slimits@|ii\rangle\\!\langle{jj}|$ is a maximally entangled operator on $\mathsf{S}^{\otimes 2}$.

## Glossary

##  Appendix A Preliminaries

###  A.1 Notation

Throughout the manuscript, we refer to an environment $\mathsf{E}$ and a
system $\mathsf{S}$: these are quantum systems with Hilbert spaces
$\mathscr{H}_{\mathsf{E}}$ and $\mathscr{H}_{\mathsf{S}}$, respectively, where
we take $\dim(\mathscr{H}_{\mathsf{S}}):=d_{\mathsf{S}}<\infty$. When
$\mathsf{E}$ is assumed to be finite, we denote
$\dim(\mathscr{H}_{\mathsf{E}}):=d_{\mathsf{E}}$. When we refer to the
composite $\mathsf{S}\mathsf{E}$, we mean the space
$\mathscr{H}_{\mathsf{E}}\otimes\mathscr{H}_{\mathsf{S}}$. When referring to
generic finite-dimensional Hilbert spaces $\mathscr{H}$, we normally use
$d=\dim(\mathscr{H})$.

We mainly use either the capital Greek letters $\Phi,\,\Lambda$, or a curly
font, e.g., $\mathcal{Q},\,\mathcal{P}$, to denote quantum channels. We use
the term quantum channel to mean Completely Positive (CP) map, and the term
quantum gate to mean a digital unitary map. One particular exception is the
symbol $\mathcal{F}_{m}$ which will be used specifically in the context of RB.
Finally, we use $\mathcal{I}$ for an identity channel and $\mathds{1}$ for an
identity operator.

We will employ a vectorized representation of density matrices (which we often
refer to simply as quantum states). Specifically, we employ

$\displaystyle|\cdot\operatorname*{\\!\rangle\\!\rangle}$
$\displaystyle:=\mathrm{vec}(\cdot),\quad\text{where}\quad\mathrm{vec}(|i\rangle\\!\langle{j}|)=|ij\rangle\quad\text{for
any}\quad\,|i\rangle,|j\rangle\in\mathscr{H},$ (12)

i.e., $|\cdot\operatorname*{\\!\rangle\\!\rangle}$ is a vectorization of
matrices sequentially stacking rows as columns, with the corresponding dual
being
$\operatorname*{\langle\\!\langle\\!}\cdot|:=|\cdot\operatorname*{\\!\rangle\\!\rangle}^{\dagger}$
where $\dagger$ is a conjugate transpose. Notice, thus, that we can now
express interior products as
$\operatorname*{\langle\\!\langle\\!}{A}|B\operatorname*{\\!\rangle\\!\rangle}=\operatorname{tr}[A^{\dagger}{B}]$.
Similarly then, for any quantum channel between bounded operator spaces,
$\Phi:\mathscr{B}(\mathscr{H})\to\mathscr{B}(\mathscr{H})$, with Kraus
operators $\phi_{\mu}$ (where $\mu$ can run on up to $d^{2}$ terms), we define

$\displaystyle|\Phi(\cdot)\operatorname*{\\!\rangle\\!\rangle}$
$\displaystyle:=\hat{\Phi}|\cdot\operatorname*{\\!\rangle\\!\rangle},\qquad\text{where}\qquad\hat{\Phi}=\sum_{\mu}\phi_{\mu}\otimes\phi_{\mu}^{*},$
(13)

i.e., $\hat{\cdot}$ turns quantum channels into matrices. Here we will always
assume _at least_ a trace non-increasing property
$\sum\phi_{\mu}^{\dagger}\phi_{\mu}\leq\mathds{1}$ (meaning the eigenvalues of
such sum of Kraus operators are positive and upper-bounded by unity) on
quantum channels and otherwise mention when a full Trace Preserving (TP)
property, for which $\sum\phi_{\mu}^{\dagger}\phi_{\mu}=\mathds{1}$, is being
considered.

###  A.2 Non-Markovianity

Here we briefly define formally what we mean by the term _non-Markovianity_ ,
as summarized in [25]. A comprehensive review and discussion on the topic can
be seen in [21].

Within classical probability theory, non-Markovianity can be described by a
stochastic process, $\\{X_{t}\\}$, where generally probabilities are
conditionally dependent on the past, i.e.

$\mathds{P}(x_{k}|x_{k-1},\ldots,x_{0})=\mathds{P}(x_{k}|x_{k-1},\ldots,x_{k-\ell}),$
(14)

for any integers $0\leq{\ell}\leq{k}$ and sequences of event outcomes $x_{i}$,
with $\mathds{P}(\cdot|\cdot)$ denoting a conditional probability. In
particular, when $\ell=1$, the process is called Markovian and when $\ell=0$
it is called random; otherwise, the process is non-Markovian with Markov order
$\ell$. The fact that the complexity in describing non-Markovian processes
increases exponentially in increasing Markov-order can be seen from joint
probabilities requiring up to $\ell$-point correlations within the respective
conditional probabilities.

Quantum mechanically, the process tensor (a.k.a. quantum comb [65, 66], causal
box [67], correlation kernel [68], process matrix [69], channel with memory
[70], or strategy [71], to mention some) framework [19, 20, 72, 45, 73, 21,
74], takes into account the invasive nature of observation to unambiguously
provide a generalization of the condition in Eq. (14), as shown in [20]. In
this case, the medium for information to be sent across timesteps is an
environment $\mathsf{E}$, part of a bipartite system $\mathsf{S}\mathsf{E}$,
with $\mathsf{S}$ being the system of interest. Then, for an initial state
$\rho$ of $\mathsf{S}\mathsf{E}$, and upon measuring a Positive Operator
Valued Measurement (POVM) $\mathcal{J}_{k}:=\\{M_{x_{n}}^{(k)}\\}_{x_{n}}$ on
system $\mathsf{S}$, we may describe the probability of observing a sequence
of quantum events $x_{k},\ldots,x_{0}$ by

$\mathds{P}(x_{k},\ldots,x_{0}|\mathcal{J}_{k},\ldots,\mathcal{J}_{0}):=\operatorname{tr}\left[M_{x_{k}}\,\rho^{(k)}\right],$
(15)

where
$\rho^{(k)}:=\operatorname{tr}_{\mathsf{E}}\left[\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{i=1}^{k-1}\left(\mathcal{U}_{i}\circ\mathcal{A}_{x_{i}}\right)\rho\right]$
is the state of system $\mathsf{S}$ at the $k$th timestep, with
$\mathcal{U}_{i}$ being dynamical maps on $\mathsf{S}\mathsf{E}$ describing
the evolution of the full system-environment between any two timestep, and
$\mathcal{A}_{x_{i}}$ being CP maps acting on system $\mathsf{S}$ alone at
timestep $i$: precisely, each
$\mathcal{J}_{i}:=\\{\mathcal{A}^{(i)}_{x_{n}}\\}_{x_{n}}$ is called an
instrument, where $\mathcal{A}^{(i)}_{x_{n}}$ is an experimental intervention
represented by a CP map with state outcome $x_{n}$, and such that
$\sum_{x_{n}}\mathcal{A}^{(i)}_{x_{n}}=\mathcal{A}^{(i)}$ is a Completely
Positive Trace Preserving (CPTP) map.

We drop the super-indices in Eq. (15) for clarity, which we may write more
succinctly as the inner product

$\displaystyle\mathds{P}(x_{k},\ldots,x_{0}|\mathcal{J}_{k},\ldots,\mathcal{J}_{0})=\operatorname{tr}\left(\Upsilon_{k}\,\Theta_{k}^{\mathrm{T}}\right),$
(16)

where $\mathrm{T}$ denotes a transpose, and $\Upsilon_{k}$ and $\Theta_{k}$
are tensors containing all dynamics $\\{\mathcal{U}_{i}\\}$ and all
interventions $\\{\mathcal{A}_{i}\\}$; in the Choi-Jamiołkowski
representation, these take the form

$\displaystyle\Upsilon_{k}$
$\displaystyle:=\operatorname{tr}_{\mathsf{E}}\left\\{\left[\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{i=1}^{k}(\mathcal{U}_{i}\otimes\mathcal{I}_{\mathsf{aux}}\circ\mathscr{S}_{i})\right]\rho\otimes\psi^{\otimes{k}}\right\\},$
(17) $\displaystyle\Theta_{k}$
$\displaystyle:=M_{x_{k}}\otimes\left[\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{i=1}^{k-1}\left(\mathds{1}_{\mathsf{A}_{i}}\otimes\mathcal{A}_{x_{i}}\right)\right]\psi^{\otimes{k}},$
(18)

where
$\mathsf{aux}:=\mathsf{A}_{1}\mathsf{B}_{1}\ldots\mathsf{A}_{k}\mathsf{B}_{k}$,
with $\mathsf{A}_{i}$, $\mathsf{B}_{i}$ being $d_{\mathsf{S}}$-dimensional
auxiliary spaces, $\mathscr{S}_{i}$ being a swap map between $\mathsf{S}$ and
$\mathsf{A}_{i}$, and $\psi=\sum|ii\rangle\\!\langle{jj}|$ being an
unnormalized maximally entangled state.

The process tensor framework thus allows to neatly separate the underlying
dynamical source for any given quantum process, including all temporal
correlations therein, from all experimentally controllable operations. This
description is entirely general as a quantum stochastic process framework [45,
21], and similarly the instruments used to describe interventions are entirely
general and can be temporally correlated themselves. Similar to the case of
quantum states, the choice of employing a Choi state representation in Eq.
(17) allows us to readily deduce properties of the process. In particular,
temporal properties get codified as spatial properties within the Choi state,
so that a Markovian process takes an uncorrelated form,

$\Upsilon^{(\mathrm{M})}:=\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{i}\mathscr{Y}_{i:i-1}\otimes\rho_{\mathsf{S}},$
(19)

with $\mathscr{Y}_{i:i-1}$ being individual Choi states of dynamics connecting
the $(i-1)$th and $i$th steps. This implies that we may quantify the non-
Markovianity of a process by simply quantifying its distinguishability from
the closest Markovian one, i.e.,
$\mathcal{N}:=\min_{\Upsilon^{(\mathrm{M})}}d(\Upsilon,\Upsilon^{(\mathrm{M})})$,
for any operationally meaningful distance measure $d(\cdot,\cdot)$. Here we
won’t focus on the multi-time processes themselves, either $\Upsilon$ or
$\Upsilon^{(\mathrm{M})}$, but rather on the consequences of each kind within
RB.

###  A.3 The Randomized Benchmarking protocol and the standard Markovian
analysis

A standard Randomized Benchmarking (RB) protocol on a system of interest,
$\mathsf{S}$, of dimension $d_{\mathsf{S}}$ proceeds as follows:

  1. 1.

Prepare an initial state $\rho$ on $\mathsf{S}$.

  2. 2.

Sample a sequence of $m$ distinct quantum gates,
$\mathcal{G}_{1},\mathcal{G}_{2},\ldots,\mathcal{G}_{m}$, uniformly at random
from a given gate set $\mathds{G}$ forming a finite group, and compile the
inverse of the sequence in
$\mathcal{G}_{m+1}=\mathcal{G}_{1}^{-1}\circ\cdots\circ\mathcal{G}_{m}^{-1}$,
where here $\circ$ denotes composition of maps (i.e.,
$\mathcal{A}\circ\mathcal{B}$ reading $\mathcal{A}$ _after_ $\mathcal{B}$). We
refer to $\mathcal{G}_{m+1}$ as an undo-gate.

  3. 3.

Apply the composition
$\mathcal{S}_{m}:=\mathcal{G}_{m+1}\circ\mathcal{G}_{m}\circ\cdots\circ\mathcal{G}_{1}$
on $\rho$.

  4. 4.

Estimate the probability
$f_{m}:=\operatorname{tr}\left[M\mathcal{S}_{m}\left(\rho\right)\right]$ via a
POVM element $M$. Generally, a real physical implementation of these steps
renders $f_{m}\neq\operatorname{tr}(M\rho)$.

  5. 5.

Repeat $n$ times the steps 1 to 4 for the same sequence length $m$, same
initial state $\rho$, same POVM element $M$, and different sets of gates
chosen uniformly at random,
$\left\\{\mathcal{G}_{i}^{(1)}\right\\}_{i=1}^{m},\,\left\\{\mathcal{G}_{i}^{(2)}\right\\}_{i=1}^{m},\ldots,\left\\{\mathcal{G}_{i}^{(n)}\right\\}_{i=1}^{m}$
from $\mathds{G}$, to obtain the corresponding probabilities
$f_{m}^{(1)},f_{m}^{(2)},\ldots{f}_{m}^{(n)}$. Estimate the average
probability $\mathcal{F}_{m}=1/n\sum_{i=1}^{n}f_{m}^{(i)}$. We refer to
$\mathcal{F}_{m}$ as an Average Sequence Fidelity (ASF).

  6. 6.

Examine the behavior of the estimated ASF, $\mathcal{F}_{m}$, over different
sequence lengths $m$, with respect to the ideal expected
$\operatorname{tr}(M\rho)$.

The protocol up to this point is independent of any assumptions for the
underlying sources of noise causing the outputs to deviate from the ideal
expected one, as well as the choices for the gate-set (which is only
restricted to be a finite group, hence contain inverses), the dimension of the
system $\mathsf{S}$ and the initial state and final measurement.

For several practical reasons, the system $\mathsf{S}$ usually consists of one
or two-qubits, i.e., $d_{\mathsf{S}}=2^{n}$ for either $n=1$ or $n=2$, the
initial state is taken to be the ground state in the computational basis
$\rho=|0\rangle\\!\langle 0|$ and $M$ the corresponding projector
$M=|0\rangle\\!\langle 0|$, so that the ASF is rather simply a survival
probability of the ground state.

###  A.4 Modeling for Markovian RB

To make the analysis tractable and connect the RB protocol with the estimation
of so-called average gate fidelities, in a manner robust to SPAM errors, some
assumptions are commonly made for the ideal gate-set and the underlying noise
affecting the physical protocol: the ideal gate-set $\mathds{G}$ is usually
taken to be the multi-qubit Clifford group 666 The Clifford group on
$n$-qubits is defined as the set of unitaries normalizing the $n$-qubit Pauli
group, $\mathds{P}_{n}^{*}$ modulo the identity, i.e.,
$\mathds{C}_{n}:=\\{G:\mathds{U}(2^{n})|P\in\pm\mathds{P}_{n}^{*}\,\Rightarrow\,GPG^{\dagger}\in\pm\mathds{P}_{n}^{*}\\}$.
The crucial property for RB is that the $n$-qubit Clifford group forms a so-
called unitary 2-design, while the main limitation to scalability comes from
the fact that $n$-qubit Clifford gates are composite gates, with both the
number of elements and elementary components scaling non-favorably in $n$. and
noise assumed to be at least _Markovian_ , i.e., uncorrelated in time, and
thus associated to each gate independently. Moreover, noise is also commonly
assumed to be independent of both _when_ a gate is applied (time-stationary)
and _which_ gate is applied (gate-independent). These can be made
mathematically concrete by modeling individual noisy gates as
$\tilde{\mathcal{G}}:=\Lambda\circ\mathcal{G}$, or equivalently, defining
noise quantum channels by $\Lambda:=\tilde{\mathcal{G}}\circ\mathcal{G}^{-1}$,
where $\tilde{\mathcal{G}}$ and $\Lambda$ are in general CP maps sending
inputs in $\mathsf{S}$ to outputs in $\mathsf{S}$. Notice that, physically,
here we are assuming that the ideal digital gates $\mathcal{G}$ are precisely
instantaneous noiseless quantum gates, while the maps $\Lambda$ can model _the
effects_ of any possible Markovian quantum noise process. Furthermore, we
could have equivalently defined $\tilde{\mathcal{G}}=\mathcal{G}\circ\Lambda$
and we would be modeling the same phenomena. Another, non-equivalent choice,
could be to model noisy gates with dynamical generators (Hamiltonians) so that
errors occur _during_ a physical (e.g., a simple pulse) but error-free gate;
this, albeit physically plausible, has limitations: a deeper discussion on
Markovian errors and their modeling can be seen in [76].

Figure 7: Modeling Markovian noise in a quantum circuit: An initial state
$\rho$, is acted on with an ideal gate $\mathcal{G}_{1}$, which itself can be
affected by a quantum channel $\Lambda_{1}$; subsequently, an ideal gate
$\mathcal{G}_{2}$ is applied and can be affected by another quantum channel
$\Lambda_{2}$, and so on, until a measurement with an element $M$ is
performed. Modeling of noisy gates is described by the composition
$\Lambda_{i}\circ\mathcal{G}_{i}$, or equivalently,
$\mathcal{G}_{i}\circ\Lambda_{i}$, for any CP and trace non-increasing map
$\Lambda_{i}$. When $\Lambda_{i}=\Lambda_{j}=\Lambda$ for all $i\neq{j}$ time-
step labels, we say that noise in the sequence is time-stationary.

With these assumptions of Markovianity, time-independence and gate-
independence, we can model noisy initial states as
$\tilde{\rho}:=\Lambda_{0}(\rho)$, RB sequences as
$\tilde{\mathcal{S}}_{m}:=\tilde{\mathcal{G}}_{m+1}\circ\tilde{\mathcal{G}}_{m}\circ\cdots\circ\tilde{\mathcal{G}}_{1}$,
and noisy probabilities as
$f_{m}:=\operatorname*{\langle\\!\langle\\!}{M}|\hat{\tilde{\mathcal{S}}}_{m}|\tilde{\rho}\operatorname*{\\!\rangle\\!\rangle}$.
Notice that measurement noise is indistinguishable from that of the initial
state. Then, this allows to show that the analytical ASF, i.e., the analytical
average of the sequences $f_{m}$ over uniformly distributed Clifford gates (or
indeed over any so-called unitary 2-design 777A unitary $t$-design is an
ensemble of unitary matrices that reproduce up to the $t$th statistical moment
of the unitary group with the uniform, so-called Haar measure; see e.g. [4])
will take an exponential form in the sequence length $m$,

$\mathcal{F}_{m}=Ap^{m}+B,\quad\text{where}\quad\,p=\frac{d_{\mathsf{S}}\mathfrak{F}_{\Lambda}-1}{d_{\mathsf{S}}-1},\quad\text{with}\quad\mathfrak{F}_{\Lambda}:=\frac{d_{\mathsf{S}}+\operatorname{tr}\left(\hat{\Lambda}\right)}{d_{\mathsf{S}}(d_{\mathsf{S}}+1)},$
(20)

and $0\leq{A,B}\leq 1$ being constants capturing SPAM errors,

$A=\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Lambda}^{\prime}|\rho_{\mathsf{S}}-\mathds{1}/d_{\mathsf{S}}\operatorname*{\\!\rangle\\!\rangle}=\operatorname{tr}[M\Lambda^{\prime}(\rho_{\mathsf{S}}-\mathds{1}/d_{\mathsf{S}})],\qquad{B}=\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Lambda}^{\prime}|\mathds{1}/d_{\mathsf{S}}\operatorname*{\\!\rangle\\!\rangle}=\operatorname{tr}[M\Lambda^{\prime}(\mathds{1}/d_{\mathsf{S}})],$
(21)

where $\Lambda^{\prime}=\Lambda\circ\Lambda_{0}$. In the case
$\rho_{\mathsf{S}}=M=|0\rangle\\!\langle 0|$, then $A=\langle
0|\Lambda^{\prime}(|0\rangle\\!\langle 0|)|0\rangle-B$, and $B=\langle
0|\Lambda^{\prime}(\mathds{1})|0\rangle/d_{\mathsf{S}}$, so that $B$ can be
seen as quantifying non-unitality of $\Lambda^{\prime}$ and $A$ as a
corresponding fidelity quantifier minus such unitality, so that for _low_ SPAM
noise (unital and close to identity), both $A,B\approx 1/d_{\mathsf{S}}$.

Detail of this result can be seen in standard RB literature, such as [4]. The
quantity ${p}\leq 1$ is sometimes labeled the _noise strength_ of $\Lambda$,
while $\mathfrak{F}_{\Lambda}\leq 1$ is the so-called _average gate-fidelity_
of $\Lambda$ with respect to the identity $\mathcal{I}$. Notice that
$\mathfrak{F}_{\Lambda}$ can equivalently be read as the average gate-fidelity
of any noisy gate, $\tilde{\mathcal{G}}_{i}$, with respect to the respective
ideal gate, $\mathcal{G}_{i}$: generically, the average gate-fidelity is
defined as
$\mathfrak{F}_{\Lambda}:=\mathbf{E}_{\psi}\left[\langle\psi|\Lambda(|\psi\rangle\\!\langle\psi|)|\psi\rangle\right]$,
where here $\mathbf{E}_{\psi}$ means uniform average over all pure states
$|\psi\rangle$ on $\mathsf{S}$, which can be shown to be equivalent to the
definition in Eq. (20). The noise strength is rigorously lower-bounded by
$1/\left(1-d_{\mathsf{S}}^{2}\right)$ and the average gate-fidelity by
$1/(d_{\mathsf{S}}+1)$, but this would be of course not sensible in practice
and $p$ is at least assumed to be positive.

Many generalizations and different versions of RB estimating various figures
of merit exist, and in particular many of the assumptions before can be
relaxed [15]. Here we will focus on the _non-Markovian_ case, i.e., where the
noise channels $\Lambda$ can be correlated with each other in time, and where
furthermore, the character of these correlations is allowed to be quantum
mechanical [78].

##  Appendix B Randomized Benchmarking for non-Markovian noise and the
reduction to the Markovian case

We now move on to considering a system-environment $\mathsf{S}\mathsf{E}$
composite, where we assume gates can only be applied directly in system
$\mathsf{S}$, and similarly state-preparation and measurement can solely be
done on $\mathsf{S}$. Following up on the definitions made in Appendix A.3 for
the Markovian case, we can generalize this to (temporally-correlated) non-
Markovian noise, by letting the noise channels act on the whole
$\mathsf{S}\mathsf{E}$ composite. That is, we now take $\Lambda$ to be a
quantum channel from $\mathsf{S}\mathsf{E}$ quantum states to
$\mathsf{S}\mathsf{E}$ quantum states, which can either be a unitary channel
or a CP map, depending on whether $\mathsf{S}\mathsf{E}$ is a closed system or
there is dissipation to a larger environment, respectively. The environment
$\mathsf{E}$ can thus be now understood as a quantum memory allowing for
quantum information to propagate between different instances of $\Lambda$. We
can further relax the time-stationary noise constraint, and associate a
distinct $\Lambda_{i}$ (acting on $\mathsf{S}\mathsf{E}$) to each ideal gate
$\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G}_{i}$ where here explicitly
$\mathcal{G}_{i}$ acts on $\mathsf{S}$ alone. Notice that, importantly, we
cannot associate noise to each individual gate $\mathcal{G}_{i}$ separately,
but we have to take into account the environment, and as we don’t have direct
access to it, we would need to track how it changes all intermediate noisy
gates throughout a given sequence. Many in-depth reviews about non-
Markovianity exist, and the topic is vast enough to be dealt with here, but
the main notions and formalism we employ can be consulted in [21]. Notice that
the RB protocol, of course, remains the same regardless of the underlying type
of noise, it is just our modeling that changes and gives a more general
analytical behavior of the ASF.

Hence we now model the global initial noisy state as
$\tilde{\rho}:=\Lambda_{0}(\rho_{\mathsf{E}}\otimes\rho_{\mathsf{S}})$ for
some ideal prepared state $\rho_{\mathsf{S}}$, and some fiducial state of the
environment $\rho_{\mathsf{E}}$. Similarly, now individual noisy gates at
time-step $i$ are taken to $\Lambda_{i}\circ\mathcal{G}_{i}$, where the
identity map on $\mathsf{E}$ is implied (henceforth, we do this when clear by
context). The final measurement remains as an operator $M$ acting solely on
$\mathsf{S}$. An ASF modeling a RB experiment under non-Markovian 888 The only
assumption we don’t explicitly relax is gate-dependence, i.e., the noise maps
$\Lambda_{i}$ at any given time-step $i$ do not depend explicitly on which
gate $\mathcal{G}_{i}$ was applied. noise would then take the form
$\mathcal{F}_{m}=\mathbf{E}_{\mathds{G}}\operatorname*{\langle\\!\langle\\!}{M}|\hat{\operatorname{tr}}_{\mathsf{E}}\hat{\Lambda}_{m+1}\hat{\mathcal{G}}_{m+1}\hat{\Lambda}_{m}\hat{\mathcal{G}}_{m}\cdots\hat{\Lambda}_{1}\hat{\mathcal{G}}_{1}|\tilde{\rho}\operatorname*{\\!\rangle\\!\rangle}$;
the average over the gate-set, $\mathbf{E}_{\mathds{G}}$, can then be computed
analytically at least either when $\mathds{G}$ is the multi-qubit Clifford
group (generally any unitary 2-design) [24] or when it is a finite group [25]
(extensions beyond finite groups or groups in general, have also been
proposed, at least within the Markovian case [80, 14]).

Figure 8: Modeling non-Markovian noise in a quantum circuit: Similar to the
case of Fig. 7, a sequence of ideal gates
$\mathcal{G}_{1},\mathcal{G}_{2},\ldots,\mathcal{G}_{m+1}$ is applied on a
system $\mathsf{S}$ with initial state $\rho_{\mathsf{S}}$ and final
measurement operator $M$. Now, however, noise is modeled as a CP $\Lambda_{i}$
associated to maps $\mathcal{I}_{\mathsf{E}}\otimes\mathcal{G}_{i}$, where
$\mathcal{I}_{\mathsf{E}}$ is an identity channel on a quantum environment
(i.e. an external, in principle inaccessible Hilbert space) $\mathsf{E}$,
initialized in some fiducial state $\rho_{\mathsf{E}}$. Temporal correlations,
i.e., non-Markovianity, between noise is mediated by $\mathsf{E}$, and SPAM
terms, $\Lambda_{0}$ and $\Lambda_{m+1}$ can correlate either preparation or
measurement on $\mathsf{S}$ and $\mathsf{E}$.

Explicitly, when $\mathds{G}$ is a finite group admitting a multiplicity-free
representation 999 A representation of a finite group $\mathds{G}$ can be
defined as a map from $\mathds{G}$ to a vector space of unitary linear
operators, e.g., that of $d$-dimensional complex matrices
$\phi:\mathds{G}\to\mathds{U}(\mathds{C}_{d\times{d}})$. We call a
representation reducible if it can be expressed as a direct sum of irreducible
representations (irreps), and we say this representation is multiplicity-free
if it contains no more than a single copy of such irreps:
$\phi=\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigoplus$}}}}_{\mu}\phi_{\mu}^{\otimes{n_{\mu}}}$
with $n_{\mu}=1$ for all labels $\mu$. This is a technical constraint that can
potentially be relaxed [63, 15] with some added but potentially non-crucial
complications., the analytical ASF for a RB experiment with $m$ gates drawn
uniformly from $\mathds{G}$, under gate-independent non-Markovian noise, takes
the form ([25])

$\mathcal{F}_{m}=\sum_{\pi\in{R}_{\mathds{G}}}\\!\operatorname*{\langle\\!\langle\\!}{M}|\hat{\operatorname{tr}}_{\mathsf{E}}\hat{\Lambda}_{m+1}\left(\hat{\mathcal{Q}}_{m,\pi}\otimes\hat{\mathcal{P}}_{\pi}\right)|\tilde{\rho}\operatorname*{\\!\rangle\\!\rangle},$
(22)

where $R_{\mathds{G}}$ is a set of labels for the spaces associated to each
irreducible representation of $\mathds{G}$, the $\hat{\mathcal{P}}_{\pi}$ are
projector operators onto these 101010 We use the notation
$\hat{\mathcal{P}}_{\pi}$ in order to not excessively conflate terms in Eq.
(22), however, all $\hat{\mathcal{P}}_{\pi}$ here are literally projection
operators, not stemming _by definition_ from a quantum channel. Of course one
can still also associate a channel $\mathcal{P}_{\pi}$ to them, which however,
does not necessarily has the action of a projection superoperator., and
$\mathcal{Q}_{m,\pi}$ is a so-called length-$m$ quality map associated to the
sequence of noise maps on the $\pi$th irreducible subspace. We define
explicitly what quality maps are below, in Eq. (23), but their qualitative
significance is essentially a generalization of the noise strength $p$, in Eq.
(20), from the standard Clifford case: indeed for the Markovian noise case,
quality _maps_ are really just quality _scalar factors_ , and with the
Clifford group there are just two irreducible subspaces, with one quality
factor corresponding to 1 and the other to $p$ (these concepts, including the
definition of a quality factor, are explained e.g., in [83]).

Noise maps $\mathcal{Q}_{m,\pi}$ can be written in superoperator form as

$\displaystyle\hat{\mathcal{Q}}_{m,\pi}=\sum_{\\{\epsilon_{j},\epsilon_{j}^{\prime}\\}_{j=0}^{m}}^{d_{\mathsf{E}}}\frac{\prod_{i=m-1}^{0}\operatorname{tr}\left(\langle\epsilon_{i+1}\epsilon_{i+1}^{\prime}|\hat{\Lambda}_{i+1}|\epsilon_{i}\epsilon_{i}^{\prime}\rangle\hat{\mathcal{P}}_{\pi}\right)}{\operatorname{tr}\left(\hat{\mathcal{P}}_{\pi}\right)^{m}}|\epsilon_{m}\epsilon_{m}^{\prime}\rangle\\!\langle\epsilon_{0}\epsilon_{0}^{\prime}|,$
(23)

where here all $\\{|\epsilon_{i}\rangle,|\epsilon_{i}^{\prime}\\}_{i=0}^{m}$
are basis vectors of $\mathsf{E}$. The derivation of this expression for
quality maps can be seen in detail in [25], specifically in Appendices B.1 and
B.3.

For the particular case of the Clifford group cited in the main text, also a
full detail of derivation can be seen in Appendix B.2 of [25]: in this case,
there are two irreducible subspaces with $\hat{\mathcal{P}}_{1}=\Psi$ and
$\hat{\mathcal{P}}_{2}=\mathds{1}-\Psi$ where
$\Psi=\sum_{i,j=1}^{d_{\mathsf{S}}}|ii\rangle\\!\langle{jj}|/d_{\mathsf{S}}$
and $\mathds{1}$ is a $d_{\mathsf{S}}^{2}$ identity, so there are two
corresponding quality factors, $\mathcal{Q}_{m,1}$ and $\mathcal{Q}_{m,2}$,
which are explicitly defined through

$\displaystyle\mathcal{Q}_{m,1}\otimes\mathcal{P}_{1}$
$\displaystyle:=\left[(\$_{\Lambda}-\Theta_{\Lambda})^{\circ{m}}\otimes\mathcal{I}_{\mathsf{S}}\right]\left(\rho-\rho_{\mathsf{E}}\otimes\frac{\mathds{1}}{d_{\mathsf{S}}}\right),\quad\text{and}\quad\mathcal{Q}_{m,2}\otimes\mathcal{P}_{2}:=\Theta_{\Lambda}^{\circ{m}}(\rho_{\mathsf{E}})\otimes\frac{\mathds{1}}{d_{\mathsf{S}}},$
(24)

where

$\displaystyle\Theta_{\Lambda}(\cdot)=\operatorname{tr}_{\mathsf{S}}\left[\Lambda\left(\cdot\otimes\frac{\mathds{1}}{d_{\mathsf{S}}}\right)\right],\quad\text{and}\quad\$_{\Lambda}(\cdot):=\sum_{\mu}\operatorname{tr}_{\mathsf{S}}(\lambda_{\mu})(\cdot)\operatorname{tr}_{\mathsf{S}}(\lambda_{\mu}^{\dagger}),$
(25)

with $\lambda_{\mu}$ being the Kraus operators of $\Lambda$, i.e.
$\hat{\Lambda}=\sum_{\mu}\lambda_{\mu}\otimes\lambda_{\mu}^{*}$. In the main
text, we express this in a simplified way as

$\displaystyle\mathcal{Q}^{\prime}_{m,A}=\Lambda_{m+1}\circ(\mathcal{Q}_{m,1}\otimes\mathcal{P}_{1}),\quad\text{and}\quad\mathcal{Q}^{\prime}_{m,B}=\Lambda_{m+1}\circ(\mathcal{Q}_{m,2}\otimes\mathcal{P}_{2}).$
(26)

Here, however, we are not interested in studying general quality maps, but
rather whether we can operationally reduce them to the Markovian quality
factors, i.e., if we can _Markovianize_ the non-Markovian ASF of Eq. (22)
through an experimental protocol embedded in RB.

Thus, before considering the question of how to do this operationally, let us
show that indeed, if the noise across all steps is uncorrelated on
$\mathsf{S}\mathsf{E}$, so that no information leaks to $\mathsf{E}$ and then
gets propagated in time, we recover the Markovian ASF. That is, let us
consider noise of the form $\Lambda=\Theta\otimes\Phi$, where $\Theta$ is a
channel acting solely from and to $\mathsf{E}$, while $\Phi$ is a channel
between $\mathsf{S}$ systems, then

$\displaystyle\hat{\mathcal{Q}}_{m,\pi}$
$\displaystyle=\sum_{\\{\epsilon_{j},\epsilon_{j}^{\prime}\\}_{j=0}^{m}}^{d_{\mathsf{E}}}\frac{\prod_{i=m-1}^{0}\langle\epsilon_{i+1}\epsilon_{i+1}^{\prime}|\hat{\Theta}_{i+1}|\epsilon_{i}\epsilon_{i}^{\prime}\rangle\operatorname{tr}\left(\hat{\Phi}_{i+1}\hat{\mathcal{P}}_{\pi}\right)}{\operatorname{tr}\left(\hat{\mathcal{P}}_{\pi}\right)^{m}}|\epsilon_{m}\epsilon_{m}^{\prime}\rangle\\!\langle\epsilon_{0}\epsilon_{0}^{\prime}|$
$\displaystyle=\frac{\operatorname{tr}\left(\hat{\Phi}_{1}\hat{\mathcal{P}}_{\pi}\right)\cdots\operatorname{tr}\left(\hat{\Phi}_{m}\hat{\mathcal{P}}_{\pi}\right)}{\operatorname{tr}\left(\hat{\mathcal{P}}_{\pi}\right)^{m}}\,\prod_{i=m-1}^{0}\hat{\Theta}_{i+1},$
(27)

for all $m$, and since
$\hat{\operatorname{tr}}\,\hat{\Theta}_{m+1}\hat{\Theta}_{m}\cdots\hat{\Theta}_{1}\hat{\Theta}_{0}|\rho_{\mathsf{E}}\operatorname*{\\!\rangle\\!\rangle}=\operatorname{tr}\left[\Theta_{m+1}\circ\cdots\circ\Theta_{0}(\rho_{\mathsf{E}})\right]=1$,
it follows that

$\displaystyle\mathcal{F}_{m}$
$\displaystyle=\sum_{\pi\in{R}_{\mathds{G}}}\prod_{i=1}^{m}\left(\frac{\operatorname{tr}\hat{\Phi}_{i}\hat{\mathcal{P}}_{\pi}}{\operatorname{tr}\hat{\mathcal{P}}_{\pi}}\right)\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Phi}_{m+1}\hat{\mathcal{P}}_{\pi}\hat{\Phi}_{0}|\rho_{\mathsf{S}}\operatorname*{\\!\rangle\\!\rangle}$
$\displaystyle=\sum_{\pi\in{R}_{\mathds{G}}}f_{1,\pi}f_{2,\pi}\cdots
f_{m,\pi}\operatorname*{\langle\\!\langle\\!}{M}^{\prime}|\hat{\mathcal{P}}_{\pi}|\rho_{\mathsf{S}}^{\prime}\operatorname*{\\!\rangle\\!\rangle},$
(28)

where we defined

$f_{i,\pi}:=\frac{\operatorname{tr}\left(\hat{\Phi}_{i}\hat{\mathcal{P}}_{\pi}\right)}{\operatorname{tr}\hat{\mathcal{P}}_{\pi}},$
(29)

which are precisely the Markovian _quality factors_ [83] mentioned before, and
where $M^{\prime}$ and $\rho_{\mathsf{S}}^{\prime}$ are the SPAM-affected
measurement and initial state. This is precisely the Markovian generalization
of Eq. (20) under time-non-stationary noise and with the gate-set being any
finite group; indeed it can be shown that it reduces to it in the time-
stationary noise case and multi-qubit Clifford group [25].

We can now proceed to show that this _Markovianization_ effect can be
achieved, approximately, by interleaving so-called ideal UDD sequences into a
non-Markovian RB experiment.

##  Appendix C Non-Markovian Randomized Benchmarking with Ideal Pulse
Universal Dynamical Decoupling

DD is in general a quantum control technique (i.e., one involving physical
control via pulses applied on a control Hamiltonian) allowing to decouple a
system of interest from undesired degrees of freedom in its evolution; as
such, it can be seen as an essential quantum error suppression technique. A
digital, so-called _ideal_ DD scheme, assumes that the control pulses applied
on system $\mathsf{S}$ are infinitely strong and narrow; i.e., that they are
achieved by a time-non-stationary Hamiltonian of the form
$H_{\mathsf{ctrl}}^{(i)}(t)=\frac{\pi}{2}\delta(t-t_{0})v_{i}$, where $\delta$
is a Dirac delta, $t_{0}$ is the time at which the pulse is applied, and
$v_{i}$ is a system $\mathsf{S}$ operator element of a group $\mathds{V}$ such
that

$\frac{1}{|\mathds{V}|}\sum_{i=1}^{|\mathds{V}|}\mathcal{V}_{i}(Z)=W\otimes\mathds{1}_{\mathsf{S}},$
(30)

for any operator $Z$ in $\mathsf{S}\mathsf{E}$ and some operator $W$ on
$\mathsf{E}$ alone, where we defined

$\mathcal{V}_{i}(\cdot):=(\mathds{1}_{\mathsf{E}}\otimes{v}_{i})(\cdot)(\mathds{1}_{\mathsf{E}}\otimes{v}_{i}^{\dagger}).$
(31)

We only consider here $\mathds{V}$ to be a unitary group of dimension
$d_{\mathsf{S}}$, i.e., the operators $v_{i}$ are unitary
$v_{i}v_{i}^{\dagger}=v_{i}^{\dagger}{v}_{i}=\mathds{1}_{\mathsf{S}}$.

In the case of so-called UDD, the goal is to apply the pulses encompassing the
decoupling group can be applied between free evolution of the whole
$\mathsf{S}\mathsf{E}$ (i.e., between the application of the random RB gates)
in order to make use of the decoupling property in Eq. (30). This can further
be done in a more elaborated way (e.g., periodically or in a concatenated
fashion [48]) over some given time interval. While many schemes of DD are
known (not necessarily universal) [34], here we only employ UDD. Similarly,
while the dimension of system $\mathsf{S}$, can in principle be kept general
(as in RB), practical considerations, such as noise and the implementation
itself, imply that one normally considers just a single-qubit (for which we
can set $\mathds{V}=\\{\mathds{1},X,Y,Z\\}$, the single-Pauli group) or quite
small-dimensional systems.

We now consider noise on the full $\mathsf{S}\mathsf{E}$ which itself might
experience both global and local dissipation to another external environment;
this model covers a broad class of dynamical processes that may be used to
model quantum noise with both a coherent (unitary) and a stochastic
(dissipative) part. A similar argument to the one we make below can be seen in
[84, 36]. We can describe this noise model, on a given timescale, as an open
system evolution on $\mathsf{S}\mathsf{E}$ as

$\Lambda^{(t)}(\cdot)=\mathrm{e}^{t\mathcal{L}}(\cdot),$ (32)

where

$\displaystyle\mathcal{L}(\rho)$
$\displaystyle:=-i[H,\rho]+\mathcal{D}(\rho),$ (33) with
$\displaystyle\quad\mathcal{D}:=\mathcal{D}_{\mathsf{E}\mathsf{S}}+\mathcal{D}_{\mathsf{E}}+\mathcal{D}_{\mathsf{S}},\quad\text{where}\quad\mathcal{D}_{\mathsf{A}}(\cdot):=\sum_{k=1}^{d_{\mathsf{A}}^{2}-1}\gamma^{(\mathsf{A})}_{k}\left(L^{(\mathsf{A})}_{k}\rho{L}^{(\mathsf{A})\dagger}_{k}-\frac{1}{2}\\{L_{k}^{(\mathsf{A})\dagger}{L}^{(\mathsf{A})}_{k},\rho\\}\right),$
(34)

for a given $\mathsf{S}\mathsf{E}$ noise total Hamiltonian $H$ and so-called
dissipation maps $\mathcal{D}_{\mathsf{A}}$, which we take in a form such that
dissipation operators $L^{(\mathsf{A})}_{k}$ are traceless and orthonormal
($\operatorname{tr}[L^{(\mathsf{A})}_{k}]=0$ and
$\operatorname{tr}\leavevmode\nobreak\
[L^{(\mathsf{A})}_{i}L^{(\mathsf{A})\dagger}_{j}]=\delta_{ij}$), and with
positive constants $\gamma^{(\mathsf{A})}_{k}$. For notation, here $[a,b]:=ab-
ba$ and $\\{a,b\\}:=ab+ba$, and we denote by
$\mathcal{D}_{\mathsf{E}\mathsf{S}}$ the global $\mathsf{S}\mathsf{E}$
dissipator term.

For small time intervals $\delta{t}:=\tau$, this implies that to the first
order in $\tau$ we have

$\displaystyle\Lambda^{(\tau)}(\rho)\approx\rho-i\tau\,[H,\rho]+\tau\,\mathcal{D}(\rho).$
(35)

Now we consider an ideal UDD sequence from $\mathds{V}$ applied with the
decoupling operators acting between small windows of free noise evolution
$\tau$. For now, suppose that the operators $v_{i}\in\mathds{V}$ themselves
are noiseless (see Remarks in Sec. C.1). Define
$\Lambda_{v_{j}}^{(\tau)}:=\mathcal{V}_{j}\circ\Lambda^{(\tau)}\circ\mathcal{V}^{\dagger}_{j}$,
then,

$\Lambda_{v_{j}}^{(\tau)}(\rho)\approx\rho-i\tau\left[\mathcal{V}_{j}(H),\,\rho\right]+\tau\,\mathcal{V}_{j}\circ\mathcal{D}\circ\mathcal{V}_{j}^{\dagger}(\rho),$
(36)

so that we may compose all the “$v$-dressed” maps,
$\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{i=1}^{|\mathds{V}|}\Lambda_{v_{i}}^{(\tau)}(\rho)$,
and approximate this as a single evolution under an average Lindbladian
$\overline{\mathcal{L}}:=\sum_{i=1}^{|\mathds{V}|}\mathcal{V}_{i}\circ\mathcal{L}$
(again, with the assumption of a small enough $\tau$), as

$\displaystyle\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.7}{$\bigcirc$}}}}_{i=1}^{|\mathds{V}|}\Lambda_{v_{i}}^{(\tau)}\approx\mathrm{e}^{\tau\overline{\mathcal{L}}},$
(37)

whose action on $\rho$ is, in turn, to the first order in $\tau$,

$\mathrm{e}^{\tau\overline{\mathcal{L}}}(\rho)\approx\rho-i\tau\left[\sum_{j}\mathcal{V}_{j}(H),\,\rho\right]+\tau\sum_{j}\mathcal{V}_{j}\circ\mathcal{D}\circ\mathcal{V}^{\dagger}_{j}(\rho).$
(38)

The decoupling property in Eq. (30) implies that within this infinitesimally
small time window $\tau$, the second term on the right-hand-side in Eq. (38),
containing the Hamiltonian, will be uncorrelated. Regarding the remaining
terms, this does not occur in general for the global $\mathsf{S}\mathsf{E}$
dissipation term. This is problematic, analytically, because the quality maps,
and hence the ASF too, will remain as general non-Markovian objects unless we
assume there is no global $\mathsf{S}\mathsf{E}$ dissipation contribution in
the noise or it is otherwise negligible for small times.

Considering thus, only local dissipation, i.e.,
$\mathcal{D}\approx\mathcal{D}_{\mathsf{E}}+\mathcal{D}_{\mathsf{S}}$, in this
(infinitesimally) small time-interval $\tau$ limit, writing
$\tilde{L}^{(\mathsf{S})}_{j}:=v_{j}{L}_{k}^{(\mathsf{S})}v^{\dagger}_{j}$, we
end up with an evolution, on an initially uncorrelated state, of the form

$\displaystyle\mathrm{e}^{\tau\overline{\mathcal{L}}}(\rho)$
$\displaystyle\approx\rho+-i\tau|\mathds{V}|\,[B,\rho_{\mathsf{E}}]\otimes\rho_{\mathsf{S}}+\tau\,\mathcal{D}_{\mathsf{E}}(\rho_{\mathsf{E}})\otimes\rho_{\mathsf{S}}+\tau\,\rho_{\mathsf{E}}\otimes\mathcal{D}_{\mathsf{S}}^{(\mathsf{dd})}(\rho_{\mathsf{S}}),$
(39)

where here $\rho=\rho_{\mathsf{E}}\otimes\rho_{\mathsf{S}}$ and we defined
$\mathcal{D}_{\mathsf{S}}^{(\mathsf{dd})}$ by

$\displaystyle\mathcal{D}_{\mathsf{S}}^{(\mathsf{dd})}(\rho_{\mathsf{S}}):=\sum_{j}\sum_{k}\gamma^{(\mathsf{S})}_{k}\bigg{(}\tilde{L}^{(\mathsf{S})}_{j,k}\rho_{\mathsf{S}}\tilde{L}^{(\mathsf{S})\dagger}_{j,k}-\frac{1}{2}\left\\{\tilde{L}^{(\mathsf{S})\dagger}_{j,k}\tilde{L}^{(\mathsf{S})}_{j,k},\rho_{\mathsf{S}}\right\\}\bigg{)}.$
(40)

Now for the quality map in RB, suppose we have $\mathsf{S}\mathsf{E}$ noise of
the form $\Lambda_{k}=\mathrm{e}^{\tau_{k}\overline{\mathcal{L}}_{k}}$ at
time-step $k$, for given small time-intervals $\tau_{k}$ and DD-averaged
Lindbladians $\overline{\mathcal{L}}_{k}$ such that $\Lambda_{k}$ acts as in
Eq. (39) with corresponding operators labeled by subscripts $k$. Let us now
label the three summands in Eq. (39) as

$\displaystyle\Theta_{k}^{(0)}\otimes\Phi_{k}^{(0)}$
$\displaystyle:=\mathcal{I},$ (41)
$\displaystyle\Theta_{k}^{(1)}\otimes\Phi_{k}^{(1)}(\rho_{\mathsf{E}}\otimes\rho_{\mathsf{S}})$
$\displaystyle:=-i[B_{k},\rho_{\mathsf{E}}]\otimes\tau_{k}|\mathds{V}|\rho_{\mathsf{S}}+\mathcal{D}_{k,\mathsf{E}}(\rho_{\mathsf{E}})\otimes\tau_{k}\rho_{\mathsf{S}},$
(42)
$\displaystyle\Theta_{k}^{(2)}\otimes\Phi_{k}^{(2)}(\rho_{\mathsf{E}}\otimes\rho_{\mathsf{S}})$
$\displaystyle:=\rho_{\mathsf{E}}\otimes\tau_{k}\,\mathcal{D}_{k,\mathsf{S}}^{(\mathsf{dd})}(\rho_{\mathsf{S}}),$
(43)

so that
$\Lambda_{k}\approx\sum_{j_{k}=0}^{2}\Theta_{k}^{(j_{k})}\otimes\Phi_{k}^{(j_{k})}$,
and thus we obtain the quality map

$\displaystyle\hat{\mathcal{Q}}_{m,\pi}$
$\displaystyle=\sum_{j_{1},\ldots,j_{m}=0}^{2}\frac{\displaystyle{\prod_{i=m-1}^{0}}\operatorname{tr}\left(\hat{\Phi}_{i+1}^{(j_{i+1})}\hat{\mathcal{P}}_{\pi}\right)\,\hat{\Theta}_{i+1}^{(j_{i+1})}}{\operatorname{tr}\left(\hat{\mathcal{P}}_{\pi}\right)^{m}}+\mathcal{O}(\tau_{k}^{2}),$
(44)

where the first term can be of order up to $\tau_{1}\tau_{2}\cdots\tau_{m}$,
with individual $\mathcal{O}(\tau_{k}^{2})$ terms, for all $k=1,2,\ldots,m$,
neglected. Upon entering this noise map in the ASF in Eq. (22), we obtain the
term

$\displaystyle\hat{\operatorname{tr}}\,\hat{\Theta}_{m+1}^{(j_{m+1})}\cdots\hat{\Theta}_{0}^{(j_{0})}|\rho_{\mathsf{E}}\operatorname*{\\!\rangle\\!\rangle}=\operatorname{tr}\left[\Theta_{m+1}^{(j_{m+1})}\circ\cdots\circ\Theta_{0}^{(j_{0})}(\rho_{\mathsf{E}})\right]:=\begin{cases}0,\quad\text{if
any}\quad j_{i}=1\\\ 1,\quad\text{otherwise}\end{cases},$ (45)

so there are no contributions from $\mathsf{E}$ at order $\tau_{k}^{1}$, and
quality factors will only take the identity and the DD-dressed dissipators
$\mathcal{D}_{k,\mathsf{S}}^{(\mathsf{dd})}$. Thus to
$\mathcal{O}(\tau_{k}^{1})$ we have a purely Markovian decay of the form

$\displaystyle\mathcal{F}_{m}$
$\displaystyle=\sum_{\pi\in{R}_{\mathds{G}}}\sum_{j_{i}\in\left\\{0,2\right\\}}\left(\frac{\operatorname{tr}\hat{\Phi}^{(j_{1})}_{1}\hat{\mathcal{P}}_{\pi}}{\operatorname{tr}\hat{\mathcal{P}}_{\pi}}\right)\cdots\left(\frac{\operatorname{tr}\hat{\Phi}^{(j_{m})}_{m}\hat{\mathcal{P}}_{\pi}}{\operatorname{tr}\hat{\mathcal{P}}_{\pi}}\right)\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Phi}^{(j_{m+1})}_{m+1}\hat{\mathcal{P}}_{\pi}\hat{\Phi}^{(j_{0})}_{0}|\rho_{\mathsf{S}}\operatorname*{\\!\rangle\\!\rangle}+\mathcal{O}(\tau_{k}^{2}),$
(46)

which can be recast as

$\displaystyle\mathcal{F}_{m}$
$\displaystyle=\sum_{\pi\in{R}_{\mathds{G}}}\left(\frac{\operatorname{tr}\hat{\Omega}_{1}\hat{\mathcal{P}}_{\pi}}{\operatorname{tr}\hat{\mathcal{P}}_{\pi}}\right)\cdots\left(\frac{\operatorname{tr}\hat{\Omega}_{m}\hat{\mathcal{P}}_{\pi}}{\operatorname{tr}\hat{\mathcal{P}}_{\pi}}\right)\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Omega}_{m+1}\hat{\mathcal{P}}_{\pi}\hat{\Omega}_{0}|\rho_{\mathsf{S}}\operatorname*{\\!\rangle\\!\rangle}+\mathcal{O}(\tau_{k}^{2}),$
(47)

where
$\Omega_{k}:=\mathcal{I}+\tau_{k}\,\mathcal{D}^{(\mathsf{dd})}_{k,\mathsf{S}}$,
and where contributions with terms $\mathcal{O}(\tau_{k}^{2})$ for any
$k=0,\ldots,m+1$, are generally non-Markovian.

For concreteness, take the multi-qubit Clifford group, which is such that
111111 Equivalently a unitary 2-design; the projectors are
$\hat{\mathcal{P}}_{1}=\Psi$ and $\hat{\mathcal{P}}_{2}=\mathds{1}-\Psi$ where
$\Psi:=\sum|ii\rangle\\!\langle{jj}|/d_{\mathsf{S}}$. See e.g. [64]. For
details on the particular case of the ASF in the Clifford case, see [25].
$\mathcal{F}_{m}=Ap_{1}p_{2}\cdots{p}_{m}+B+\mathcal{O}(\tau_{k}^{2})$, with

$\displaystyle p_{k}$
$\displaystyle=\frac{\operatorname{tr}(\hat{\Omega}_{k})-1}{d_{\mathsf{S}}^{2}-1}=1+\tau_{k}\frac{\operatorname{tr}\left(\hat{\mathcal{D}}_{k,\mathsf{S}}^{(\mathsf{dd})}\right)}{d_{\mathsf{S}}^{2}-1}=1-\tau_{k}\frac{d_{\mathsf{S}}|\mathds{V}|\sum\gamma_{k}\operatorname{tr}\left(L^{(\mathsf{S})}_{k}L_{k}^{(\mathsf{S})\dagger}\right)}{d_{\mathsf{S}}^{2}-1}=1-\tau_{k}\frac{d_{\mathsf{S}}|\mathds{V}|\sum_{k_{\mathsf{S}}}\gamma_{k_{\mathsf{S}}}}{d_{\mathsf{S}}^{2}-1},$
(48) $\displaystyle A$
$\displaystyle=\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Omega}_{m+1}\hat{\Omega}_{0}|\rho_{\mathsf{S}}-\frac{\mathds{1}}{d_{\mathsf{S}}}\operatorname*{\\!\rangle\\!\rangle},\qquad
B=\operatorname*{\langle\\!\langle\\!}{M}|\hat{\Omega}_{m+1}\hat{\Omega}_{0}|\frac{\mathds{1}}{d_{\mathsf{S}}}\operatorname*{\\!\rangle\\!\rangle}\,$
(49)

where $A$, $B$ are standard Markovian SPAM noise contributions.

Finally, in the time-stationary case, $\Omega_{0}=\cdots=\Omega_{m+1}=\Omega$,
the ASF would further reduce to an exponential plus higher order terms,
$\mathcal{F}_{m}=Ap^{m}+B+\mathcal{O}(\tau^{2})$, with the single factor $p$
corresponding to the quality factor of the map
$\Omega=\mathcal{I}+\tau\,\mathcal{D}_{\mathsf{S}}^{(\mathsf{dd})}$. This is
presented in the main text as
$\mathcal{F}_{m}=Ap_{\tau_{\mathsf{dd}}}^{m}+B+\mathcal{O}(\tau_{\mathsf{dd}}^{2})$,
where

$\displaystyle p_{\tau_{\mathsf{dd}}}$
$\displaystyle=1-\tau_{k}\frac{d_{\mathsf{S}}|\mathds{V}|\sum_{k_{\mathsf{S}}}\gamma_{k_{\mathsf{S}}}}{d_{\mathsf{S}}^{2}-1},$
(50) $\displaystyle A$
$\displaystyle=\operatorname{tr}\left[M\Omega^{2}(\rho_{\mathsf{S}}-\mathds{1}/d_{\mathsf{S}})\right]\approx\operatorname{tr}\left[M(\rho_{\mathsf{S}}-\mathds{1}/d_{\mathsf{S}})\right]+2\tau_{\mathsf{dd}}\operatorname{tr}\left[M\mathcal{D}_{\mathsf{S}}^{(\tau_{\mathsf{dd}})}(\rho_{\mathsf{S}}-\mathds{1}/d_{\mathsf{S}})\right],$
(51) $\displaystyle B$
$\displaystyle=\operatorname{tr}\left[M\Omega^{2}(\rho_{\mathsf{S}})\right]\approx\operatorname{tr}\left[M(\rho_{\mathsf{S}})\right]+2\tau_{\mathsf{dd}}\operatorname{tr}\left[M\mathcal{D}_{\mathsf{S}}^{(\tau_{\mathsf{dd}})}(\rho_{\mathsf{S}})\right],$
(52)

with first-order $\tau_{\mathsf{dd}}$ approximations in $A$, $B$. This means
that when such approximation holds and the underlying noise has the form
described above, the data of a RB experiment, originally displaying non-
Markovian deviations, will fall to an exponential plus higher-order time
corrections if the ideal UDD pulses are applied between short time-windows
$\tau$ (or $\tau_{\mathsf{dd}}$); the average gate-fidelity extracted in such
a protocol would then contain that of $\Omega$ with respect to the identity
$\mathcal{I}$ to the first time-order.

###  C.1 Remarks

####  C.1.1 Noise model

The main reason to model the noise maps $\Lambda_{k}$ as a dynamical map is
twofold: one is to set a time scale between the action of the control gates,
parametrized by the time-intervals $\tau_{k}$, and the other is to express
both the generators of non-Markovianity (the unitary part) and Markovianity
(the dissipating part) within the global $\mathsf{S}\mathsf{E}$ noise
dynamics. This can encompass a broad class of physically motivated noise
models.

The only assumption made in order to obtain Eq. (47) was that either there is
no global $\mathsf{S}\mathsf{E}$ dissipation contribution in the noise
dynamics, or that it is negligible to the first order in $\tau$, as discussed
below Eq. (38). Notice that if there is no dissipation at all, local or
otherwise, to the first order in $\tau$, the ASF would be trivial (i.e., look
noiseless) and all higher-order contributions to the decay are generally non-
Markovian or from the local Hamiltonian. Numerically, in such case while it
may still be seen that the ASF approximates a decaying exponential to a large
extent –depending mainly on how fast the DD sequences are applied–, only
higher-orders would describe this decay.

Finally, the presence or absence of dissipator terms has to be explicitly
stated, as it is really the full $\Lambda^{(\tau)}$ map that generates noise.
For example, one could think of phrasing Result 3 in the main text for a
dilated Hamiltonian $H^{\prime}$ acting in a larger Hilbert space
$\mathsf{S}\mathsf{E}\mathsf{E}^{\prime}$ so that the Lindbladian
$\mathcal{L}^{\prime}(\cdot)=-i[H^{\prime},\cdot]$ contains no explicit
dissipators whatsoever; the point is that the actual original
$\Lambda^{(\tau)}$ will have to be written in terms of a partial tracing of
$\mathsf{E}^{\prime}$ and some fiducial initial state
$\rho_{\mathsf{E}^{\prime}}$, i.e.,
$\Lambda^{(\tau)}(\cdot)=\operatorname{tr}_{\mathsf{E}^{\prime}}\left[\mathrm{e}^{\tau\mathcal{L}^{\prime}}(\rho_{\mathsf{E}^{\prime}}\otimes\cdot)\right]$,
so that the dynamics could in general still contain global
$\mathsf{S}\mathsf{E}$ dissipators, albeit hidden implicitly through
$H^{\prime}$. In other words, the presence or absence of dissipation in the
dynamics as stated for Result 3 is really a statement of non-unitarity
(incoherence) or unitarity (full coherence) on $\mathsf{S}\mathsf{E}$,
respectively, of the noise map $\Lambda^{(\tau)}$.

####  C.1.2 Non-ideal dynamical decoupling

The assumption of instantaneous, infinitely strong DD pulses has been studied
extensively [48, 34], with the main consequence being that decoupling is not
achieved even at the first small-time-order; however, as long as the width of
the pulses is sufficiently small (relative to the inverse magnitude of the
largest singular value of the evolution Hamiltonian), decoupling is
approximately achieved. So in a sense, finite-width DD would still be
effective in Markovianizing RB, as long as pulse widths are sufficiently
small.

Similarly, one can relax the assumption of noiseless DD operators in a
straightforward way: as long as the noise is approximately local, time-
stationary and gate-independent, at each time-step of the RB protocol, this
would amount to a relabeling of the noise maps
$\Lambda_{i}\to\Lambda_{i}^{\prime}$, and Eq. (48) would hold at least up to
the second equality with a different map
$\mathcal{D}_{\mathsf{S}}^{\prime(\mathsf{dd})}$ accounting for noise in the
DD.

####  C.1.3 Single-qubit subsystem

For a single-qubit system $\mathsf{S}$, the Pauli group
$\mathds{V}=\\{\mathds{1},X,Y,Z\\}$ can be employed with either the so-called
$XY4$ or the $XZ4$ sequences. Letting $\mathcal{X}(\cdot):=X(\cdot)X$, and
similarly for the remaining Pauli operators, the ideal $XY(Z)4$ sequence takes
the form

$\mathcal{Y}(\mathcal{Z})\circ\Lambda^{(\tau)}\circ\mathcal{X}\circ\Lambda^{(\tau)}\circ\mathcal{Y}(\mathcal{Z})\circ\Lambda^{(\tau)}\circ\mathcal{X}\circ\Lambda^{(\tau)},$
(53)

or equivalently,
$\Lambda^{(\tau/2)}\mathcal{Y}(\mathcal{Z})\Lambda^{(\tau)}\mathcal{X}\Lambda^{(\tau)}\mathcal{Y}(\mathcal{Z})\Lambda^{(\tau)}\mathcal{X}\Lambda^{(\tau/2)}$,
where we removed the composition symbols for simplicity. These are arguably
the simplest UDD sequences; the reason why they implement the property in Eq.
(30) is because, up to a phase, they are both equivalent to interleaving all
of the single-qubit Pauli operators among the noise, which amounts to
averaging its generators over the subsystem qubit to a first-time order, as
shown in Appendix C. Other variants of the $XY(Z)4$ sequences achieving
higher-order decoupling exist, such as Concatenated DD [48], but we will not
employ them here.

####  C.1.4 Numerical example and the decoupling time $\tau_{\mathsf{dd}}$

The numerical example we consider in the main text is a
1$\mathsf{S}$+1$\mathsf{E}$ qubit model with noise dynamics of the form
$\Lambda^{(t)}=\mathrm{e}^{t\mathcal{L}}$, where
$\mathcal{L}(\cdot):=-i[H,\cdot]+\mathcal{D}_{\mathsf{S}}(\cdot)$ for
$H=JXX+h_{x}(XI+IX)+h_{y}(YI+IY)$, and we choose the parameters
$\tau_{\mathsf{rb}}=0.03$, $\rho_{\mathsf{S}}=M=|0\rangle\\!\langle 0|$,
$J=1.7$, $h_{x}=1.47$, $h_{y}=-1.05$ and Lindblad terms on the
$\mathsf{S}$-dissipator $L_{0}=X$, $\gamma_{0}=0.002$, $L_{1}=Z$ and
$\gamma_{1}=0.007$. We approximated the action of the exponential in the
Lindblad evolution up to $\mathcal{O}(t^{10})$.

In particular, here we discuss the statement that _how well DD Markovianizes
–whether it outputs decays closed to the purely exponential
$\mathcal{O}(\tau_{\mathsf{dd}})$ decay–, depends mainly on whether the time-
scale $\tau_{\mathsf{dd}}$ is small relative to that of the noise,
$\tau_{\mathsf{rb}}$_. This can be assessed by the separation between the
full-time-order ASF incorporating DD and the corresponding first-time-order.
This is easily tractable in our numerical example by inspecting this
separation of the outputs in Fig. 5 and/or other samples. This is done in Fig.
9, where the convergence of the ASF to the first time-order is clear, with the
separation in sequence length being linear and suppressed in smaller
$\tau_{\mathsf{dd}}$.

Figure 9: Behavior of the Markovianized full ASF with respect to the first
time-order with $XY4$ in $\tau_{\mathsf{dd}}$: We consider a relative
difference, $\Delta\mathcal{F}_{m}$, of the full numerical ASF (i.e., to all
time-orders) with single interleaved $XY4$ sequences with pulses on time-
intervals $\tau_{\mathsf{dd}}$, without uncertainties, with respect to the
analytical exponential first time-order, as given by Eq. (6); in (a) we plot
the values of $\Delta\mathcal{F}_{m}$ with respect to $\tau_{\mathsf{dd}}$
(with $\tau_{\mathsf{dd}}\approx\tau_{\mathsf{rb}}$ and
$\tau_{\mathsf{dd}}\ll\tau_{\mathsf{rb}}$) and sequence length $m$, showing it
generally decreases together with $\tau_{\mathsf{dd}}$ and increases in $m$,
while in (b) it is shown that $\Delta\mathcal{F}_{m}$ increases linearly in
most cases by fitting a linear function with a rate proportionally small in
$\tau_{\mathsf{dd}}$.

The behavior in larger, $\tau_{\mathsf{dd}}>2\tau_{\mathsf{rb}}$ can quickly
become unreliable as uncertainties increase and the ASF drops quickly as the
sequence length increases, so generally, in practice the timescale of the RB
sequences, $\tau_{\mathsf{rb}}$, can serve as the reference timescale with
which $\tau_{\mathsf{dd}}$ should be minimized.

##  Appendix D Subsystem Pauli-twirled noise

Other Markovianizing strategies different from DD could also be conceivable.
In particular, techniques tailoring noise into a so-called _Pauli channel_
have been shown to reduce the impact of non-Markovian noise [39, 35] in a
statistically significant way.

Let us consider for now a closed $n$-qubit system. Let $\mathds{P}_{n}$ be the
$n$-qubit Pauli group, i.e., that generated by $n$-fold Kronecker products of
single-qubit Pauli operators. The so-called twirl,
$\mathcal{T}_{\mathds{P}_{n}}$, over $\mathds{P}_{n}$, of a CP map between
$n$-qubits, $\Phi(\cdot):=\sum_{\mu}\phi_{\mu}(\cdot)\phi_{\mu}^{\dagger}$,
can be written as

$\displaystyle\mathcal{T}_{\mathds{P}_{n}}\left(\hat{\Phi}\right)$
$\displaystyle:=\frac{1}{4^{n}}\sum_{P\in\mathds{P}_{n}}P^{\dagger}\phi_{\mu}{P}\otimes{P}^{\mathrm{T}}\phi_{\mu}^{*}{P}^{*},$
(54)

where $\hat{\Phi}:=\sum_{\mu}\phi_{\mu}\otimes\phi_{\mu}^{*}$ is the
vectorized form of $\Phi$, as defined in Appendix A.1. Hence, expanding the
Kraus operators in the Pauli basis,
$\phi_{\mu}=\sum_{i}\alpha_{i}^{\mu}{P}_{i}$, we can rewrite the Pauli-twirled
channel as

$\displaystyle\hat{\Phi}^{\mathsf{P}}:=\mathcal{T}_{\mathds{P}_{n}}\left(\hat{\Phi}\right)$
$\displaystyle=\frac{1}{4^{n}}\sum\alpha_{i}^{\mu}\alpha_{j}^{*\mu}P_{k}^{\dagger}{P}_{i}{P}_{k}\otimes{P}_{k}^{\mathrm{T}}P_{j}^{*}{P}_{k}^{*}=\sum{p}_{i}P_{i}\otimes{P}_{i}^{*},$
(55)

with $p_{i}:=\sum_{\mu}|\alpha_{i}^{\mu}|^{2}$. In general, any quantum
channel of the form of Eq. (55) is called an _Pauli channel_. Such Pauli
channel $\Phi^{\mathsf{P}}$ has Kraus operators $\sqrt{p_{i}}P_{i}$, where
$p_{i}=\sum_{\mu}|\alpha_{i}^{\mu}|^{2}$ define a probability distribution
satisfying $p_{i}\geq 0$ and $\sum_{i}p_{i}\leq 1$ (saturated if the channel
$\hat{\Phi}$ is also TP). Colloquially, then, we can say that Pauli-twirling a
quantum channel, has the effect of tailoring it into a Pauli channel. In the
specific case where we think of quantum channels as modeling quantum noise, we
may refer to $p_{i}$ as the _Pauli error-rates_ , i.e., the probability
distribution of any Pauli error from happening.

Throughout, we will reserve the zeroth Pauli to denote the identity,
$P_{0}:=\mathds{1}$.

A particularly useful representation when dealing with multi-qubit systems is
the so-called Pauli Transfer Matrix (PTM) representation; in such
representation, Pauli-twirled channels take a manifestly diagonal form.
Explicitly, the PTM representation of an $n$-qubit quantum channel $\Phi$, is
also a vectorized form (i.e., matrix-form) of the channel, but expanded in the
Pauli basis,

$R_{\Phi}:=2^{-n}\sum_{Q,P\in\mathds{P}_{n}}\operatorname*{\langle\\!\langle\\!}{Q}|\Phi(P)\operatorname*{\\!\rangle\\!\rangle}\,|Q\operatorname*{\\!\rangle\\!\rangle}\\!\operatorname*{\langle\\!\langle\\!}{P}|,$
(56)

that is, the PTM of $\Phi$ is defined by the elements

$(R_{\Phi})_{ij}=2^{-n}\operatorname*{\langle\\!\langle\\!}{P_{i}}|\Phi(P_{j})\operatorname*{\\!\rangle\\!\rangle}=2^{-n}\operatorname{tr}\left[P_{i}\,\Phi(P_{j})\right],\qquad{P}_{i},P_{j}\in\mathds{P}_{n}.$
(57)

In the case of a Pauli channel, $\Phi^{\mathsf{P}}$, the corresponding PTM
matrix $R_{\Phi^{\mathsf{P}}}$ is diagonal, with its (real) elements being
related to the probabilities $p_{i}$ in Eq. (55) by a so-called Walsh-Hadamard
transform (see e.g. [12]).

Remarkably, a Pauli-twirled noise channel has the same average gate-fidelity
as the bare noise channel; this can be seen directly in Eq. (20), as the
average gate-fidelity only takes into account the diagonal error rates. This,
in turn, implies that any technique estimating average gate-fidelities will
remain reliable under Pauli noise, and not greatly overestimate gate-quality
by missing out on all other non-diagonal terms: indeed, the diamond norm for
quantum channels with a given fidelity is minimized by Pauli channels [86, 87,
54] and can be explicitly computed in such case through the trace-norm of the
vector of Pauli error rates [4]. Pauli-twirling can be operationally
implemented through protocols such as Pauli Frame Randomization [88, 39] or
Randomized Compiling [53, 35].

We can thus now consider our more general open-system scenario where our
subsystem $\mathsf{S}$ is a $n_{s}$-qubit system, with the limiting case of
perfect twirling of some noise channel $\Lambda$ with the Pauli group solely
on subsystem $\mathsf{S}$. We refer to this as the $\mathsf{S}$-Pauli-twirled
channel, $\Lambda^{\mathsf{P}}$, associated to $\Lambda$. For simplicity, let
us assume $\mathsf{E}$ is also multi-qubit; otherwise Pauli operators should
be replaced with a corresponding $\mathsf{E}$ operator basis. Let
$\hat{\Lambda}=\sum_{\nu}\lambda_{\nu}\otimes\lambda_{\nu}^{*}$, and expand
the Kraus operators $\lambda_{\nu}$ in the Pauli basis as
$\lambda_{\nu}:=\sum_{\mu,i}\alpha_{\mu,i}^{\nu}{P}_{\mu}\otimes{P}_{i}$.
Define
$\sum_{\nu}\alpha_{\mu,i}^{\nu}\alpha_{\sigma,j}^{*\nu}:=\chi_{\mu\sigma,ij}$.
Then we can write 121212 Strictly, the ordering in the vectorized (hat
$\hat{\mathcal{X}}$) representation, as we defined it, is
$\mathsf{E}\mathsf{S}\otimes\mathsf{E}\mathsf{S}$; we are slightly abusing
notation for clarity of presentation.

$\displaystyle\hat{\Lambda}^{\mathsf{P}_{\mathsf{S}}}:=\left(\mathcal{I}_{\mathsf{E}}\otimes\mathcal{T}_{\mathds{P}_{n_{\mathsf{S}}}}\right)\hat{\Lambda}$
$\displaystyle=\sum\chi_{\mu\sigma,ij}\left(P_{\mu}\otimes{P}_{\sigma}^{*}\right)\otimes\mathcal{T}_{\mathds{P}_{n_{\mathsf{S}}}}\left(P_{i}\otimes{P}^{*}_{j}\right)$
$\displaystyle=\sum\chi_{\mu\sigma,ii}\left(P_{\mu}\otimes{P}_{\sigma}^{*}\right)\otimes\left(P_{i}\otimes{P}^{*}_{i}\right),$
(58)

which means that twirling solely over subsystem $\mathsf{S}$ generates a
Pauli-like channel on system $\mathsf{S}$ that retains the dependence on the
environment $\mathsf{E}$ only on the local Pauli error rates. The
corresponding PTM is not diagonal, despite nullifying many entries on the
$\mathsf{S}$ subspace, as exemplified in Fig. 10. If the environment is traced
out, then the reduced channel is a Pauli channel with $\mathsf{E}$-dependent
Pauli error rates,
$\operatorname{tr}_{\mathsf{E}}[\Lambda^{\mathsf{P}_{\mathsf{S}}}(\cdot)]\propto\sum{p}_{i}(\rho^{\prime}_{\mathsf{E}})P_{i}\otimes{P}_{i}$
for some reduced state of the environment $\rho^{\prime}_{\mathsf{E}}$; in
multi-time processes, however, non-Markovianity is still carried through this
environmental dependence.

Figure 10: Effect of twirling a subsystem with the Pauli group: (Left) PTM of
the model used in the main text, $\Lambda=\exp(\tau\mathcal{L})$ where
$\mathcal{L}(\cdot):=-i[H,\cdot]+\mathcal{D}_{\mathsf{S}}(\cdot)$ for fixed
$\tau=0.03$, $H=JXX+h_{x}(XI+IX)+h_{y}(YI+IY)$, and chosen parameters
$\rho_{\mathsf{S}}=M=|0\rangle\\!\langle 0|$, $J=1.7$, $h_{x}=1.47$,
$h_{y}=-1.05$ and Lindblad terms $L_{0}=X$, $\gamma_{0}=0.002$, $L_{1}=Z$,
$\gamma_{1}=0.007$; (Right) PTM of the corresponding $\mathsf{S}$-Pauli
twirled channel, $\Lambda^{\mathsf{P}_{\mathsf{S}}}$.

###  D.1 ASF under $\mathsf{S}$-Pauli-twirled noise

Consider the limiting case where Pauli twirling can be implemented effectively
perfectly, via Randomized Compiling or otherwise, alongside RB for non-
Markovian noise, so that quality maps $\mathcal{Q}_{\pi,m}$ are defined in
terms of maps $\hat{\Lambda}_{i}^{\mathsf{P}_{\mathsf{S}}}$ of the form of Eq.
(58).

We now denote time-step indices as superscripts in $\chi$ and as a labeling
for each sum index, i.e.,

$\hat{\Lambda}^{\mathsf{P}_{\mathsf{S}}}_{i}=\sum\chi_{\mu_{i}\sigma_{i},k_{i}k_{i}}^{(i)}\left(P_{\mu_{i}}\otimes{P}_{\sigma_{i}}^{*}\right)\otimes\left(P_{k_{i}}\otimes{P}_{k_{i}}^{*}\right),$
(59)

then, plugging a set of
$\\{\Lambda_{1}^{\mathsf{P}_{\mathsf{S}}},\ldots,\Lambda_{m}^{\mathsf{P}_{\mathsf{S}}}\\}$
noise maps into a quality map of length $m$ (defined in Eq. (23)), denoting by
$\\{|\epsilon_{i}\rangle,|\epsilon_{i}^{\prime}\\}_{i=0}^{m}$ basis vectors of
$\mathsf{E}$, we have

$\displaystyle\hat{\mathcal{Q}}_{m,\pi}$
$\displaystyle=\sum_{\\{\epsilon_{j},\epsilon_{j}^{\prime}\\}_{j=0}^{m}}^{d_{\mathsf{E}}}\frac{\prod_{i=m-1}^{0}\operatorname{tr}\left(\langle\epsilon_{i+1}\epsilon_{i+1}^{\prime}|\hat{\Lambda}^{\mathsf{P}_{\mathsf{S}}}_{i+1}|\epsilon_{i}\epsilon_{i}^{\prime}\rangle\hat{\mathcal{P}}_{\pi}\right)}{\operatorname{tr}\left(\hat{\mathcal{P}}_{\pi}\right)^{m}}|\epsilon_{m}\epsilon_{m}^{\prime}\rangle\\!\langle\epsilon_{0}\epsilon_{0}^{\prime}|$
$\displaystyle=\sum\frac{\prod_{i=m-1}^{0}\chi^{(i+1)}_{\mu_{i+1}\sigma_{i+1},{k}_{i+1}{k}_{i+1}}\langle\epsilon_{i+1}|P_{\mu_{i+1}}|\epsilon_{i}\rangle\\!\langle\epsilon_{i+1}^{\prime}|P^{*}_{\sigma_{i+1}}|\epsilon_{i}^{\prime}\rangle\operatorname{tr}\left[\left(P_{k_{i+1}}\otimes{P}^{*}_{k_{i+1}}\right)\hat{\mathcal{P}}_{\pi}\right]}{\operatorname{tr}\left(\hat{\mathcal{P}}_{\pi}\right)^{m}}|\epsilon_{m}\epsilon_{m}^{\prime}\rangle\\!\langle\epsilon_{0}\epsilon_{0}^{\prime}|$