# Impact of dynamics, entanglement, and Markovian noise on the fidelity of
few-qubit digital quantum simulation

Max D. Porter<EMAIL_ADDRESS>Fusion Energy Sciences Program, Lawrence
Livermore National Laboratory  Ilon Joseph  Fusion Energy Sciences Program,
Lawrence Livermore National Laboratory

(June 9, 2022)

###### Abstract

For quantum computations without error correction, the dynamics of a
simulation can strongly influence the overall fidelity decay rate as well as
the relative impact of different noise processes on the fidelity. Theoretical
models of Markovian noise that include incoherent Lindblad noise, gate-based
errors, and stochastic Hamiltonian noise qualitatively agree that greater
diffusion and entanglement throughout Hilbert space typically increase the
fidelity decay rate. Simulations of the gate-efficient quantum sawtooth map
support these predictions, and experiments performed on three qubits on the
IBM-Q quantum hardware platform at fixed gate count qualitatively confirm the
predictions. A pure depolarizing noise model, often used within randomized
benchmarking (RB) theory, cannot explain the observed effect, but gate-based
Lindblad models can provide an explanation. They can also estimate the
effective Lindblad coherence times during gates, and find a consistent
$2\text{--}3\times$ shorter effective $T_{2}$ dephasing time than reported for
idle qubits. Additionally, the observed error per CNOT gate exceeds IBM-Q’s
reported value from RB by $3.0\times$ during localized dynamics and
$4.5\times$ during diffusive dynamics. This demonstrates the magnitude by
which RB understates error in a complex quantum simulation due to dynamics and
crosstalk.

##  1 Introduction

###  1.1 Motivation

Quantum computers may eventually become indispensable to scientific computing
due their ability to accelerate many calculations of interest. While achieving
this goal will require error correction, the noisy intermediate scale quantum
(NISQ) devices that are available today provide a unique platform for
exploring the dynamics of many-body quantum systems. The ubiquitous presence
of interactions with the environment generates noise which adds complexity to
the interpretation of the results. Depending on context, noise can affect the
dynamics in different ways. It may completely wash out the dynamics of
interest or it can potentially be used to measure key signatures of the
dynamics [1]. In this work we ask the question: how does the dynamics
influence the effects of noise? We show that varying the dynamics of a chaotic
quantum map alters the fidelity decay due to noise in three different noise
models, and demonstrate that these effects are observable today using
experiments on the IBM-Q Quantum Lab.

Quantum simulation is one of the main applications of interest for achieving
near term quantum advantage [2]. Simulating chaotic dynamics is of particular
interest due to the provable difficulty of simulating chaos. Recent claims of
quantum supremacy relied on the difficulty of simulating chaotic quantum
circuits [3, 4] and among paths to quantum advantage simulating chaotic
dynamics may be the most qubit-efficient [5]. Measuring the decay of fidelity
can also provide an exponentially efficient measure of the classical Lyapunov
exponent [6, 7] and chaotic decoherence [8]. Measuring dynamical localization
of classically chaotic systems may yield a similar speedup [9]. The fidelities
of quantized versions of classically chaotic Hamiltonian systems can decay
differently (often faster) than for integrable dynamics [10, 11, 1]. Quantum
chaos can magnify the effect of Trotter errors [12] while quantum localization
can bound Trotter errors [13].

Both chaotic diffusion and dynamical localization can be observed in the
quantum sawtooth map, which also has a particularly efficient algorithm for
running on gate-based quantum devices. The transition from diffusion to
localization occurs when the ratio of diffusion strength to effective $\hbar$
is small and quantum interference dominates. Refs. [14] and [15] initiated the
use of the quantum sawtooth map to characterize experimental noise by using
the degree of localization as a test of device fidelity. We extend this by
studying the phase transition to diffusion and the Loschmidt echo fidelity
throughout this transition, opening the door to probing the interaction
between quantum map dynamics and experimental noise.

###  1.2 Summary of results

The main results are: (1) the dynamics being simulated affect the fidelity
decay rate (2) due to the amount of entanglement and superposition that is
generated. (3) Theoretical analysis quantifies how a randomly entangled state
combines the effects of relaxation and dephasing errors, while (4) a gate-
based Lindblad model can explain experimental observations and (5) allows one
to measure effective $T_{1}$ and $T_{2}$ times for quantum gates. These
effective coherence times more accurately describe the duration of complex
quantum simulations, and the gate-based Lindblad model easily scales to many
qubits.

The effect of the dynamics on the fidelity decay is demonstrated through
theory, simulation, and experiment. The results of different noise models are
compared via numerical fit in Sec. 4.2.2. The CNOT average error per gate for
different dynamics is extracted from experiment and reported in Table 2. The
effective $T_{1}$ and $T_{2}$ times are also extracted from the model fits and
reported in Table 3.

Our gate-based dynamical Lindblad model has two main components. First,
theoretical expressions are derived for the fidelity decay of a noisy
Hamiltonian simulation during localized dynamics, chaotic/diffusive dynamics,
and the middling case of unentangled superposition. (The term diffusive
dynamics will be preferred since the dynamics are not clearly chaotic until a
sufficient number of qubits is used [1].) To derive these expressions, each
type of dynamics is approximated as an appropriate steady state that gradually
decoheres under Lindblad evolution. The resulting expressions reveal differing
dependencies on the effective $T_{1}$ and $T_{2}$ times due to different
degrees of random entanglement and superposition. The second component of this
model is to connect it to the gate picture relevant to IBM-Q and many other
platforms. Rather than qubits always decaying at a constant rate, they are
assumed to decay more rapidly when acted on by a gate. On IBM-Q, the most
error-prone gates are the two-qubit CNOT gates and, hence, we focus our
attention on these gates alone. This retains the expected error dependence on
the number of CNOT gates.

In reality errors come in many varieties not described by a Lindblad model.
The most comprehensive description of Markovian gate errors requires gate set
tomography [16, 17] to find a consistent set of process matrices for all
gates. Yet the task of measuring all such processes is not scalable to many
qubits and is not possible with the limited runtime on the open-access IBM-Q
system. A reduced-order Lindlbad model can still capture some of the key
aspects of a process matrix better than the randomized benchmarking (RB) and
single-qubit metrics currently used by IBM-Q. When a depolarizing process,
e.g. determined by RB, is used to model the evolution [18], it provides errors
that are independent of the density matrix (except for the fully mixed state),
and hence the dynamics. Even the simplest Lindblad model captures more detail
than this.

We briefly consider a stochastic error model that also predicts dynamics-
dependent fidelity decay [19]. In this model a parameter of the Hamiltonian is
perturbed stochastically during each discrete step of the chaotic quantum
sawtooth map, though a static Hamiltonian perturbation model would have
achieved a similar effect [20]. Parameter noise models correspond to lower-
frequency noise that interacts more strongly with the dynamics, such as
coherent error that varies from gate to gate, while fixed perturbation
corresponds to zero frequency, namely static coherent error. Parameter noise
effects may be qualitatively similar to measured pink ($1/f$) noise dephasing
processes [21]. Both types of noise predict an algebraic fidelity decay rate
when the dynamics being simulated are localized, or more generally regular
[19, 22, 1]. This is attributed to the slower separation of non-chaotic
semiclassical path integral trajectories. However, through simulation we find
this effect requires more than three qubits and so is not yet experimentally
achievable with the quantum map and hardware used here. In the chaotic regime
the parameter noise model also allows one to measure the Lyapunov exponent,
though at least six qubits are required to observe it [6, 1]. We find that
these parameter noise effects are compatible with the presence of Lindblad
noise and might be observable in the near future.

###  1.3 Rationale for gate-based Lindblad model

Working with a complete decoherence model of quantum gates is quite
challenging because there are many terms that are potentially active [23]. For
a linear trace-preserving quantum process model, the number of terms is
$N^{4}-N^{2}$, where $N$ is the number of states. Clearly, for all-to-all
connectivity, a complete characterization is not scalable to many qubits.

Performing quantum process tomography [24] for a single entangling gate
already requires a relatively large number of measurements. While there are
software packages such as pyGSTi [16, 17] that allow one to automatically
characterize a complete set of two-qubit gates, many of the commercial quantum
hardware platforms available today do not allow enough run time to perform a
complete characterization of even a single entangling gate. Moreover, because
the complete protocol takes so long to run, the hardware characteristics may
drift significantly over the duration of the experiments.

To make headway under such circumstances, one must attempt to utilize
plausible simplifications of the number and types of decoherence processes
under consideration. Assume that for each gate, one can identify a simplified
set of decoherence processes that are dominant. One could then perform a
restricted form of process tomography to characterize each gate. One could
even characterize an entire gate set – assuming that the number of independent
decoherence processes is small enough that this procedure is scalable.

Presumably, a detailed understanding of the physical processes that mediate
decoherence mechanisms would yield a relatively concise mathematical
description of the dominant decoherence processes. However, this could still
potentially be quite complicated, for instance including non-Markovian
processes, and would require custom experiments to diagnose and verify.
Moreover, this could potentially require a physical model for the environment
as well as for the quantum computing hardware, making the model non-Markovian
and more complex than the process tomography paradigm typically admits.

Yet simply understanding the way that the physical system imposes a
mathematical structure on the decoherence processes is potentially quite
valuable. For example, theoretical studies have proven that the Lindblad
master equation [23] is the appropriate infinitesimal generator of the
temporal evolution of a quantum system coupled to an incoherent process that
is both trace-preserving and completely positive. Hence, this implies that the
Lindbladian form should be considered primary and that integrating the master
equation over a finite time interval should be used to determine the
associated process matrix. Clearly, a sparse representation in terms of
Lindblad operators will typically yield a rather dense process matrix. Thus,
from the physical perspective, it is much more likely for the Lindblad
representation to be sparse than it is for the associated process matrix.

Another plausible physical assumption is that the dominant noise processes are
from one- and two-qubit interactions. This reduces the number of processes
exponentially, from $O(N^{4})$ (the elements of the infinitesimal process
matrix) when including up to all-qubit processes to just $O(\log^{2}(N))$ (the
number of qubit pairs) for up to two-qubit processes. In the infinitesimal
process matrix this corresponds to an exponential reduction in the number of
elements. This can be seen by considering the relationship between the
Lindblad equation, which describes the time evolution of the density matrix,
and the infinitesimal process matrix, which describes the coupling of density
matrix elements. Most process matrix elements of a large system couple density
matrix elements between which more than two qubits have differing states. So
limiting to one- and two-qubit processes treats most process matrix elements
as zero, leaving an exponentially small number of elements active.

Two-qubit processes could be further limited to fewer physically plausible
channels. For instance photon exchange between qubits mediated by the
environment could be an important channel, since it requires only one photon
and is an extension of single-qubit relaxation and excitation processes [25].
A stronger restriction is limiting interactions to nearest neighbors in a
sparsely connected device, reducing the number of processes to $O(\log(N))$.
This also strengthens the assumption that no more than two-qubit processes
need to be considered.

Two-qubit dephasing is a physically plausible channel that we found in our
numerical studies to have very similar qualitative effects as the standard
single-qubit dephasing. Presumably, the reason is that each multi-qubit
dephasing process can be understood as single-qubit dephasing in an alternate
basis. Understanding whether or not such non-standard dephasing terms are
important for describing actual quantum hardware performance could be an
interesting subject for future work.

The greatest physical simplification that can plausibly be applied to a
Lindblad model is that only the single-qubit relaxation and dephasing
processes are important, averaged over all qubits for simplicity, but that
they are enhanced over their measured values during gate operation. (For
typical quantum hardware platforms, the temperature is so far below the
excitation energy that excitation processes can be considered subdominant.)
While this may be a vast oversimplification, these processes are clearly
important because they are universally required to model experimental data in
practice. The physical meaning of this model is simple: the physical processes
that actuate the gate renormalize the interaction between the qubit and the
environment, and, hence, renormalize the characteristic single qubit
relaxation and dephasing rates. Note, however, that the predictions are
nontrivial because the associated time-integrated process matrix has a
relatively dense form that depends on the products of the rates with the
overall gate time.

While many other decoherence processes are likely needed for a truly accurate
description of the quantum hardware evolution, we have found that this
simplest Lindblad model is sufficient for modeling the experimental data
described in this work. Our work here has the restricted goal of attempting to
determine the single-qubit decoherence enhancement for the CNOT gates that
clearly dominate the IBM-Q error budget. Moreover, given the limited
experimental data that one is able to obtain using today’s platforms, the data
may not contain enough information to robustly determine additional free
parameters over and above those contained in this simple model. In the future,
it could be interesting to explore a more comprehensive approach to the
determination of the relevant decoherence processes. Perhaps one could even
attempt a version of gate set tomography [16] that deduces the relevant
relaxation and dephasing enhancements for every qubit and every gate in the
set.

###  1.4 Overview of contents

In Sec. 2.1 and Sec. 2.2 the quantum sawtooth map is introduced and the
conditions for dynamical localization are described. In Sec. 2.3 the noise
models under consideration are described in detail. Simplified analytical
models are compared to numerical results and scaling with number of qubits is
explored. Section 2.3.2 and Sec. 2.3.3 consider the Lindblad noise model, Sec.
2.3.4 considers the qiskit Aer noise model, and Sec. 2.3.5 consides the
parameter noise model. Presented next is the implementation of the quantum
sawtooth map as an algorithm. Sec. 3.1 gives its gate decomposition while Sec.
3.2 discusses the details of native gates, device connectivity, and gate depth
scaling. In Sec. 4 experimental results are presented, compared to IBM-Q
reported metrics, and fits using various models are performed. Lastly Sec. 5
provides a summary of our results.

##  2 Theoretical analysis

###  2.1 Quantum sawtooth map

The quantum sawtooth map (QSM) is defined by the time-periodic Hamiltonian

$H_{QSM}=\frac{\hat{J}^{2}}{2}-\sum_{n}K\frac{\hat{\theta}^{2}}{2}\delta(t-n)\text{
for }\theta\bmod 2\pi$ (1)

which is derived from the classical Hamiltonian by quantizing to get
$\hat{J}=\hbar\hat{p}$ and discretizing to get $\hat{\theta}=2\pi\hat{q}/N$ to
yield momentum and position operators $\hat{p}$ and $\hat{q}$ respectively,
with eigenvalues $-N/2\leq p,q<N/2$. The quantum evolution propagator over one
period is then

$\displaystyle U_{QSM}$
$\displaystyle=\hat{\mathcal{T}}e^{-i\int_{0}^{1}H_{QSM}dt/\hbar}$
$\displaystyle=U_{kin}U_{pot}$ (2) $\displaystyle U_{pot}(\hat{q})$
$\displaystyle=e^{ik(\beta\hat{q})^{2}/2}$ $\displaystyle U_{kin}(\hat{p})$
$\displaystyle=e^{-i\hbar{\hat{p}}^{2}/2}$

where $\hat{\mathcal{T}}$ is the time-ordering operator, $k\equiv K/\hbar$ is
the quantum kicking parameter, and $\beta\equiv 2\pi/N$ for $N$ basis states.
In the context of quantum computing $N=2^{n}$ for $n$ qubits. This single-
period propagator is often called a Floquet operator, but since it corresponds
to a classical map it is also known as a quantum map. Periodicity matching
between the classical and quantum systems gives $\hbar=2\pi L/N$ for positive
integer $L$ [1]. Note that when N is a power of two, making $L$ even would
allow $L/N$ to simplify and reduce the pseudorandomness of the eigenphases of
$U_{pot}$ and $U_{kin}$, disrupting the chaotic dynamics, so $L$ is assumed
odd.

Note a partial symmetry between the phases of $U_{kin}$ and $U_{pot}$:
$\hbar=L*2\pi/N$ while $-k\beta^{2}=-K/L*2\pi/N$. Typically the qubits are
mapped to the $p$-basis, but mapping instead to the $q$-basis would swap the
roles of $L$ and $-K/L$.

###  2.2 Localization and initial conditions

When the classical diffusion coefficient $D_{K}$ is small compared to
$\hbar^{2}$, the QSM is predicted after the Heisenberg time to reach a steady
state that is exponentially localized around an initial momentum state
$\ket{p_{0}}$ as

$P_{p}=|\innerproduct{p}{\psi}|^{2}\approx\frac{1}{\ell}e^{-\frac{2|p-p_{0}|}{\ell}}$
(3)

with localization length [7]

$\ell\approx D_{K}/\hbar^{2}.$ (4)

The Heisenberg time or “break” time [7] is the time to resolve the energy
levels, defined as the inverse mean energy level spacing [26, 27].

The presence of broken cantori in the classical system can slow diffusion at
small $K$, so there are two regimes given by

$\displaystyle D_{K}\approx\begin{cases}(\pi^{2}/3)K^{2}\qquad\text{ for
}K>1\\\ 3.3K^{5/2}\quad\text{ for }0<K<1.\\\ \end{cases}$ (5)

Due to the periodicity and finite size of the system, localization only occurs
if the localization length is small enough to have a global maximum at its
central peak. This occurs when

$\displaystyle k$ $\displaystyle<k_{\text{loc}}$
$\displaystyle\equiv\begin{cases}\sqrt{\frac{3}{\ln(2)\pi^{2}}N}\approx
0.66N^{1/2}\quad\text{ for }K>1\\\
\left(\frac{1}{3.3\sqrt{2\pi}\ln(2)}\frac{N^{3/2}}{L^{1/2}}\right)^{2/5}\approx
0.50N^{3/5}L^{-1/5}\\\ \qquad\qquad\qquad\qquad\quad\text{ for
}0<K<1\end{cases}$ (6) $\displaystyle=\max(0.66N^{1/2},0.50N^{3/5}L^{-1/5})$

[1] with diffusion occuring otherwise. These two dynamical regimes are
demonstrated in Fig. 1.

Figure 1: Exact noiseless simulations of the quantum sawtooth map (QSM),
showing the localized case $k=0.1$ (blue) and the diffusive case $k=4.55$
(red) for $t=1,2,4,8$.

When localization is strong ($\ell\ll N$) the above formula for average $\ell$
does not uniformly apply to all initial conditions, with $\ell$ instead
depending on initial condition in a manner dependent on $L$. Less localized
groups of Floquet eigenstates appear at $L$ equally spaced locations in
momentum space, reducing $\ell$ for initial momentum states in the vicinity of
these groups. Once $L\sim N$ this transitions to more uniform $\ell$ with just
$\ket{p=0}$ strongly localized. This is how Refs. [14, 15] used $N=8,L=7$ to
obtain strong localization at $\ket{p=0}$, despite other states being less
localized. In this study we fix $L=1$ to keep the state-dependent effect
constant. When demonstrating strong localization in Fig. 1 we show a single
$\ket{p\neq 0}$, while fidelity such as in Fig. 3 is averaged over all initial
conditions to avoid state-dependent effects.

It is also worth noting that the initial conditions are chosen to be momentum
eigenstates rather than Floquet eigenstates. In the strongly localized limit,
Floquet eigenstates are complex superpositions of $\ket{p}$ and $\ket{-p}$. As
the quantum map causes the phase of each Floquet state to evolve at a rate
proportional to its quaisenergy, pairs of states will fall out of phase,
causing an initial state $\ket{p}$ to evolve to $\ket{-p}$ after a phase
difference $\pi$. However, in the localized limit the quasienergies of these
Floquet state pairs are close together, making this process long and
irrelevant for the short time scales discussed in this paper. The fact that we
reverse the map to calculate fidelity further reduces the relevance of this.

###  2.3 Noise models and fidelity

####  2.3.1 Types of noise and fidelities

Present day quantum devices are impacted by many different types of noise.
This motivates studies of the types of noise that are present and how the
errors impact algorithms of interest. On the IBM-Q platform errors occur
primarily during two-qubit gates which in the QSM algorithm contribute about
$10\times$ more to the total error over single-qubit gates, as discussed in
Sec. 3.2. The types of error that occur may be incoherent Markovian relaxation
and dephasing error ($T_{1}$ and $T_{2}$ respectively), coherent error, multi-
qubit incoherent errors, or something else. Here three models of noise in a
quantum map are considered for understanding the effects of errors: (1) A
Lindblad master equation model of incoherent Markovian noise (both approximate
theory and full simulation), (2) the IBM-Q Qiskit Aer gate-based Kraus
simulator model of incoherent Markovian noise, and (3) a unitary parameter
noise model of stochastic Hamiltonian error.

The overall impact of noise can be studied by its effect on the rate of
fidelity decay of our quantum system. When the noise is unitary, such as
parameter noise, with total evolution $U_{\sigma}$ and noise amplitude
$\sigma$, and a pure state $\ket{\psi}$ is used as an initial condition, the
pure state fidelity is measured by

$f(t)=|\bra{\psi}U_{\sigma^{\prime}}^{-t}U_{\sigma}^{t}\ket{\psi}|^{2}$ (7)

which is also known as the Loschmidt echo. $\sigma^{\prime}$ indicates the
same magnitude as $\sigma$ while being statistically independent. For non-
unitary noise, such as Lindblad noise, one must use the more general density
matrix formulation

$\displaystyle
f(\rho,\sigma)=\left(\Tr\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right)^{2}$ (8)

for ideal (initial) and noisy (final) density matrices $\rho$ and $\sigma$,
respectively. For an initial pure state $\rho=\outerproduct{\psi}{\psi}$ this
simplifies to

$\displaystyle f$ $\displaystyle=\bra{\psi}\sigma\ket{\psi}$ (9)
$\displaystyle=\sum_{i,j}\innerproduct{\psi}{i}\bra{i}\sigma\ket{j}\innerproduct{j}{\psi}=\sum_{i,j}\sigma_{i,j}\rho_{i,j}^{*}.$
(10)

Contrary to previous studies, noise here occurs during both forward and
backward evolution, to connect simulation to experiment. This increases the
physical time per map step by a factor of two, and therefore the observed
fidelity decay rate by the same factor. In all cases we average the fidelity
over the $N$ computational basis state initial conditions $\ket{p}$ in order
to study the average dynamics.

####  2.3.2 Effects of dynamics on Lindblad noise

The Lindblad master equation is the most general type of completely positive
trace-preserving (CPTP) Markovian master equation. It can be written as

$\displaystyle\frac{\partial\rho}{\partial
t}=-i[H,\rho]+\sum_{i}\nu_{i}\left(L_{i}\rho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\rho\\}\right).$ (11)

in dimensionless units where $t=t_{\text{phys}}/T_{\text{step}}$ has units in
the number of map steps, $H=H_{\text{phys}}T_{\text{step}}/\hbar$, and
$\nu_{i}=\nu_{i,\text{phys}}T_{\text{step}}$.

The time $t$ describes forward-only noise, where $\rho$ is evolved forward
with Lindblad noise and then its fidelity is calculated by comparison to a
noiseless evolution. In experiment this would require time-intensive full
state tomography or a qubit-expensive swap test. Instead all results in this
paper use forward-and-back time evolution in which $\rho$ is evolved forward
with noise and then backward ($H\rightarrow-H$) with a different realization
of the noise. This doubles the physical time spent decohering per map step so
that forward-and-back time $t_{\text{fb,phys}}=2t_{\text{phys}}$ at the end of
each map step. Converting this to map step time by
$t_{\text{fb}}=t_{\text{fb,phys}}/(2T_{\text{step}})=t$ results in a fidelity
decay $f_{\text{fb}}(t_{\text{fb}})$ that is different from the forward-only
fidelity $f(t)$. For many processes, like standard Lindblad decay, the doubled
effect of noise doubles the decay rate per map step such that
$f_{\text{fb}}(t_{\text{fb}})=f(2t)$.

Equations in this paper still use forward-only time fidelity $f(t)$ for
simplicity, but when comparing to simulations and experiment that use forward-
and-back time and noise, $f(2t)$ is used. Parameter noise has an exception
that limiting behaviors such as Lyapunov rate and algebraic rate decays appear
to follow $f_{\text{fb}}(t_{\text{fb}})\approx f(t)$ (see Sec. 2.3.5).

Figure 2: Theoretical Lindblad fidelity evolution for $\nu_{1}=0.1$ and
$\nu_{2}=0.2$: (a) comparing the decay of a diffusive state for three, six,
and nine qubits; (b) comparing the decay of localized, superposition, and
diffusive states for six qubits. The average fidelity plateaus at the
uniformly mixed limit $1/N$ when no information about the original state
remains.

We use a simple model of unitary dynamics modified by the non-unitary
decoherence processes of $2n$ single-qubit Lindblad operators $L_{i}$, causing
relaxation to the ground state at rate $\nu_{1}$ and pure dephasing at rate
$\nu_{2}$ for each qubit. These relate to relaxation time $T_{1}$ and total
dephasing time $T_{2}$ via

$\displaystyle 1/T_{1}$ $\displaystyle=$ $\displaystyle\nu_{1}$ $\displaystyle
1/T_{2}$ $\displaystyle=$ $\displaystyle\nu_{1}/2+\nu_{2}/2$ (12)

as can be seen by comparing Appendix C to Ref. [28]. For calculating the
fidelity, the ideal expected density matrix, which is also the initial density
matrix, will be denoted $\rho$, and the noisy evolving density matrix will be
denoted $\sigma$.

The main interaction between dynamics and Lindblad noise is that different
types of dynamics have different typical density matrices that are influenced
by each decoherence effect to different degrees. These include on-diagonal
relaxation (of computational basis states), off-diagonal relaxation (of
superpositions of basis states), and off-diagonal dephasing. In Appendix C, we
analytically determine the decay of fidelity for pure states in the two
dynamical limits of being highly localized and highly diffusive, as well as
the case of uniform superposition, all three of which are summarized in this
section. Briefly, the fidelity of localized pure states is affected only by
on-diagonal relaxation due to their lack of off-diagonal terms, while the
fidelity of diffusive pure states is affected only by off-diagonal relaxation
and dephasing due to a cancellation of on-diagonal effects. The key finding
here is that the $\nu_{2}$ dependence differs between the two limiting cases,
with a continuum of decay rates for the less localized or less entangled cases
in between.

Starting with the fully localized case which has $k=0$, initial conditions
$\ket{\psi}=\ket{p}$ of computational basis states, and fidelity
$f=\sigma_{p,p}$, the evolution of $\ket{\psi}$ due to the diagonal
Hamiltonian has no effect on the density matrix $\sigma$ and therefore on the
fidelity. The Lindblad evolution acts alone, causing relaxation of each qubit
towards the ground state at rate $\nu_{1}$. If one averages the fidelity over
all possible initial states $\ket{p}$, the result is

$\displaystyle f_{L}=$
$\displaystyle\left(1+e^{-2\nu_{\text{single}}t}\right)^{n}/2^{n}$ (13)

so that the initial effective decay rate is

$\displaystyle f_{L}$ $\displaystyle\approx
1-\nu_{\text{eff}}t+O(\nu_{\text{eff}}^{2}t^{2})$ $\displaystyle\approx
e^{-\nu_{\text{eff}}t}$ (14) $\displaystyle\nu_{\text{eff}}$
$\displaystyle=n\nu_{\text{single}}$

where

$\displaystyle\nu_{\text{single}}$ $\displaystyle=$
$\displaystyle\nu_{1}/2=1/(2T_{1})$ (15)

for the localized case. The factor of $1/2$ derives from the average $1/2$
chance of each qubit starting in the excited state. The factor of $n$ due to
$n$ qubits decaying is an important aspect of interpreting measured $T_{1}$
and $T_{2}$ times. Dynamics that are less than fully localized will start to
show effects of the diffusive case described below.

Between the localized and diffusive cases, another interesting case is the
decay of a uniform superposition state, where all $n$ qubits are in $\ket{+}$
states which are unentangled relative to the computational basis. The fidelity
is a product of the single-qubit fidelities. This may look like a simple model
of a diffused state since the probability has spread to all states equally,
but a general diffused state would also be entangled with random, independent
phases on each of the $2^{n}$ states. The uniform superposition state is then
a useful test of the relative roles of spreading probability versus
entanglement in quantum state diffusion. Only half of the single-qubit
fidelity decays, similar to the averaging factor of one half in the localized
case. The same Eq. 13 and Eq. 14 apply with

$\displaystyle\nu_{\text{single}}$ $\displaystyle=$
$\displaystyle(\nu_{1}+\nu_{2})/4=1/(2T_{2}).$ (16)

In the diffusive case the Hamiltonian evolution becomes important, but for the
QSM we have found it can be abstracted away to simplify the problem and still
produce predictions that agree well with simulation. To simplify, note that
chaotic mixing occurs on faster time scales than Lindblad decay. Then the
density matrix can be expected to quickly reach a randomly entangled pure
state of uniformly spread probability, after which the Hamiltonian evolution
has little qualitative effect aside from rapidly changing the precise values
of the phases. Fidelity decay during diffusive dynamics can then be
approximated as the decay of this random pure state under Lindblad evolution.
Averaging over initial conditions contributes to an averaging over these
random phases, causing the phase-dependent terms in $f$ to vanish. Inspecting
$f=\sum_{i,j}\sigma_{i,j}\rho_{i,j}^{*}$ shows that only $\sigma_{i,j}$ terms
that retain their initial phase from $\rho_{i,j}$ survive due to phase
cancellation in $f$. Surprisingly, this causes the qubits to effectively
disentangle for off-diagonal terms $\sigma_{i,j}\rho_{i,j}^{*}$, in a manner
describable by two pieces: a non-trace-preserving single-qubit decay and a
trace-preserving correction. Along the $n$-qubit diagonal the flat
$\rho_{i,i}=1/2^{n}$ causes relaxation effects across
$\sigma_{i,i}\rho_{i,i}^{*}$ to cancel out. This leaves only effects from
dephasing and off-diagonal relaxation, where the random phase averaging
removes the gain to lower states but not the loss from higher states. The
exact fidelity expression for the diffusive case is provided in Eq. C.3.2 in
Appendix C, but its initial decay rate using Eq. 14 is simply

$\displaystyle\nu_{\text{single}}$ $\displaystyle=$
$\displaystyle\nu_{1}/2+\nu_{2}/4$ $\displaystyle=$ $\displaystyle
1/(4T_{1})+1/(2T_{2})\geq 1/(2T_{1}),$

which is strictly faster than both the localized and superposition cases when
the decay rates are non-zero. Random entanglement combines the effects of full
dephasing and partial relaxation for a decay rate that is larger than for any
unentangled case when $\nu_{1}\sim\nu_{2}$.

Figure 2 compares the full analytic fidelity evolution for each dynamical
case. Each fidelity has the expected feature of starting with an exponential
decay that gradually relaxes to the uniformly mixed value of $1/N$. In
superconducting qubits the $\nu_{2}$ rate often dominates, so three types of
dynamics are compared for that case in Fig. 2(b). These fidelity decay rates
do not account for the gate implementation used in experiment, which is partly
rectified in the next section.

####  2.3.3 Gate-based Lindblad model

The theoretical expressions for the localized and diffusive dynamical cases
can now be modified to a gate-based form to better match experimental
observations that the majority of errors occur during certain specific gates.
Since two-qubit gates are the dominant source of error on IBM-Q, one can model
circuit error to lowest order as solely being due to the two-qubit gates
causing an enhanced Lindblad decay for each of the target qubits. Averaged
over the possible two-qubit subsystems on which the gates act during a
circuit, this is approximately described by the dynamics-dependent expressions
previously derived, with number of qubits set to two. If the gates are
performed serially then the serial gate-based fidelity $f_{\text{GB, S}}(t,n)$
can be given in terms of the dynamics-dependent expressions for the Lindblad
fidelity $f_{\text{L}}(t^{\prime},n^{\prime})$ as

$\displaystyle
f_{\text{GB,S}}(t,n)=(f_{\text{L}}(1/M,2))^{Mt}(1-1/2^{n})+1/2^{n}$ (18)

where $M$ is the number of gates per map step, $t=1$ is the time to complete a
map step, and $1/M$ is the time per gate. Note the late-time (average)
fidelity should always approach $1/2^{n}=1/N$. For large $M$ the Lindblad
expression simplifies to $f_{\text{L}}(t\ll 1,n)\approx
e^{-n\nu_{\text{single}}t}$, so

$\displaystyle f_{\text{GB,S}}(t,n)\approx
e^{-2\nu_{\text{single}}t}(1-1/2^{n})+1/2^{n}.$ (19)

This model will be seen in Sec. 4 to qualitatively match experimental results
for $n=3$. Compared to a model that assumes three decaying qubits, when fit to
experimental data this gate-based model measures $\nu_{1}$ and $\nu_{2}$ that
are about $3/2$ larger.

Performing some gates in parallel would increase the average number of
simultaneously targeted qubits to $n_{\text{eff}}\leq n$ and decrease gate
depth to $D=2M/n_{\text{eff}}$. Then the fidelity would be

$\displaystyle
f_{\text{GB,P}}(t,n_{\text{eff}})=(f_{\text{L}}(1/D,n_{\text{eff}}))^{Dt}(1-1/2^{n})+1/2^{n}.$

For large $M$ this simplifies to

$\displaystyle f_{\text{GB,P}}(t,n_{\text{eff}})\approx
e^{-n_{\text{eff}}\nu_{\text{single}}t}(1-1/2^{n})+1/2^{n}.$

Since the physical map step time $T_{\text{step}}$ is $n_{\text{eff}}/2$ times
shorter for the parallel case compared to serial, the two cases yield the same
exponential rate. This just reflects the fact that total fidelity is the
product of gate fidelities when gate error is much greater than idle error.

Figure 3: Comparing Lindblad evolutions from full simulation (solid) to the
theoretical steady state approximation (dashed) for the fully localized
(blue), semi-localized (purple, $k=0.1$) and fully diffusive (red, $k=10.0$)
cases. Using $\nu_{1}=0.1$ and $\nu_{2}=0.2$. (a) $n=3$ ($n_{\text{eff}}=2$)
and (b) $n=6$ ($n_{\text{eff}}=6$).

To test the theoretical predictions of Eq. 18 and 2.3.3, we use QuTiP’s [29]
master equation solver. To partially mimic a circuit model we simulate the
four unitary steps of Eq. 24 sequentially, rather than simulating the whole
Hamiltonian at once. This turns out to be especially important for the
localized case, as described below. For the three-qubit simulation we also
mimic two-qubit gates by only applying Lindblad operators to two qubits at a
time, alternating the qubit pairs [0,1] and [1,2] to mimic the near-linear
connectivity on many IBM-Q devices. This matches the effective two-qubit decay
rate from theory.

The simulation results in Fig. 3 require a modified interpretation of the
theory, as shown for three qubits with serial gates and six qubits with
parallel gates. In the localized case, a straightforward Hamiltonian
simulation would keep the state localized during the whole map step and would
fit well to the theoretical localized decay. But decomposing the evolution to
substeps (or further to gates) reveals that the quantum Fourier transforms
(QFTs) take the localized state to and from a highly entangled state during
half of each map step. This highly entangled state is not fully randomly
entangled like a fully diffusive state, yet the simulation results in Fig. 3
for $k=0.1$ fit well to a half-localized, half-diffusive model where the QFT-
transformed state is assumed to decay at the same rate as a diffusive state,
forming a map step fidelity of $f_{\text{semi-loc}}\approx
f_{\text{loc}}(t/2)f_{\text{dif}}(t/2)\approx(f_{\text{loc}}f_{\text{dif}})^{1/2}$.
This implies a semi-localized decay rate given by

$\displaystyle\nu_{\text{single}}=\nu_{1}/2+\nu_{2}/8.$ (22)

With this modified theoretical prediction, the difference between localized
simulation and semi-localized theory is small. One cause of the remaining
difference in Fig. 3 is the small non-zero $k$, which localizes to several
neighboring states rather than just one, creating some diffusion and
entanglement that enhance the fidelity decay rate. For $n=3$ this accounts for
$\geq 50\%$ of the difference for $t_{\text{fb}}\geq 5$, but for $n=6$ it
explains only $12\%$ of the difference at $t_{\text{fb}}=5$. This suggests the
fidelity difference between a randomly diffused state and the QFT of a
localized state grows as $n$ and $t_{\text{fb}}$ increase.

In the diffusive case the difference between theory and simulation is
generally smaller despite neglect of the exact Hamiltonian evolution. This
suggests that the effect of the precise Hamiltonian evolution on the fidelity
is subdominant.

While this model only includes error from two-qubit gates, it could be
modified to include other gates. The largest correction would come from idle
gates for idling qubits, which have their own decay rates $\nu_{1}$ and
$\nu_{2}$. Idle gates are important because CNOT gate error comes primarily
from the duration of the gate, implying that any qubit idling in parallel with
a CNOT is experiencing a significant fraction of the error. For more details
see Sec. 3.2.2 for gate specifics and Table 3 for experimental fits.

####  2.3.4 IBM Aer gate-based Lindblad simulator

IBM offers several options for artificial noisy simulators. One option is
their device-specific noise models [30] which are calibrated to each device’s
reported errors. Reported $T_{1}$ and $T_{2}$ times and per-gate error rates
from randomized benchmarking are used as input for a combined relaxation,
dephasing, and depolarization error model. However, because these models are
not customizable, they cannot be fit to experimental data.

We instead use a custom IBM Aer noise model. Closest to our gate-based Linblad
model is a simple Aer model based on the function thermal_relaxation_error
with parameters $T_{1}$ and $T_{2}$ and an effective temperature that we set
to zero (the default option) [31]. This function applies a single-qubit
Lindblad decay to each qubit targeted by a gate based on the duration of the
gate, after the gate is done. When $T_{2}>T_{1}$ it applies Kraus operators,
but when $T_{2}<T_{1}$ as is common it probabilistically applies single-qubit
gates instead. We modify this function to always apply Kraus operators. This
model will be used in Fig. 9 to compare the relative importance of specifying
the full gate dynamics as in Aer versus simplifying to the four unitary
substeps Eq. 24 as in the gate-based Lindblad model.

####  2.3.5 Effects of parameter noise

Hamiltonian noise models have been have been studied widely in the quantum
maps literature [32, 33, 6, 34, 35, 7, 20, 36, 37, 19], with noise in a
parameter being a common variant. Hamiltonian noise models use unitary
Markovian errors representing low-frequency or coherent noise. In this
section, we study the model of Ref. [6, 35, 7] in which Hamiltonian evolution
is perturbed stochastically at each discrete map step by $k\rightarrow
k+\Delta k$. The noise $\Delta k$ is drawn from a normal distribution with
standard deviation $\sigma$, where the PDF is given by

$\displaystyle p(\Delta k)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{\Delta
k^{2}}{2\sigma^{2}}}.$ (23)

This model is discussed more thoroughly in Ref.[1]. This produces a similar
effect as the fixed Hamiltonian perturbation discussed in Ref. [19], which
itself resembles static coherent noise. Coherent noise can be divided into a
component that is consistent over many gates, which can often be fixed through
calibration, and a component that varies, which is more difficult to address.
The latter resembles parameter noise.

When this parameter noise model is used, the initial conditions
$\ket{p}=\ket{0},\ket{-N/2}$ are excluded from simulations because their
symmetry with respect to the map causes their noise-averaged fidelity to
asymptotically approach $2/N$ rather than $1/N$. This effect can be misleading
on small systems so is avoided here for clarity.

In the presence of random parameter noise, the fidelity of a chaotic map
initially decays as $e^{-An_{g}t}$ with Fermi golden rule $A\sim\sigma^{2}$
for noise magnitude $\sigma$ and $n_{g}$ gates [20], until the golden rule
breaks down for large $\sigma$ [1]. In the presence of static unitary noise
the fidelity initially decays tangent to $e^{-An_{q}n_{g}^{2}t}$ with
$A\sim\sigma^{2}$ for $n_{q}$ qubits, but at some fraction of the Heisenberg
time this transitions to a faster Gaussian decay of $\ln(f)\sim-t^{2}$ [20].
Quantum dynamical effects, such as Lyapunov and algebraic decay described
below, do not immediately set in, so they do not affect the above initial
decay rates.

It was previously shown that in the diffusive regime a minimum of six qubits
is needed to observe the fidelity decaying at the Lyapunov exponent rate [1].
However in the localized regime an algebraic decay rate is predicted instead
[19, 22, 1]. This effect is shown in the QSM in Fig. 4(b), close to the
predicted $f\propto 1/t$ at intermediate times for dimension $d=1$ and no time
auto-correlation of the noise. Since algebraic decay has less stringent
constraints than Lyapunov rate decay, it might be observable on fewer than six
qubits, though no fewer than four as discussed below.

For parameter noise $f_{\text{fb}}(t_{\text{fb}})$ still indicates forward-
and-back noise and its Fermi golden rule decay has double the decay rate with
respect to map steps, $f_{\text{fb,FGR}}(t_{\text{fb}})\approx
f_{\text{FGR}}(2t)$. But its Lyapunov rate decay has the dependence
$f_{\text{fb,Lyap}}(t_{\text{fb}})\approx e^{-\lambda
t_{\text{fb}}}=e^{-\lambda t}\approx f_{\text{Lyap}}(t)$, showing no increased
rate [1]. A similar effect cannot be tested for in algebraic decay, since
after fitting the initial fidelities $f_{\text{alg}}(t\gg 1)=f_{0}t_{0}/t$ and
$f_{\text{alg}}(2t\gg 1)=f_{0}\,2t_{0}/2t$ are equivalent.

Figure 4: Parameter noise simulations for (a) n=3 and (b) n=6 on log-log plots
to show the algebraic decay rate for localized dynamics. Values of the kick
$k$ and Gaussian noise $\sigma$ are chosen to avoid $k_{\text{eff}}=k+\Delta
k<0$ while keeping the number of map steps small. The resulting localized
values of $k$ are further from zero due to this. Using 1000 and 100
realizations of the stochastic noise plus 6 and 62 initial conditions for (a)
and (b) respectively.

Since the diffusive and localized cases have different qualitative behavior
under parameter noise, a fidelity gap between the two should appear. This
might be confused with the effect of Lindblad noise in a sufficiently noisy
experiment, so it is important to check when algebraic decay might occur.
Simulations with only parameter noise show a gap is only observed on four or
more qubits. On three qubits the gap is in fact reversed, with diffusive
fidelity decaying slightly more slowly than localized fidelity. This is
demonstrated in Fig. 4(a) where a reversed gap of about $3\%$ (absolute)
appears. Interpreting this is complicated by the use of a Guassian parameter
noise distribution which in the case shown causes $2\%$ of map steps to have
$k_{\text{eff}}=k+\Delta k<0$, where anomalous diffusion occurs. However the
reversed gap persists when anomalous diffusion is avoided with small noise
$\sigma\ll k$, so that is not the main cause. A reversed gap or no gap is
observed consistently when comparing $0\leq k\leq 1.9\approx k_{\text{loc}}$
to $k=20$ at noise values $0\leq\sigma\leq 1.0$. Comparing the decay rates to
$1/t$ suggests that the localized case is still close to algebraic decay, but
the diffusive case is limited by the system size to an even slower decay. This
result for three qubits suggests the predictions of Lindblad noise are more
relevant than parameter noise in present day experiments.

Dynamical localization displays Poisson-like energy spacing statistics, and so
indicates locally integrable quantum dynamics, which is the source of the
algebraic fidelity decay [19]. But a dynamical fidelity gap can also occur in
locally integrable quantum systems such as Ref. [11] (see their Fig. 3).

####  2.3.6 Combining Lindblad and parameter noise

The theory of Appendix C demonstrates that Lindblad noise by itself does not
cause the same qualitative behavior as the parameter noise model; e.g.
algebraic decay or exponential decay at the Lyapunov rate. Combining Lindblad
and parameter noise can generate these effects, but if the parameter noise is
small relative to the Lindblad noise the effects may be difficult to observe.

If the noise processes are statistically independent, then Lindblad and
parameter noise fidelities should multiply together. The simulations of the
combined noise model in Fig. 5 support this conclusion for both localized and
chaotic dynamics. In the chaotic case the rate slows for $1\leq t\leq 4$ to
the product of the Lindblad diffusive rate and the Lyapunov rate, matching
pure parameter noise simulations. If the two fidelities instead summed then
the slowest exponential rate, namely the Lindblad rate, would dominate by
outlasting the others, but that would predict a $3\times$ slower decay than
observed. The localized case also supports multiplying fidelities, as during
$4\leq t_{\text{fb}}\leq 7$, near the start of $1/t$ decay in Fig. 4, the
decay remains faster than either the $1/t$ or semi-localized (Eqs. 2.3.3 and
22) analytic decays, rather than one rate dominating. By $t=8$ the $1/N$ term
begins to dominate and slow the decay.

Figure 5: Comparing combined Lindblad and parameter noise simulations (solid)
to theory (dashed) for the semi-localized (purple, $k=1.8$) and diffusive
(red, $k=10.0$) cases. Theory (dashed) includes $1/t$ (blue), Lindblad semi-
localized (purple) and the product of Lindblad diffusive and Lyapunov rate
decays (red). All results use $n=7,\nu_{1}=0.025,\nu_{2}=0.05,\text{ and
}\sigma=0.9$. For Lyapunov rate decay, theory assumes
$f_{\text{fb,Lyap}}(t_{\text{fb}})=e^{-\lambda t_{\text{fb}}}$ as observed.

In Ref. [1] three bounds were derived on the observability of Lyapunov decay
under parameter noise. However if dynamical and Lindblad noise do have
multiplicative effects on the fidelity decay, then one can predict how this
noise combination would cause these bounds to change. Bound (1) should remain
unchanged, as Lindblad noise likely affects both Fermi golden rule decay and
Lyapunov decay equally, causing no change in their competition. Bound (2)
should also remain unchanged, as the transition to localization likely occurs
at the same parameter value. But bound (3) would tighten, since it describes
the limited time steps available to observe Lyapunov decay, which the addition
of Lindblad noise would certainly reduce.

##  3 Implementation

###  3.1 QSM algorithm

There is a natural mapping of the QSM to a qubit-based quantum computer. The
$N$ momentum eigenstates can be mapped to the $2^{n}$ qubit states when
$N=2^{n}$. The unitary $U_{QSM}$ can then be implemented exactly in four steps
[1, 7, 9], written compactly as

$U_{QSM}=U_{kin}U_{QFT}^{-1}U_{pot}U_{QFT}$ (24)

where the operators $U_{kin}=U_{phase}(\hbar)$ and
$U_{pot}=U_{phase}(-k\beta^{2})$ are many-qubit diagonal phase operators in
the position and momentum bases respectively, defined in Eq. 2.1, and
$U_{QFT}$ is the quantum Fourier transform used to alternate between the
momentum and position bases. The diagonal phase operators can be implemented
exactly due to a method for decomposing order-$P$ polynomial terms in the
Hamiltonian into $P$-qubit gates, given in Appendix B [9]. A similar yet
approximate algorithm exists for any quantum map whose Hamiltonian has
separable kinetic and potential energy terms each with a convergent power
series expansion, such as the standard map [9] or kicked Harper model [38]. An
exact representation of trigonometric terms requires ancilla qubits. For the
$U_{QFT}$ steps we use the standard algorithm from IBM’s library, which is
exact, efficient, and requires no ancilla qubits [39, 40].

Using the derivation in Appendix B, the efficient circuit decomposition for
the QSM is

$\displaystyle U_{QFT}=$
$\displaystyle\prod_{j_{1}=0}^{n/2-1}\texttt{SWAP}_{j_{1},n-1-j_{1}}\prod_{j_{2}>j_{1}}^{n-1}\texttt{CP}_{j_{1},j_{2}}(\pi/2^{j_{2}})\prod_{j_{1}=0}^{n-1}\texttt{H}_{j_{1}}$
$\displaystyle U_{pot}=$
$\displaystyle\prod_{j_{2}>j_{1}}^{n-1}\texttt{CP}_{j_{1},j_{2}}(k\beta^{2}2^{j_{1}+j_{2}})\crossproduct$
$\displaystyle\;\prod_{j_{1}=0}^{n-1}\texttt{P}_{j_{1}}(k\beta^{2}2^{2j_{1}-1}-k\beta^{2}N2^{j_{1}-1})$
(25) $\displaystyle U_{kin}=$
$\displaystyle\prod_{j_{2}>j_{1}}^{n-1}\texttt{CP}_{j_{1},j_{2}}(-\hbar
2^{j_{1}+j_{2}})\crossproduct$
$\displaystyle\;\prod_{j_{1}=0}^{n-1}\texttt{P}_{j_{1}}(-\hbar
2^{2j_{1}-1}+\hbar N2^{j_{1}-1})$

where H is a Hadamard gate and P and CP are one-qubit phase and two-qubit
controlled-phase gates respectively. This algorithm is shown for three qubits
in Fig. 6. Note the SWAP gates from $U_{QFT}$ have been eliminated just by
reversing the order of qubits during $U_{pot}$.

Figure 6: Circuit for a single forward map iteration of the three-qubit QSM
algorithm from Eq. 3.1, before conversion to hardware connectivity and
transpilation to native gates. Two-qubit CPHASE gates are used. $U_{pot}$ and
$U_{kin}$ steps use PHASE and CPHASE gates, while $U_{QFT}$ steps use CPHASE
and H gates. Using $k=4.55$.
Device, # qubits
(connectivity) |  (a)
# CNOTS
(raw Qiskit) |  (b)
# CNOTs
(Eq. 3.1+Qiskit) |  (c)
# physical
single-qubit gates
(Eq. 3.1+Qiskit) |  (d)
fidelity $f(1)$
(Eq. 3.1+Qiskit)
---|---|---|---|---
ibmq_manila, n=2 |  4 |  4 |  12-16 |  0.88-0.90
ibmq_manila, n=3 (linear) |  136 |  66 |  28-32 |  0.27-0.37
ibmq_5_yorktown, n=3 (triangular) |  82 |  38 |  26-30 |  0.125 (1/N)
Table 1:  Native gate counts and fidelity for executing each forward-and-back
iteration of the quantum sawtooth map (QSM) experimentally on IBM-Q devices.
To calculate forward-only gate counts as for Fig. 7, divide by two. (a) CNOT
gate count when Qiskit transpiler attempts direct gate decomposition of the
QSM unitary. (b) CNOT gate count when using the efficient algorithm Eq. 3.1
plus transpiler optimization on linear qubit connectivity. (c) Physical
single-qubit gate count, not including virtual RZ gates. Range is over initial
condition and dynamical map parameter $k$. (d) Fidelity as measured by the
one-step Loschmidt echo, partly from Fig. 8. Range is over $k$, varied from
diffusive to localizing dynamics, after averaging over initial conditions.
Experiment on ibmq_5_yorktown performed on October 22, 2020 at 4:41pm EST.
Fidelities include measurement error, and at full decoherence reach $1/N$.

###  3.2 IBM-Q Implementation

####  3.2.1 Circuit optimization

Despite the relative simplicity of the QSM algorithm, executing the circuit in
Eq. 3.1 on IBM-Q’s state-of-the-art hardware quickly accumulates error.
Optimizing the circuit to reduce error can greatly improve the results.

We only attempt to reduce error by decreasing the gate count of the two-qubit
CNOT gate, the largest-error gate. This is motivated by the goal of
benchmarking, which would be complicated by the use of other noise mitigation
techniques such as dynamical decoupling. The focus on CNOT gates is motivated
by the small total contribution of single-qubit gate errors. According to
IBM-Q’s reported RB error estimates, each CNOT gate produces about $30\times$
more error than physical single-qubit gates, namely SX (square-root-of-not)
and X gates. (Z gates are virtual with presumably even smaller error [41].)
This can partially be attributed to the gate duration, with a CNOT gate taking
on average about $10\times$ longer to execute than SX and X gates, which have
identical duration. Combining error per gate with the relative gate counts in
Table 1 suggests single-qubit gates contribute no more than $4\%$ of the total
error in this context. This remains true for increasing qubit count $n$ per
Eq. 3.1.

To determine the base CNOT gate count of the circuit Eq. 3.1 on IBM-Q, first
logic gates are decomposed to native hardware gates. The simple relations for
decomposing are: PHASE($\phi$) = RZ($\phi$) up to an unimportant global phase,
where RZ($\phi$) is a Z-rotation by angle $\phi$; H = RZ($\pi/2$) SX
RZ($\pi/2$); $\texttt{SWAP}_{01}$ =
$\texttt{CNOT}_{01}\texttt{CNOT}_{10}\texttt{CNOT}_{01}$; and lastly the
CPHASE gate is given by
$\texttt{CPHASE}_{01}(\phi)=\texttt{RZ}_{1}(\phi/2)\texttt{CNOT}_{01}\texttt{RZ}_{1}(2\pi-\phi/2)\texttt{CNOT}_{01}\texttt{RZ}_{0}(\phi/2)$.
These decompositions are performed automatically by Qiskit’s transpile
function.

Reducing the CNOT gate count is done with the gate count optimizer that is
automatically applied by the transpile function. This does not guarantee
optimal results but provides a useful improvement. For this the barriers in
Fig. 6 are removed, though barriers between map steps are retained for
scalability with number of map steps. Using the above decompositions of each
gate on Fig. 6, the base CNOT count on three fully-connected qubits is
$12\times 2=24$ CNOT gates. Using transpile per Appendix A applies gate
transformations stochastically to search for gate reductions. Out of 100
attempts the best optimization reduces gate count from $24\rightarrow 19$
CNOTs, for a $21\%$ reduction. Smaller gate counts were found, but were for
incorrect circuits due to a bug described in Appendix A, or used
approximations that did not generalize to other values of $k$.

Retaining consistent gate count across $k$ is crucial to isolating the effect
of dynamics on fidelity decay. This required avoiding specific values of $k$,
such as $2.0$ and $4.5$, and avoiding $k\ll 2\pi/\beta^{2}$ where some CP
gates in $U_{\text{pot}}$ approach identity and the transpiler may eliminate
them from the circuit.

Sparse device connectivity increases gate counts, which allows for further
gate count optimization. For two unconnected qubits, a $\texttt{CNOT}_{02}$
gate is replaced with
$\texttt{SWAP}_{01}\texttt{CNOT}_{12}\texttt{SWAP}_{01}$, transforming one
gate to seven. However each $\texttt{CPHASE}_{02}$ gate only needs one pair of
SWAPs despite its base decomposition to two CNOT gates. So Fig. 6 before
optimization transforms due to linear connectivity from $24\rightarrow 48$
CNOT gates. Gate optimization over 100 transpile attempts reduces from
$48\rightarrow 33$ CNOTs, a $31\%$ reduction. The connectivity-caused CNOTs
have partially canceled with algorithmic CNOTs, reducing the cost of sparse
connectivity.

####  3.2.2 Circuit scaling

The scaling of circuit error is not very relevant to the three-qubit
experiments in this paper, but it is useful for estimating the feasibility of
larger simulations, e.g. that could be used to observe the Lyapunov exponent
[1]. It is discussed here for completeness in understanding the QSM system as
a future benchmarking tool. Gate error is assumed constant for increasing
number of qubits, without accounting for error sources like crosstalk which
are known to increase the average error per CNOT gate [18]. Moreover,
additional scaling may depend on the device under consideration.

When scaling to more than three qubits, performing CNOT gates in parallel
becomes a crucial tool in reducing total error. The degree to which the total
circuit error depends on gate count versus gate depth of parallel gates
depends on how the error of a single CNOT gate $\epsilon_{\text{CNOT}}$
compares to idle error of the same duration $\epsilon_{\text{IDLE}}$. If
$\epsilon_{\text{CNOT}}=\epsilon_{\text{IDLE}}$, then performing two CNOTs in
parallel has identical error as performing a single CNOT with the other qubits
idle, meaning performing gates in parallel is a big gain and gate depth
determines circuit error. If
$\epsilon_{\text{CNOT}}\gg\epsilon_{\text{IDLE}}$, then two CNOTs in parallel
have similar error as two in series, meaning parallelization does little and
gate count determines circuit error. The model in Sec. 2.3.3 assumed
$\epsilon_{\text{CNOT}}\gg\epsilon_{\text{IDLE}}$ to determine gate-specific
error to lowest order. When analyzing circuit scaling, different hypotheses
about this relation can lead to a range of results. This combines with
uncertainty about how parallelized a circuit is for a total effectiveness of
parallelization. This effectiveness reduces error relative to a serial circuit
by a factor between one and $n/2$. The relationship between
$\epsilon_{\text{CNOT}}$ and $\epsilon_{\text{IDLE}}$ is experimentally
estimated in Sec. 4.2.2.

To predict the scaling of circuit error due to CNOT gate count and gate depth
to many qubits, the previous insights can be combined. Base CNOT gate count
from Eq. 3.1 scales as $O(n^{2})$ for $n$ qubits. Worse-case device
connectivity (linear) increases this, as the average CNOT requires an extra
$O(n)$ SWAPs, each requiring three CNOTs, to connect qubits separated by an
average of $O(n)$ connections. Gate count optimization, e.g. using the
transpile function or similar, may reduce this by up to $O(n)$.
Parallelization effectiveness decreases circuit error by between $O(1)$ and
$O(n)$ regardless of connectivity. In summary, a circuit with all-to-all
connectivity has between $O(n)$ and $O(n^{2})$ circuit error scaling due to
CNOT gates, and one with linear connectivity has between $O(n)$ and $O(n^{3})$
scaling. These ranges depend on the effectiveness of gate count optimization
and parallelization. They could be refined by using IBM’s transpiler to see
how well CNOT gates parallelize in the QSM across different topologies.

For the gate count optimization performed by transpile, it is unknown whether
an optimal solution can be found in a scalable fashion. Finding an optimal
gate reduction on all-to-all qubit connectivity is an NP-complete problem
[42]. In fact even for few qubits with linear connectivity transpile struggles
to find optimal results. In Table 1 this is demonstrated by comparing the CNOT
gate count from the known algorithm Eq. 3.1 plus optimization in column (b) to
a naive approach of compiling the raw unitary $U_{QSM}$ directly in column
(a). In each case about one hundred iterations of transpile is attempted. For
$n=3$ the naive approach yields over twice the gate count. Whether this
suboptimal optimization can provide a scalable benefit is unclear, though
focusing on optimizing few-qubit sub-circuits in large algorithms may be the
most practical use case.

Since these scalings are all polynomial, they would all achieve an exponential
speedup relative to the $O(n2^{n})$ operations required by a classical
computer using the fast Fourier transform [7]. However one should keep in mind
the scaling of the measurement step, since measuring the entire final state
would take an exponential number of runs and destroy the speedup. The speedup
is preserved only for measuring collective properties like the Lyapunov
exponent [1], localization length, or anamolous diffusion exponent [7].

####  3.2.3 Choosing a device

At the time of writing IBM-Q freely provides six five-qubit devices, each
containing only linear three-qubit sublattices. Among these we choose to test
ibmq_manila based on its small reported gate error (often $<6.5\mathrm{E}{-3}$
on some three-qubit sublattice) and large quantum volume (32). The now-retired
five-qubit device (ibmq_5_yorktown) had a bow-tie shape with two triangular
sublattices but much larger gate error ($>2.0\mathrm{E}{-2}$), so we performed
experiments on it to explore the trade off between connectivity and gate
error. In Table 1 ibmq_5_yorktown and ibmq_manila are compared, albeit 1.5
years apart, with the latter giving much better performance. Our experiments
on ibmq_santiago contemporary to the ibmq_5_yorktown experiments gave similar
results as the more recent ibmq_manila, indicating that the time lapse has not
had a large effect. A similar result comparing ibmq_5_yorktown and
ibmq_santiago was also found in Ref. [15].

##  4 Experimental results

###  4.1 Dynamical localization

One goal of simulating the QSM on present day hardware is to assess the
hardware’s ability to execute complex dynamical simulations. In Fig. 7 are the
results of simulating the QSM “forward-only” on the ibmq_manila device in both
the localized and diffusive regimes. The effect of dynamics is clearly
apparent, and the localized state’s probability is an informal metric of the
fidelity of the quantum hardware. The raw data shows the localized state
retains the maximum probability among the eight states through $t=7$. However
Fig. 8’s results for simulating forward-and-back retain fidelity through at
least $t_{\text{fb}}=5$, suggesting the forward-only localized case may retain
information about its initial state through at least $t=10$.

The metric of the localized state’s probability was used more thoroughly in
Ref. [15]. Fig. 7 does not easily compare to Ref. [15] because different map
parameters were used, namely $L=1$ in this work and $L=7$ in Ref. [15]. This
significantly changes the eigenstates, independently of $k$, and therefore
changes the degree of localization. Different parameters were used in the
present work to achieve more strongly localized dynamics, as explained in Sec.
2.2.

Figure 7: Dynamics of the QSM for three qubits on the ibmq_manila device,
showing localization ($k=0.1$) and diffusion ($k=4.55$) for $t=1,2,4,8$ and
initial condition $p=-2$ $(\psi=\ket{010})$ for best localization. Compare to
Fig. 1. Hereafter 8192 experimental shots were taken for each $k$, $t$, and
initial condition, giving statistical uncertainty $1/\sqrt{N_{\text{shots}}}$.
Experiment performed on April 11, 2022 at 10:19pm EST. ibmq_manila is re-
calibrated every 1-2 hours to adjust for drift that can increase error.

###  4.2 Fidelity and dynamics

A more quantitative metric of hardware ability is the Loschmidt echo from Sec.
2.3.1. It measures fidelity as the ability to evolve the system forwards then
backwards in time and recover the original state. By choosing simple initial
conditions $\ket{\psi}=\ket{p}$ and working in the momentum basis, the initial
states are just computational basis states and measuring their final fidelity
is just $\bra{p}\sigma\ket{p}=\text{Prob}(\ket{p})$.

By varying the parameter $k$ of the QSM the interaction between experimental
noise and dynamics can be measured, as shown in Fig. 8. As predicted by the
Lindblad noise model in Sec. 2.3.2, but not the parameter noise model in Sec.
2.3.5 for $n=3$, localized and diffusive dynamics have different fidelity
decay rates during the entire simulation. This result shows the existence of
this effect in the presence of hardware noise and suggests the Lindblad model
can capture this important aspect of three-qubit complex quantum simulation.

Note that for $n=3,L=1$ as used here, the predicted transition to full
diffusion should occur at $k_{\text{loc}}\approx 1.87$. In the experiment, the
largest three $k$ values have indistuigishable decays up to a $1.5\%$ absolute
difference despite $k=1.0$ being below the transition threshold. Further
resolution of the observed transition value $k_{\text{loc}}$ requires more
statistics and larger system size. The smallest three $k$ values show a
gradual transition from the strongly localized $k=0.1$ to the weakly localized
$k=0.45$.

|  CNOT error (localized) |  CNOT error (diffusive)
---|---|---
ibmq_manila reported $n=2$ |  5.79E-3 |  5.79E-3
ibmq_manila experiment $n=3$ |  1.76E-2 |  2.59E-2
ratio of experiment:reported |  3.0 |  4.5
Table 2: Comparison of IBM-Q’s reported RB gate error to error extracted from
a three-qubit experiment with localized ($k=0.1$) or diffusive ($k=4.55$)
dynamics. Experimental error $\epsilon$ is calculated from fidelity decay
$f(t)$ via $f(1)=(f(0)-1/2^{n})(1-\epsilon)^{66}+1/2^{n}$.

####  4.2.1 Gate error

The most direct metric for comparing our simulation results to reported
metrics from IBM-Q is the CNOT gate error. This does not rely on any models
from Sec. 2.3, as it is a single-parameter fit (gate error) for each $k$. In
Table 2 error fits are reported, using just the first time step $f(1)$ and the
SPAM error $f(0)$. The fidelity dependence on dynamics is codified here as
gate error dependence on dynamics, with a factor of $1.5\times$ in gate error
between the extreme dynamical cases.

More interestingly, this range of observed errors is $3.0\text{--}4.5\times$
worse than the error reported by IBM-Q on their online “Systems” screen at the
time the experiment was performed [43]. Their reported error comes from
standard two-qubit RB with a depolarizing noise model [18]. What causes such a
difference between these metrics?

The first suspect is the different dynamics caused by RB Clifford gates and
the QSM pseudorandom substeps. Recent research suggests Clifford gates have
different quantum scrambling properties than general unitary dynamics: out-of-
time-order correlators (OTOCs), which are measures of quantum chaos and
scrambling and close relatives of the Loschmidt echo fidelity, reach very
different asymptotic values under Clifford and non-Clifford unitary evolution
[44, 45]. This may relate to the Clifford group being a 2-design on qudits
(and a 3-design on qubits) [44, 46, 47]. Since OTOCs for quantum chaotic
systems only grow exponentially until the Ehrenfest time $\tau_{E}$ [48], the
QSM which has $\tau_{E}\sim 1$ (see Appendix C) saturates quickly. While this
is suggestive of a relation between dynamics and fidelity, quantifying the
effect is beyond the scope of this paper.

Coherent noise on the hardware could also interact with the dynamics to alter
the fidelity decay rate. In Ref. [20] it was shown that in a quantum chaotic
map with static Hamiltonian perturbations the fidelity transitions from Fermi
golden rule decay to a ballistic Gaussian decay at a fraction of the
Heisenberg time. This could further separate the QSM from RB.

Another, better-known suspect is the extra crosstalk from adding a third qubit
relative to two-qubit RB. The role of third qubit crosstalk has been
investigated with simultaneous RB [18], where average CNOT error per gate was
found to increase from $1.63\mathrm{E}{-2}$ for two-qubit RB to
$2.70\mathrm{E}{-2}$ for 2+1-qubit simultaneous RB, a factor increase due to
crosstalk of $1.66$. This is surely relevant here, and the factor may vary
enough across devices or experiments to wholly explain the error difference of
$3.0\text{--}4.5$. A closer study would be useful for teasing out the size of
crosstalk’s contribution.

The effect of dynamics on CNOT error due to Lindblad noise can more fully be
captured by a CNOT process matrix tomography experiment. The diffusion and
entanglement determine which elements of the process matrix are more strongly
enacted in an experiment and therefore which types of error accumulate. The
two-parameter gate-based Lindblad model primarily adds a physically-motivated
heuristic for the fidelity dependence on dynamics observed in Fig. 8. It
provides a balance of utility and efficiency, but sits within the process
matrix picture.

The effect of crosstalk goes beyond a static process matrix, as it causes the
process matrix to depend on the number of qubits. Both active and spectator
qubits can increase error in the CNOT process matrix, meaning even a two-qubit
process matrix must be measured simultaneously to a variety of gates applied
to neighboring qubits in order to correctly predict the fidelity of many-qubit
algorithms [18]. The limited connectivity of IBM-Q devices becomes desirable
here, as it reduces the number of neighbors per qubit which should limit the
magnitude of crosstalk when scaling to many-qubit algorithms. If crosstalk
error can be made to saturate at nearest neighbor effects then both
characterization and NISQ applications become more feasible.

Considering the contribution of each source to the $3.0\text{--}4.5\times$
larger error is crucial to predicting NISQ quantum advantage. The better each
contribution is understood, the more predictive models of error can be
constructed.

Lastly, the depolarizing noise model determined through RB is clearly unable
to capture or explain the fidelity dependence on dynamics. Its single
parameter $\alpha_{2Q}$ of two-qubit gate error measures an average tendency
towards the state $\rho=I/N$, rather than capturing important details of how
different density matrix elements contribute different rates of decay.
Additionally, its focus on averaging over unitary errors, while mathematically
convenient, is perhaps less appropriate than single-qubit decoherent errors
for describing present day superconducting quantum devices.

Figure 8: Average fidelity of the QSM on the ibmq_manila device for varying $k$. Localization occurs below $k_{\text{loc}}\approx 1.87$. Data is averaged over all eight initial computational basis states. Statistical uncertainty per data point is $1/\sqrt{8192*8}\approx 0.4\%$. Number of CNOT gates per forward-and-back step is $M_{\text{CNOT}}=66$. The absolute fidelity gap at $t_{\text{fb}}=1$ between the most localized and most diffusive cases is 10.6%. Experiment performed on January 28, 2022 at 2:44pm EST. Figure 9: Numerical fits of several models to the data in Fig. 8 for the extreme conditions of localized (purple, $k=0.1$) and diffusive (red, $k=4.55$) dynamics. Models described in the main text. |  $T_{1,\text{phys}}$ ($\mu$s) |  $T_{2,\text{phys}}$ ($\mu$s) |  $\nu_{1}$ |  $\nu_{2}$
---|---|---|---|---
ibmq_manila reported $n=1$ |  143 |  37.4 |  0.081 |  0.537
theory fit to experiment $n=3$ |  34.6 |  14.4 |  0.334 $\pm$ 0.016 |  1.271 $\pm$ 0.068
Lindblad simulation fit $n=3$ |  90.1 |  14.3 |  0.128 $\pm$ 0.059 |  1.486 $\pm$ 0.185
Aer simulator fit $n=3$ |  250 |  13.4 |  0.046 $\pm$ 0.019 |  1.68 $\pm$ 0.10
Table 3: Parameter fits from Fig. 9. All values are averaged over three qubits
connected in a line on ibmq_manila. $T_{\text{1,phys}}$ and
$T_{\text{2,phys}}$ values are converted from simulations using
$T_{\text{1,phys}}=T_{\text{step}}/\nu_{1}$ and
$T_{\text{2,phys}}=T_{\text{step}}*2/(\nu_{1}+\nu_{2})$ where
$T_{\text{step}}=33*350\mu s$. Relaxation $\nu_{1}$ and pure dephasing
$\nu_{2}$ are dimensionless decay rates per single-direction map step forwards
or backwards in time. The analytic theory applies continuous Lindblad decay
for two qubits for each of $66$ gates per map step. The Lindblad simulation
continuously applies Lindblad operators to two qubits in alternating pairs
during $8$ algorithm substeps. The modified Aer simulator applies Lindblad
decay to one and two qubits as Kraus operators after each gate.

####  4.2.2 Fitting theory to experiment

To compare the experiment of Fig. 8 to the Lindblad theory in Sec. 2.3.2, it
is convenient to focus on the extreme cases of $k=0.1$ and $k=4.55$. An
account of continuously changing $k$ in a Lindblad model and its effect on
fidelity would require a careful application of localization length to the
Lindblad model, which is beyond the scope of this paper.

Fig. 9 shows numerical fits to experimental fidelity of extreme $k$ for the
gate-based Lindblad theory, the gate-based Lindblad simulations, and an IBM-Q
Aer simulator model, all described in Sec. 2.3. Also included is the default
Aer model from the function AerSimulator.from_backend(FakeManila()), which is
designed to fit the error measured by RB.

The default Aer model underestimates error, but the factor difference between
its error and observed error is only $1.5$ rather than the factor of
$3.0\text{--}4.5$ from Table 2. This is because it draws from a previous
device calibration when the relevant gates were $2.0\text{--}3.0\times$ worse,
artificially improving its accuracy. On average it is likely to be less
accurate for the three-qubit QSM than shown here.

The Lindblad theory fit uses the gate-based model Eq. 19 with rates Eq. 22 for
localized dynamics and Eq. 2.3.2 for diffusive dynamics. The Lindblad
simulation fit applies single-qubit relaxation and dephasing of uniform rates
$\nu_{1}$ and $\nu_{2}$ continuously to two out of three qubits, alternating
between qubits 0 and 1 or 1 and 2. This is applied during the four unitary
steps per Eq. 24 evolved forward-and-back with QuTiP’s mesolve function. The
Aer simulator fit applies two-qubit and single-qubit gates followed by Kraus
operators on the targeted qubits that are determined from decoherence
parameters $T_{1,\text{phys}}$ and $T_{2,\text{phys}}$.

These models of fidelity all have just two parameters: $\nu_{1}$ for
relaxation and $\nu_{2}$ for dephasing, which relate to the physical
decoherence times through Eq. 2.3.2 and
$T_{1}=T_{1,\text{phys}}/T_{\text{step}},T_{2}=T_{2,\text{phys}}/T_{\text{step}}$.
$T_{\text{step}}$ is the time to complete a single forward or backward map
step on hardware, with $T_{\text{step}}=33*350\mu s$ to match $33$ CNOT gates
of average duration $\approx 350\mu s$ on each of the forward and backward
simulation steps. This only considers time spent in CNOT gates, in accordance
with the single-gate-based Lindblad model. The Aer simulator model uses the
same $T_{\text{step}}$ to convert between parameters, despite including
single-qubit gates decaying at the same rate as the CNOT gates, since the
single-qubit gates contribute $<5\%$ to the total duration of gates. Since it
directly fits $T_{1,\text{phys}}$ and $T_{2,\text{phys}}$, this extra error is
in the converted $\nu_{1}$ and $\nu_{2}$ values. Measurement error is
neglected in all models, which requires additional parameters to describe the
probability of each state being misclassified, and which has small effect on
the large exponential decay rates measured here.

The most striking result in Fig. 3 is the consistent values of $T_{2}$ across
models, as it is also the dominant error source. $T_{2}$ measures the
decoherence of entanglement, and the QSM with diffusive dynamics produces
randomly entangled states (see Appendix C), so it is well-suited to benchmark
this noise process.

The variance in $\nu_{2}$ is larger, but seems to be matched by variance in
$\nu_{1}$ to produce a consistent $T_{2}$. This might be explained by
differing weights on $\nu_{1}$ and $\nu_{2}$ across models. In Sec. 2.3.2 the
theoretical model suggests that during diffusive dynamics, averaging over
random phases of entangled states increases the $\nu_{1}$ contribution to the
fidelity decay. This results in $\nu_{1}$ dominating the decay shared by the
localized and diffusive cases and $\nu_{2}$ being the sole cause of difference
between the cases. However three qubits is a small system and has poor random
phase averaging, which could alter the relative roles of $\nu_{1}$ and
$\nu_{2}$. This might explain the difference between the analytic theory
result and the results of the two noise simulations. While the simulation fits
are more trustworthy in this case, for larger system sizes the simplified
theoretical model should converge more closely to the simulation results.

The average $T_{2}$ fit of the three models is $14.0\mu s$, which is
$2.7\times$ smaller than the reported single-qubit idle $T_{2}$ time from
IBM-Q. This is mostly attributable to larger pure dephasing rate $\nu_{2}$,
with an average value $1.48$ that is $2.8\times$ larger than reported. The
$T_{1}$ and $\nu_{1}$ fit values are less consistent, though this is due
mostly to the small values of $\nu_{1}$ being close to zero and therefore
being similar in magnitude to the error bars, making them difficult to
resolve. In particular the observed $T_{1}$ from the Aer simulator fit is
larger than the reported idle value, casting some doubt on the model’s
accuracy with respect to $T_{1}$.

The difference between observed and reported values for $T_{2}$ is less
surprising than the difference in CNOT gate error. The value of measuring this
effective CNOT $T_{2}$ is to provide a more realistic coherence time for
digital quantum simulations and explain heuristically the fidelity dependence
on dynamics. It is also a useful test of the question posed in Sec. 3.2.2: How
does the error of a CNOT gate $\epsilon_{\text{CNOT}}$ compare to the error of
an idle gate of the same duration $\epsilon_{\text{IDLE}}$? Since the
effective CNOT $T_{2}$ time dominates the error and is $2.7\times$ shorter
than the idle time, this suggests $\epsilon_{\text{CNOT}}\sim
3\times\epsilon_{\text{IDLE}}$. This is right in between the limits of
$\epsilon_{\text{CNOT}}=\epsilon_{\text{IDLE}}$ and
$\epsilon_{\text{CNOT}}\gg\epsilon_{\text{IDLE}}$, suggesting the benefit of
gate parallelization must be treated carefully when scaling to larger systems.

Note that IBM-Q reports a $T_{2}$ value measured by the Hahn echo $T_{2E}$.
This experiment is less sensitive to “inhomogeneous broadening”, that is low-
frequency fluctuations that can cause large differences from trial to trial
[28]. Since $T_{2E}\geq T_{2}$, comparing our smaller $T_{2}$ fits from
running an algorithm to IBM-Q’s larger reported $T_{2E}$ may overstate the
difference between CNOT $T_{2}$ and idle $T_{2}$, which would suggest a
difference factor less than $2.7$.

##  5 Conclusion

In this study an important relationship between the dynamics of a digital
quantum simulation and its fidelity under Markovian noise was illustrated.
Incoherent Lindblad noise was shown to produce a particularly strong effect in
the few-qubit experiments that are currently possible. As the dynamics are
changed from localized, to diffusive but unentangled, to diffusive and
entangled, the corresponding fidelity decays at an effective $T_{1}$ rate, an
effective $T_{2}$ rate, and a combined rate that is faster than either,
respectively.

A Lindblad model with single-qubit noise processes was used to predict this
effect, and several related models were used to fit an experimental
demonstration. They contain far less information than full process tomography,
but more than a depolarizing noise model. This type of physically motivated
model is ideal for capturing the effects of dynamics at minimal experimental
cost.

To connect to the circuit model, Lindblad noise was assumed to only occur
during CNOT gates, which have the largest error. This retains the expected
fidelity dependence on gate count while expanding the gate error model from
one parameter to two ($T_{1}$ and $T_{2}$). When further details of the
unitary gate decomposition were included, it was found that the dynamics
within a simulation time step did have a large impact on fidelity. However the
difference between decomposing a time step into four unitary substeps versus
the 60+ single- and two-qubit hardware gates was relatively small.

From experiment, first the average error per CNOT gate was estimated. Due to
the fidelity dependence on dynamics, the average error varies by a factor of
$1.5\times$ over the dynamics. The observed error is $3.0-4.5\times$ greater
than reported by IBM-Q’s method of randomized benchmarking with a depolarizing
noise model. This discrepancy is largely due to the increase in crosstalk
error from two qubits to three. However, how much the different dynamics
between RB and the QSM contribute to this error difference is not yet clear.

All models that were fit to experiment agree well on the effective $T_{2}$
time, which is dominant, but diverge on the $T_{1}$ time, which has a weaker
effect and is difficult to resolve. The effective CNOT $T_{2}$ time is
$2.7\times$ shorter than that of an idle qubit, suggesting a $\sim 3\times$
larger error than an idle gate for diffused states.

Retrieving useful experimental fidelities for complex Hamiltonian simulation
was enabled by the highly efficient gate decomposition of the quantum sawtooth
map (QSM). It combines the exactness of quantum maps, which do not require
Trotterization, with the efficient algorithm for low-order polynomial
Hamiltonians. On IBM-Q the QSM required only 33 CNOT gates per evolution time
step and 66 CNOT gates per fidelity time step. This enabled seven time steps
of clearly localized dynamics and five time steps of forward-and-back fidelity
decay.

A model of low-frequency parameter noise was also studied that is relevant to
future experiments on more qubits. This model was previously used to show a
minimum of six qubits is required to observe the Lyapunov exponent in the QSM.
Here it was combined with Lindblad noise to show that the Lyapunov exponent
can still be observed, but will require seven or more qubits depending on the
strength of Lindblad noise.

The existence of a large fidelity gap between different dynamics in digital
quantum simulations also raises the question of whether arbitrary algorithms
may have dynamics that affect their fidelity. In particular, the degrees of
entanglement, superposition, and randomness during an algorithm could suggest
particular fidelity dependence on the effective $T_{1}$ and $T_{2}$. This
could be crucial for anticipating quantum advantage by better extrapolating
limited device characterization to predict the fidelities of NISQ algorithms.

##  6 Acknowledgements

The authors thank the Quantum Leap group at LLNL for stimulating discussions
and ideas that improved the manuscript, including Vasily Geyko, Frank R.
Graziani, Stephen B. Libby, and Yuan Shi. We also thank Jonathan L. DuBois,
Kristin M. Beck, Alessandro R. Castelli, and Yaniv J. Rosen of the Quantum
Coherent Device Physics group at LLNL, Kyle A. Wendt of the Nuclear Data and
Theory group at LLNL, and Robert Tyler Sutherland at UT San Antonio for their
insights into quantum hardware and noise processes.

This work was performed by LLNL under the auspices of the U. S. DOE under
Contract DE-AC52-07NA27344 and was supported by the DOE Office of Fusion
Energy Sciences “Quantum Leap for Fusion Energy Sciences” project FWP-SCW1680
and by LLNL Laboratory Directed Research and Development project 19-FS-078.

## References

  * [1] Max D Porter and Ilon Joseph.  Observability of fidelity decay at the lyapunov rate in few-qubit quantum simulations.  arXiv preprint arXiv:2110.07767, 2021.
  * [2] Ryan Babbush, Jarrod R McClean, Michael Newman, Craig Gidney, Sergio Boixo, and Hartmut Neven.  Focus beyond quadratic speedups for error-corrected quantum advantage.  PRX Quantum, 2(1):010103, 2021.
  * [3] Sergio Boixo, Sergei V Isakov, Vadim N Smelyanskiy, Ryan Babbush, Nan Ding, Zhang Jiang, Michael J Bremner, John M Martinis, and Hartmut Neven.  Characterizing quantum supremacy in near-term devices.  Nature Physics, 14(6):595–600, 2018.
  * [4] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando GSL Brandao, David A Buell, et al.  Quantum supremacy using a programmable superconducting processor.  Nature, 574(7779):505–510, 2019.
  * [5] Ryan Babbush.  Google quantum summer symposium 2021: Google’s perspective on the viable applications of early fault-tolerant quantum computers.  https://www.youtube.com/watch?v=-fcQt5C2XGY&list=PLpO2pyKisOjL7JdCjzMeOY1w3TnwTkBT-&index=16, 2021\.  Accessed: 2021-09-27.
  * [6] Giuliano Benenti and Giulio Casati.  Quantum-classical correspondence in perturbed chaotic systems.  Physical Review E, 65(6):066205, 2002.
  * [7] Giuliano Benenti, Giulio Casati, and Simone Montangero.  Quantum computing and information extraction for dynamical quantum systems.  Quantum Information Processing, 3(1):273–293, 2004.
  * [8] David Poulin, Robin Blume-Kohout, Raymond Laflamme, and Harold Ollivier.  Exponential speedup with a single bit of quantum information: Measuring the average fidelity decay.  Physical Review Letters, 92(17):177906, 2004.
  * [9] Bertrand Georgeot and Dima L Shepelyansky.  Exponential gain in quantum computing of quantum chaos and localization.  Physical Review Letters, 86(13):2890, 2001.
  * [10] Asher Peres.  Stability of quantum motion in chaotic and regular systems.  Physical Review A, 30(4):1610, 1984.
  * [11] Nathan K Lysne, Kevin W Kuper, Pablo M Poggi, Ivan H Deutsch, and Poul S Jessen.  Small, highly accurate quantum processor for intermediate-depth quantum simulations.  Physical Review Letters, 124(23):230501, 2020.
  * [12] Lukas M Sieberer, Tobias Olsacher, Andreas Elben, Markus Heyl, Philipp Hauke, Fritz Haake, and Peter Zoller.  Digital quantum simulation, trotter errors, and quantum chaos of the kicked top.  npj Quantum Information, 5(1):1–11, 2019.
  * [13] Markus Heyl, Philipp Hauke, and Peter Zoller.  Quantum localization bounds trotter errors in digital quantum simulation.  Science Advances, 5(4):eaau8342, 2019.
  * [14] Michael K Henry, Joseph Emerson, Rudy Martinez, and David G Cory.  Localization in the quantum sawtooth map emulated on a quantum-information processor.  Physical Review A, 74(6):062317, 2006.
  * [15] Andrea Pizzamiglio, Su Yeon Chang, Maria Bondani, Simone Montangero, Dario Gerace, and Giuliano Benenti.  Dynamical localization simulated on actual quantum hardware.  Entropy, 23(6):654, 2021.
  * [16] Robin Blume-Kohout, John King Gamble, Erik Nielsen, Jonathan Mizrahi, Jonathan D Sterk, and Peter Maunz.  Robust, self-consistent, closed-form tomography of quantum logic gates on a trapped ion qubit.  arXiv preprint arXiv:1310.4492, 2013.
  * [17] Erik Nielsen, Kenneth Rudinger, Timothy Proctor, Antonio Russo, Kevin Young, and Robin Blume-Kohout.  Probing quantum processor performance with pygsti.  Quantum Science and Technology, 5(4):044002, 2020.
  * [18] David C McKay, Sarah Sheldon, John A Smolin, Jerry M Chow, and Jay M Gambetta.  Three-qubit randomized benchmarking.  Physical Review Letters, 122(20):200502, 2019.
  * [19] Philippe Jacquod and Cyril Petitjean.  Decoherence, entanglement and irreversibility in quantum dynamical systems with few degrees of freedom.  Advances in Physics, 58(2):67–196, 2009.
  * [20] Klaus M Frahm, Robert Fleckinger, and Dima L Shepelyansky.  Quantum chaos and random matrix theory for fidelity decay in quantum computations with static imperfections.  The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics, 29(1):139–155, 2004.
  * [21] Eliot Kapit.  The upside of noise: engineered dissipation as a resource in superconducting circuits.  Quantum Science and Technology, 2(3):033002, 2017.
  * [22] Fernando M Cucchietti.  The loschmidt echo in classically chaotic systems: Quantum chaos, irreversibility and decoherence.  arXiv preprint quant-ph/0410121, 2004.
  * [23] Heinz-Peter Breuer, Francesco Petruccione, et al.  The theory of open quantum systems.  Oxford University Press on Demand, 2002.
  * [24] Isaac L Chuang and Michael A Nielsen.  Prescription for experimental determination of the dynamics of a quantum black box.  Journal of Modern Optics, 44(11-12):2455–2467, 1997.
  * [25] Kyle Wendt.  Private communication, January 2022.
  * [26] Dima Shepelyansky.  Scholarpedia: Ehrenfest time and chaos.  doi:10.4249/scholarpedia.55031, 2020.  Accessed: 2022-05-20.
  * [27] Jan Šuntajs, Janez Bonča, Tomaž Prosen, and Lev Vidmar.  Quantum chaos challenges many-body localization.  Physical Review E, 102(6):062144, 2020.
  * [28] Philip Krantz, Morten Kjaergaard, Fei Yan, Terry P Orlando, Simon Gustavsson, and William D Oliver.  A quantum engineer’s guide to superconducting qubits.  Applied Physics Reviews, 6(2):021318, 2019.
  * [29] J Robert Johansson, Paul D Nation, and Franco Nori.  Qutip: An open-source python framework for the dynamics of open quantum systems.  Computer Physics Communications, 183(8):1760–1772, 2012.
  * [30] Device backend noise model simulations.  https://qiskit.org/documentation/tutorials/simulators/2_device_noise_simulation.html, 2021\.  Accessed: 2021-11-01.
  * [31] IBM.  Source code for thermal relaxation error.  https://qiskit.org/documentation/_modules/qiskit/providers/aer/noise/errors/standard_errors.html#thermal_relaxation_error, 2021\.  Accessed: 2022-05-05.
  * [32] Philippe Jacquod, Peter G Silvestrov, and Carlo WJ Beenakker.  Golden rule decay versus lyapunov decay of the quantum loschmidt echo.  Physical Review E, 64(5):055203, 2001.
  * [33] Giuliano Benenti, Giulio Casati, Simone Montangero, and Dima L Shepelyansky.  Efficient quantum computing of complex dynamics.  Physical Review Letters, 87(22):227901, 2001.
  * [34] Giuliano Benenti, Giulio Casati, Simone Montangero, and Dima L Shepelyansky.  Eigenstates of an operating quantum computer: hypersensitivity to static imperfections.  The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics, 20(2):293–296, 2002.
  * [35] Giuliano Benenti, Giulio Casati, Simone Montangero, and Dima L Shepelyansky.  Dynamical localization simulated on a few-qubit quantum computer.  Physical Review A, 67(5):052312, 2003.
  * [36] Wen-ge Wang, Giulio Casati, and Baowen Li.  Stability of quantum motion: Beyond fermi-golden-rule and lyapunov decay.  Physical Review E, 69(2):025201, 2004.
  * [37] Thomas Gorin, Tomaž Prosen, Thomas H Seligman, and Marko Žnidarič.  Dynamics of loschmidt echoes and fidelity decay.  Physics Reports, 435(2-5):33–156, 2006.
  * [38] Benjamin Lévi and Bertrand Georgeot.  Quantum computation of a complex system: The kicked harper model.  Physical Review E, 70(5):056218, 2004.
  * [39] Michael A Nielsen and Isaac L Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2010.
  * [40] Quantum fourier transform.  https://qiskit.org/textbook/ch-algorithms/quantum-fourier-transform.html, 2021\.  Accessed: 2021-11-01.
  * [41] David C McKay, Christopher J Wood, Sarah Sheldon, Jerry M Chow, and Jay M Gambetta.  Efficient z gates for quantum computing.  Physical Review A, 96(2):022330, 2017.
  * [42] Adi Botea, Akihiro Kishimoto, and Radu Marinescu.  On the complexity of quantum circuit compilation.  In Eleventh annual symposium on combinatorial search, 2018.
  * [43] IBM.  IBM-Q Systems.  https://quantum-computing.ibm.com/services?services=systems&systems=yours, 2022\.  Accessed: 2022-01-28.
  * [44] Daniel A Roberts and Beni Yoshida.  Chaos and complexity by design.  Journal of High Energy Physics, 2017(4):1–64, 2017.
  * [45] Lorenzo Leone, Salvatore FE Oliviero, and Alioscia Hamma.  Isospectral twirling and quantum chaos.  Entropy, 23(8):1073, 2021.
  * [46] Zak Webb.  The clifford group forms a unitary 3-design.  arXiv preprint arXiv:1510.02769, 2015.
  * [47] Huangjun Zhu.  Multiqubit clifford groups are unitary 3-designs.  Physical Review A, 96(6):062336, 2017.
  * [48] Koji Hashimoto, Keiju Murata, and Ryosuke Yoshii.  Out-of-time-order correlators in quantum mechanics.  Journal of High Energy Physics, 2017(10):1–31, 2017.
  * [49] Max D Porter.  Quantum computing stack exchange question.  https://quantumcomputing.stackexchange.com/questions/23458/qubit-identities-get-swapped-in-ibm-qiskit, 2021\.  Accessed: 2022-05-05.

##  Appendix A Optimizing performance of a complex algorithm on IBM-Q

Despite using an efficient algorithm, various difficulties can harm
performance. Here are some insights for avoiding pitfalls:

  * •

Optimize for qubit connectivity: Circuits that are transpiled to devices with
sparse connectivity will add SWAP gates to connect distant qubits. This allows
for partial cancellation between the algorithm CNOT gates and the swapping
CNOT gates. The automated approach to this is to use a stochastic transpiler
as offered by Qiskit. To enable this level of optimization, first set the
transpile function’s option optimization_level from =1 (the default) to =3
(the maximum). Next perform a trial-and-error optimization over the transpiler
seed by iterating over the option seed_transpiler=seed and looking for the
lowest resulting gate count of CNOTs using circ.count_ops() on the transpiled
circuit. Check that the transpiled circuit produces the expected outcome in a
noiseless Aer simulator, due to the bug discussed below. Consider saving the
result as seed behavior may change. This procedure does not directly scale to
many qubits, but on few-qubit subspaces it may be useful at NISQ scales, and
may eventually be automated for real applications.

  * •

Optimize algorithms: Look for small optimizations that can be done by hand.
For example the QFT algorithm and its inverse require SWAP gates to reverse
the qubit ordering. But this can be canceled with the argument do_swaps=False
and replaced with a free, manual reversal of qubit ordering for gates in
between the QFT pair.

  * •

Optimize coherence times: Noise properties of backend devices can drift and
worsen from their reported values. IBM-Q routinely recalibrates their systems,
adjusting qubit and gate pulse properties to reduce error. Each device is
currently recalibrated between hourly and weekly, with time since last
calibration reported on each device’s status screen online. Timing experiments
to occur soon after calibration can help achieve optimal performance for
devices that calibrate less than hourly. Choosing the connected set of qubits
with the lowest-error gates is also important and should be repeated after a
calibration that re-characterizes them, though these occur less often.

  * •

Avoid bugs: A major bug in Qiskit (as of this publication) causes qubit
identities to be swapped at the end of a circuit when transpiled without a
measurement, as the measurements absorb swaps even when absent [49]. This
causes problems when repeating an algorithm step as done here, as it is far
more efficient to transpile a single iteration (without measurement) and
repeat it rather than transpile a many-iteration circuit all at once. Since
this error occurs stochastically, the transpiled circuit should be carefully
checked for correct final qubit identities. Seeds that produce incorrect
identities can be avoided, or corrective SWAP gates may be added and even re-
transpiled to undo the error.

##  Appendix B Gate decomposition of polynomial Hamiltonian evolution

The first step in converting the operators of Eq. 2.1 into a series of common
two-qubit gates is knowing how to implement $\hat{p}$ or $\hat{q}$ (which have
identical form in their respective bases). Representing a single-body state of
quantum number $p^{\prime}=0,...,N-1$ on $n=\log(N)$ qubits can be
accomplished via the binary representation [9]

$\ket{p^{\prime}}\rightarrow
p^{\prime}=\sum_{j=0}^{n-1}\alpha_{j}2^{j}\rightarrow\ket{\alpha_{n-1}...\alpha_{0}}$
(26)

where $\alpha_{j}=0$ or $1$ for each qubit $j$ depending on the value of
$p^{\prime}$. One can express diagonal operators similarly. The operator
$U_{p^{\prime}}=e^{-i\hbar\hat{p}^{\prime 2}/2}$ can be expressed in the
$p^{\prime}$ basis by noting

$\displaystyle p^{\prime 2}$
$\displaystyle=\sum_{j_{1},j_{2}}\alpha_{j_{1}}\alpha_{j_{2}}2^{j_{1}+j_{2}}\text{\quad
and so}$ (27) $\displaystyle e^{-i\hbar p^{\prime 2}/2}$
$\displaystyle=\prod_{j_{1},j_{2}}e^{-i\hbar\alpha_{j_{1}}\alpha_{j_{2}}2^{j_{1}+j_{2}-1}}$
$\displaystyle=\prod_{j_{1}}e^{-i\hbar\alpha_{j_{1}}^{2}2^{2j_{1}-1}}\prod_{j_{1}\neq
j_{2}}e^{-i\hbar\alpha_{j_{1}}\alpha_{j_{2}}2^{j_{1}+j_{2}-1}}$ (28)
$\displaystyle e^{-i\hbar\hat{p}^{\prime 2}/2}$
$\displaystyle\rightarrow\left[\begin{pmatrix}1&\\\ &e^{-i\hbar
2^{2n-3}}\end{pmatrix}\otimes\dots\begin{pmatrix}1&\\\ &e^{-i\hbar
2^{-1}}\end{pmatrix}\right]\times$
$\displaystyle\quad\,\,\,\left[\bigotimes^{n-2}I\otimes\begin{pmatrix}1&&&\\\
&1&&\\\ &&1&\\\ &&&e^{-2i\hbar 2^{0}}\end{pmatrix}\right]\times\dots.$ (29)

This describes PHASE gates on every qubit $j_{1}$ and CPHASE gates on every
pair of qubits $j_{1},j_{2}$ with phase $-\hbar 2^{j_{1}+j_{2}-1}$ (applied
twice, for $j_{1}>j_{2}$ and $j_{1}<j_{2}$). Since the desired operator
$U_{kin}$ is shifted by $N/2$ from $U_{p^{\prime}}$, because
$p=-N/2,...,(N-1)/2$, it gets modified to

$\displaystyle U_{kin}\equiv U_{phase}(\hbar)$
$\displaystyle=e^{-i\hbar(\hat{p}^{\prime}-N/2)^{2}/2}$ $\displaystyle\sim
e^{-i\hbar\hat{p}^{\prime 2}/2}e^{i\hbar N\hat{p}^{\prime}/2}$ (30)

with unobservable global phase neglected. This contributes further PHASE gates
described by

$e^{i\hbar N\hat{p}^{\prime}/2}=\prod_{j}e^{i\hbar N\alpha_{j}2^{j-1}}.$ (31)

This describes the $U_{kin}$ gate, but the $U_{pot}$ gate is identical in its
own basis with $\hbar$ replaced by $-k\beta^{2}$ and $\hat{p}$ by $\hat{q}$
with the same $q^{\prime}=0,...,N-1$:

$\displaystyle U_{pot}\equiv U_{phase}(-k\beta^{2})$
$\displaystyle=e^{ik\beta^{2}(\hat{q}^{\prime}-N/2)^{2}/2}$ $\displaystyle\sim
e^{ik\beta^{2}\hat{q}^{\prime 2}/2}e^{-ik\beta^{2}N\hat{q}^{\prime}/2}.$ (32)

Switching between bases as in Eq. 24 is needed to execute both $U_{kin}$ and
$U_{pot}$ with this efficient decomposition. The full algorithm is converted
to gates in Eq. 3.1.

##  Appendix C Fidelity decay of localized and diffusive dynamics under
Lindblad noise

###  C.1 Fidelity of pure states

The quantum sawooth map exhibits dynamics that range from localized to
diffusive (chaotic). To study its interaction with Lindblad noise from a
theoretical standpoint, we can make some reasonable simplifying assumptions.
The main assumption is that decay of an evolving density matrix is
approximately the decay of a pure state with distribution and entanglement
typical of the evolving state. In the localized limit this is exactly true. In
the diffusive case this requires that the initial evolution to a typical
diffusive state is fast relative to Lindblad decay and that after a typical
diffusive state is reached the ongoing Hamiltonian evolution can be neglected.
The goal is to calculate the fidelity of the decaying state relative to the
initial pure state over time. The pure state assumption allows a simplified
form. First the general fidelity of corrupted density matrix $\sigma$ relative
to pure state density matrix $\rho$ is

$\displaystyle f=\left(\Tr\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right)^{2}.$
(33)

But a pure state has the nice property that

$\displaystyle\rho^{2}$
$\displaystyle=\ket{\psi}\innerproduct{\psi}{\psi}\bra{\psi}=\rho=\sqrt{\rho}$
(34)

which simplifies the fidelity to

$\displaystyle f$ $\displaystyle=\bra{\psi}\sigma\ket{\psi}$
$\displaystyle=\sum_{i,j}\sigma_{i,j}\rho_{j,i}=\sum_{i,j}\sigma_{i,j}\rho_{i,j}^{*},$
(35)

the Frobenius inner product of the density matrices. It will be useful to
consider the elements of the Hadamard element-wise product

$\displaystyle(\sigma\circ\rho^{*})_{ij}=\sigma_{i,j}\rho_{i,j}^{*}$ (36)

which has elements that sum to the fidelity.

###  C.2 Localized dynamics

The quantum sawtooth map has a limit of localized dynamics at $k\rightarrow
0$. In this limit the Hamiltonian becomes a diagonal phase matrix
$H=\text{diag}(H_{0}\,H_{1})$. This only evolves the off-diagonal elements
$\rho_{i,j}$ of the density matrix $\rho$, as seen on a single-qubit example

$\displaystyle[H,\rho]$
$\displaystyle=\begin{pmatrix}0&(H_{0}-H_{1})\rho_{01}\\\
(H_{1}-H_{0})\rho_{10}&0\end{pmatrix}.$ (37)

If we assume an initial state in the computational basis
$\rho=\text{diag}(0...1...0)$ then the Hamiltonian evolution has no effect on
$\rho$.

####  C.2.1 Single-qubit decomposition

The effect of Lindblad noise on the localized state can be trivially
decomposed into single-qubit effects, since both the Lindblad operators and
density matrix decompose.

The Lindblad equation Eq. 11 in the main text describes the noisy time
evolution of initial state $\rho$, with Lindblad collapse operators for
$T_{1}$ relaxation and $T_{\varphi}$ pure dephasing processes [28] on qubit
$j$ being

$\displaystyle L_{1,j}$ $\displaystyle=I\otimes\dots\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\dots\otimes I$ (38) $\displaystyle L_{2,j}$
$\displaystyle=I\otimes\dots\begin{pmatrix}0&0\\\ 0&1\end{pmatrix}\dots\otimes
I$ (39)

respectively. For a single unentangled qubit, the solution $\sigma$ to the
Lindblad equation when $H$ has no effect evolves as

$\displaystyle\dot{\sigma}=\begin{pmatrix}\nu_{1}\sigma_{11}&-\frac{\nu_{1}+\nu_{2}}{2}\sigma_{01}\\\
-\frac{\nu_{1}+\nu_{2}}{2}\sigma_{10}&-\nu_{1}\sigma_{11}\end{pmatrix}.$ (40)

with decay rates that obey $\nu_{1}\equiv 1/T_{1}$, $\nu_{2}\equiv
1/T_{\varphi}$, and $\nu_{1}/2+\nu_{2}/2=1/T_{2}$. Solving this differential
equation is simple, as three of the elements depend only on themselves,
yielding exponential decays, and the $\sigma_{00}$ element determined by trace
preservation $\Tr(\sigma)=1$.

Localization simplifies both $\sigma$ and the fidelity calculation. If the
qubit starts in $\ket{\psi}=\ket{1}$, then

$\displaystyle\sigma=\begin{pmatrix}1-e^{-\nu_{1}t}&0\\\
0&e^{-\nu_{1}t}\end{pmatrix}$ (41)

and $f=\sigma_{11}=e^{-\nu_{1}t}$. If it starts in $\ket{\psi}=\ket{0}$ then
$f=1$. No off-diagonal terms means no $\nu_{2}$ dependence.

####  C.2.2 Average fidelity

It is useful to average the localized fidelity over all initial conditions for
comparison to the diffusive case later. The fidelity for $n$ qubits decomposes
into the product of single-qubit fidelities since both $\rho$ and $L_{i,j}$
decompose. $j$ excited qubits then decay at an exponential rate $-j\nu_{1}$.
Averaging over initial conditions includes $n$ states decaying at $-\nu_{1}$,
($n$ choose 2) states decaying at $-2\nu_{1}$, and so on. These sum to

$\displaystyle f_{ave}$ $\displaystyle=\frac{1}{2^{n}}\sum_{j=0}^{n}{n\choose
j}e^{-j\nu_{1}t}$ $\displaystyle=\left(\frac{1+e^{-\nu_{1}t}}{2}\right)^{n}$
$\displaystyle\approx 1-\frac{n}{2}\nu_{1}t+O(n^{2}\nu_{1}^{2}t^{2})$
$\displaystyle\approx e^{-n(\nu_{1}/2)t}.$ (42)

The single-qubit decay rate is seen to be $-\nu_{1}/2$, as expected.

###  C.3 Diffusive dynamics

In the other limit of the quantum sawtooth map, large $k$ leads to chaotic,
diffusive dynamics. Hamiltonian evolution now has an effect, but will be
simplified by two assumptions.

The first assumption is that the fast Hamiltonian evolution reaches a fully
diffused state quickly enough to be considered immediate, so that the analysis
begins with a randomly entangled state. To justify this, there are two
relevant time scales for the initial wave function spreading: the Ehrenfest
time $\tau_{E}$ and the Heisenberg time $\tau_{H}$ (see Sec. 2.2). The
Ehrenfest time is how long the quantum dynamics closely match the classical
dynamics [26]. For a chaotic system this implies rapid spreading at the rate
$e^{\lambda t}$. This must end after reaching the finite system size $N$. This
implies an Ehrenfest time of $\tau_{E}\sim\ln(N)/\lambda$. To avoid
localization $K\gtrsim N^{-2/5}$ when $L=1$, implying $\lambda\gtrsim
N^{-1/5}$ and a bound on Ehrenfest time of $\tau_{E}\lesssim\ln(N)N^{1/5}$.
For small system sizes this is clearly small, and simulations suggest
$\tau_{E}\sim 1$.

Between the Ehrenfest time and Heisenberg time a much slower diffusion occurs.
But after the Heisenberg time the effects of quantum interference are fully
realized and dynamical localization occurs, halting diffusion. The Heisenberg
time for localization is $\tau_{H}\sim\ell$ [7, 26], which if extended to the
diffusive regime suggests $\tau_{H}\gtrsim N$. Since the Ehrenfest time
represents the bulk of the spreading and occurs on a short time scale, this
justifies the assumption of an immediate fully diffused state.

The second assumption is that Hamiltonian evolution after full diffusion has
little effect on the fidelity decay and can be neglected to simplify the
Lindblad evolution. This assumption is supported by numerical simulation
results in Fig. 3.

A “typical” diffused state will be considered, with random features to be
averaged over. This diffused state is statistically independent of the initial
state before Hamiltonian evolution, so an explicit averaging over initial
conditions will not be needed. While both amplitudes and phases should be
random, a further simplification is to assume that the amplitudes are nearly
uniform after the diffusive process occurs, leaving only random phases. This
assumption is justified by ergodicity of the dynamics.

After these simplifications, a diffusive state is characterized by two
features: statistically uniform probability distribution and random
entanglement.

####  C.3.1 Unentangled state

To test uniform probability and random entanglement independently, first
consider a uniform state that is not entangled. Such a state is constructed by
initializing in the ground state and applying a Hadmard gate to each qubit. In
this case one can again decompose to single-qubit evolution Eq. 40, but acting
on initial state

$\displaystyle\ket{\psi}$ $\displaystyle=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$
$\displaystyle\rho$ $\displaystyle=\frac{1}{2}\begin{pmatrix}1&1\\\
1&1\end{pmatrix}.$ (43)

The Lindblad solution for initial condition $\rho$ is

$\displaystyle\sigma=\frac{1}{2}\begin{pmatrix}2-e^{-\nu_{1}t}&e^{-\frac{\nu_{1}+\nu_{2}}{2}t}\\\
e^{-\frac{\nu_{1}+\nu_{2}}{2}t}&e^{-\nu_{1}t}\end{pmatrix}.$ (44)

The fidelity is again the product of single-qubit fidelities, with each qubit
contributing
$f=\frac{1}{2}\sum_{i,j}\sigma_{i,j}=\frac{1}{2}(1+e^{-\frac{\nu_{1}+\nu_{2}}{2}t})$
where exponentials on the diagonal cancel. The $n$-qubit fidelity for any
initial condition is

$\displaystyle f$
$\displaystyle=\left(\frac{1+e^{-\frac{\nu_{1}+\nu_{2}}{2}t}}{2}\right)^{n}$
$\displaystyle\approx e^{-(\nu_{1}/4+\nu_{2}/4)t}$ (45)

which is the same form as the averaged localized case Eq. C.2.2, but with
$\nu_{1}\rightarrow\frac{\nu_{1}+\nu_{2}}{2}$.

####  C.3.2 Entangled states

In a highly entangled diffusive state, the phases of each $n$-qubit basis
state would be random, preventing decomposition to single qubits. This means

$\displaystyle\ket{\psi}$
$\displaystyle=\frac{1}{\sqrt{2^{n}}}\Big{(}\ket{\dots
00}+e^{i\phi_{1}}\ket{\dots 01}$
$\displaystyle\qquad\qquad+e^{i\phi_{2}}\ket{\dots 10}+\dots\Big{)}$ (46)
$\displaystyle\rho$
$\displaystyle=\frac{1}{2^{n}}\begin{pmatrix}1&e^{-i\phi_{1}}&e^{-i\phi_{2}}&\dots\\\
e^{i\phi_{1}}&1&e^{i(\phi_{1}-\phi_{2})}\\\
e^{i\phi_{2}}&e^{i(\phi_{2}-\phi_{1})}&1&\\\ \vdots&&&\ddots\end{pmatrix}$
(47)

The fidelity $f=\sum_{i,j}\sigma_{i,j}\rho_{i,j}^{*}$ shows these phases
cancel out at $t=0$ when $\sigma(t=0)=\rho$, giving $f=4^{n}/4^{n}=1$. But how
does an $n$-qubit $\sigma$ evolve?

Similarly to Eq. 40, the equation for a given $\dot{\sigma}_{i,j}$ will have a
relaxation term $-\nu_{1}\sigma_{i,j}$ for each qubit that is 1 in both states
$i$ and $j$, a gain term $\nu_{1}\sigma_{i+2^{k},j+2^{k}}$ for each qubit that
is 0 in both states $i$ and $j$ where $k$ is that qubit’s index in the list of
qubits, and a dephasing term $-\frac{\nu_{1}+\nu_{2}}{2}\sigma_{i,j}$ for each
qubit that is 0 in $i$ or $j$ but 1 in the other. With only the gain term
coupling the matrix elements, and only along the same diagonal, it is
straightforward if time-consuming to solve the set of equations on each
diagonal of $\sigma$. Then the elements of $\sigma\circ\rho^{*}$ can be summed
to calculate the fidelity.

##### Two qubits

It is instructive to work out the fidelity for two qubits. The main diagonal
terms have the simplest contribution, since
$\sum\sigma_{i,i}\rho_{i,i}=\sum\sigma_{i,i}/2^{n}=1/2^{n}$ due to trace
preservation $\Tr(\sigma)=1$. The precise dynamics on the diagonal are
irrelevant for the diffusive fidelity, as was seen even in the unentangled
state.

For off-diagonal terms, many high-index elements of $\sigma$ do not have gain
terms, making them simple exponential decays. For two qubits,

$\displaystyle\dot{\sigma}_{2,3}$
$\displaystyle=-\frac{\nu_{1}+\nu_{2}}{2}\sigma_{2,3}-\nu_{1}\sigma_{2,3}$
$\displaystyle\sigma_{2,3}$
$\displaystyle=\frac{1}{2^{n}}e^{i(\phi_{2}-\phi_{3})}e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}$
(48)

where the indices correspond to binary states, i.e. $2=10_{2},3=11_{2}$. This
matrix element describes qubit 0 in superposition and qubit 1 in state 1,
producing the first and second terms above, respectively. In contrast, a low-
index state depends on decay from high-index states:

$\displaystyle\dot{\sigma}_{0,1}$
$\displaystyle=-\frac{\nu_{1}+\nu_{2}}{2}\sigma_{0,1}+\nu_{1}\sigma_{2,3}$

If we insert the ansatz $\sigma_{0,1}=Ae^{at}+Be^{bt}$ into the equation, then
we find

$\displaystyle\dot{\sigma}_{0,1}$
$\displaystyle=a(Ae^{at}+Be^{bt})+B(b-a)e^{bt}$ $\displaystyle a$
$\displaystyle=-\frac{\nu_{1}+\nu_{2}}{2}$ $\displaystyle b$
$\displaystyle=-\frac{3\nu_{1}+\nu_{2}}{2}$ $\displaystyle B$
$\displaystyle=\frac{\nu_{1}}{2^{n}}e^{i(\phi_{2}-\phi_{3})}\frac{1}{b-a}=-\frac{1}{2^{n}}e^{i(\phi_{2}-\phi_{3})}$

with initial condition

$\displaystyle\sigma_{0,1}(t=0)$
$\displaystyle=A+B=\frac{1}{2^{n}}e^{-i\phi_{1}}$ $\displaystyle A$
$\displaystyle=\frac{1}{2^{n}}\left(e^{-i\phi_{1}}+e^{i(\phi_{2}-\phi_{3})}\right)$

giving the solution

$\displaystyle\sigma_{0,1}=\frac{1}{2^{n}}\Big{[}\left(e^{-i\phi_{1}}+e^{i(\phi_{2}-\phi_{3})}\right)e^{-\frac{\nu_{1}+\nu_{2}}{2}t}-$
$\displaystyle
e^{i(\phi_{2}-\phi_{3})}e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}\Big{]}.$ (49)

In the unentangled limit where all $\phi_{i}\rightarrow 0$, the fidelity
contributions of the two elements $\sigma_{2,3}$ and $\sigma_{0,1}$ partially
cancel, eliminating the faster decay. But with random phase entanglement their
combined contribution generally doesn’t cancel, instead giving

$\displaystyle\sigma_{0,1}\rho_{0,1}^{*}$
$\displaystyle+\sigma_{2,3}\rho_{2,3}^{*}=$
$\displaystyle\frac{1}{(2^{n})^{2}}\Big{[}$
$\displaystyle\left(1+e^{i(\phi_{1}+\phi_{2}-\phi_{3})}\right)e^{-\frac{\nu_{1}+\nu_{2}}{2}t}+$
$\displaystyle\left(1-e^{i(\phi_{1}+\phi_{2}-\phi_{3})}\right)e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}\Big{]}.$
(50)

In the fidelity this is further combined with its complex conjugate, resulting
in cosines of the phases.

So the effect of random phases from diffusive entanglement does not change
which fidelity decay rates are possible, but rather to change the relative
weight of each term. This lesson will generalize to $n$ qubits.

Since Eq. C.1 says fidelity is a sum over matrix elements of
$\sigma\circ\rho^{*}$ (Eq. 36), it is helpful to show this full matrix for two
qubits:

$\displaystyle\sigma\circ\rho^{*}=\frac{1}{4^{2}}\begin{pmatrix}\begin{matrix}4-4e^{-\nu_{1}t}+\\\
e^{-2\nu_{1}t}\end{matrix}&\begin{matrix}\left(1+e^{i\phi_{123}}\right)e^{-\frac{\nu_{1}+\nu_{2}}{2}t}+\\\
\left(-e^{i\phi_{123}}\right)e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}\end{matrix}&\begin{matrix}\left(1+e^{i\phi_{123}}\right)e^{-\frac{\nu_{1}+\nu_{2}}{2}t}+\\\
\left(-e^{i\phi_{123}}\right)e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}\end{matrix}&e^{-\frac{2\nu_{1}+2\nu_{2}}{2}t}\\\
c.c.&2e^{-\nu_{1}t}-e^{-2\nu_{1}t}&e^{-\frac{\nu_{1}+\nu_{2}}{2}t}&e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}\\\
c.c.&c.c.&2e^{-\nu_{1}t}-e^{-2\nu_{1}t}&e^{-\frac{3\nu_{1}+\nu_{2}}{2}t}\\\
c.c.&c.c.&c.c.&e^{-2\nu_{1}t}\end{pmatrix}$ (51)

where $\phi_{123}\equiv(\phi_{1}+\phi_{2}-\phi_{3})$. The diagonal matrix
elements only have terms for relaxation to lower states and gain from higher
states. Off-diagonal terms also have dephasing due to superposition qubits.
The number of superposition qubits for an element can be counted from its
indices, i.e. element $(\sigma\circ\rho^{*})_{2,3}$ measures the superposition
between states $2=10_{2}$ and $3=11_{2}$, so only qubit 0 is in superposition
between $0$ and $1$. This contributes one factor of
$e^{-\frac{\nu_{1}+\nu_{2}}{2}t}$, which combines with the decay term
$e^{-\nu_{1}t}$ from qubit 1. Off-diagonal gain terms are the only terms for
which the random phase factors of $\sigma$ and $\rho^{*}$ do not cancel,
making them the only way phase dependence enters the fidelity.

#####  $n$ qubits

An $n$-qubit average fidelity can be found based on this two-qubit example.
First, when $\nu_{1}=0$ the coupling between elements disappears, eliminating
the phase-dependent terms in fidelity. This is just the unentangled case with
the solution Eq. C.3.1.

When $\nu_{1}\neq 0$ the fidelity $f=\sum_{i,j}\sigma_{i,j}\rho_{i,j}^{*}$
generally depends on the phases. However we can calculate the average fidelity
over random phases $\phi_{i}$. The motivation for this phase averaging is
twofold. First, the $2^{n}$ initial conditions of basis states $\\{\ket{p}\\}$
are averaged over, with each producing statistically independent phases.
Second, even if the phases have some probability distribution with finite
variance and zero mean, the number of independent phases that sum in the
fidelity grows quickly with number of qubits. Two-qubit fidelity has only one
independent phase $\phi_{123}$ per initial condition and time step, but three-
qubit fidelity has four, each appearing with multiple decay rates. As the
number of qubits grows, most decay rate coefficients sum over many independent
phases, producing an additional effective averaging. Between these two
averaging effects over the uniform distribution of phases, the variance of the
mean should quickly shrink towards zero.

Some phase-dependent terms will always be controlled by just a few phases,
specifically $\sim n$ phases rather than one due to $n$-fold symmetry in
permuting qubit identitites. This could break the assumption of phase
averaging. However, in the worst case a single such poorly-averaged decay rate
coefficient will have fidelity contribution $\sim n/4^{n}$, so the effect of
each will be small, as well as independent.

Assuming then that averaging over each phase gives an approximately correct
total fidelity, one can use $\langle e^{i\phi_{i}}\rangle=0$ to eliminate
phase dependence in the fidelity. This does not reduce to the unentangled
case. Off-diagonal elements of $\sigma\circ\rho^{*}$ lose their gain terms but
not their relaxation and dephasing terms. This means probability flow from
high to low states retains the loss from high but not the gain to low.
Equivalently, all coupling terms in the differential equations for $\sigma$
can be dropped, greatly simplifying each matrix element to a single term each.
Diagonal elements are not changed by phase averaging since they have no phase
dependence, so they retain their gain terms.

The general n-qubit average fidelity can be constructed from this
understanding. The fidelity is

$\displaystyle
f=\frac{1}{2^{n}}+\frac{1}{4^{n}}\sum_{k=1}^{n}e^{-k\frac{\nu_{1}+\nu_{2}}{2}t}{n\choose
k}2^{k}\sum_{l=0}^{n-k}e^{-l\nu_{1}t}{n-k\choose l}.$ (52)

This construction goes as follows: First count the diagonal terms as the sum
$\frac{1}{2^{n}}$ due to trace preservation. Then for each off-diagonal term
$\sigma_{i,j}$, count the number of qubits $k$ in superposition between states
$i$ and $j$. For $k$ qubits in superposition, the effects of dephasing combine
for a decay rate contribution of $-k\frac{\nu_{1}+\nu_{2}}{2}$. The number of
qubit permutations with $k$ qubits in superposition is ($n$ choose $k$), and
there are $2^{k}$ ways the superposition qubits can be flipped between matrix
index contribution ordering $0,1$ and $1,0$, i.e. matrix element
$0,3=000_{2},011_{2}$ and $2,1=010_{2},001_{2}$ have the same dephasing rate.
These two factors count the number of matrix elements with that dephasing
rate. The non-superposition qubits can each be in state 0 or 1, but phase
averaging has eliminated gain to 0 while preserving relaxation from 1. This
requires a sum over number of excited qubits. Elements with $l$ excited qubits
have additional exponential fidelity decay at a rate $-l\nu_{1}$, and there
are ($(n-k)$ choose $l$) permutations of these excited qubits. Lastly, since
the diagonal elements must be excluded from these sums for not following the
same rules after averaging, the index $k=0$ for zero superposition qubits must
be removed. It is helpful to subtract this term after the fact.

One validation of Eq. 52 is that it can be easily modified to describe the
unentangled case of Eq. C.3.1. When all phases $\phi_{i}$ are zero, instead of
the gain terms averaging to zero they perfectly cancel the decay terms. To
reflect this, remove the effect of decay in Eq. 52 via the replacement
$e^{-l\nu_{1}t}\rightarrow 1$. The sum over $\l$ then evaluates to $2^{n-k}$
and simplifying recovers Eq. C.3.1.

The full Eq. 52 simplifies to

$\displaystyle f$
$\displaystyle=\frac{1}{4^{n}}\left(1+e^{-\nu_{1}t}+2e^{-\frac{\nu_{1}+\nu_{2}}{2}t}\right)^{n}$
$\displaystyle\quad-\frac{1}{4^{n}}\left(1+e^{-\nu_{1}t}\right)^{n}+\frac{1}{2^{n}}$
(53) $\displaystyle\approx
1-n\left(\nu_{1}/2+\nu_{2}/4\right)t+\frac{1}{2^{n}}n(\nu_{1}/2)t$
$\displaystyle\quad+O(n^{2}(\nu_{1}/2+\nu_{2}/4)^{2}t^{2})$ (54)
$\displaystyle\approx e^{-n\left(\nu_{1}/2+\nu_{2}/4\right)t}.$ (55)

Comparing this averaged diffusive entangled decay to the diffusive unentangled
decay Eq. C.3.1, the effect of $\nu_{1}$ on the initial rate approximately
doubled but the effect of $\nu_{2}$ remains the same. The extra $\nu_{1}$
contribution comes from relaxation terms that no longer cancel with gain
terms.

The diffusive fidelity Eq. C.3.2 can be understood more intuitively by noting
the phase averaging decouples all matrix elements except the $n$-qubit
diagonal. A single decoupled, phase-averaged qubit within the $n$-qubit
system, ignoring the trace preservation from the $n$-qubit diagonal, can be
considered as having a fidelity matrix $\sigma\circ\rho^{*}$ of

$\displaystyle\frac{1}{4}\begin{pmatrix}1&e^{-\frac{\nu_{1}+\nu_{2}}{2}t}\\\
e^{-\frac{\nu_{1}+\nu_{2}}{2}t}&e^{-\nu_{1}t}\end{pmatrix}$ (56)

where the gain term has been removed relative to Eq. 44. The fidelity
$f=\sum_{i,j}\sigma_{i,j}\rho_{i,j}^{*}$ of this to the $n$th power for $n$
qubits forms the first term in Eq. C.3.2. The latter two terms are corrections
for the $n$-qubit diagonal, since the diagonal has no phases to average as is
clear from initial state Eq. 47. The second term removes the erroneous
diagonal term and the third replaces it with the appropriate $1/2^{n}$ from
trace preservation.