#  Benchmarking quantum simulators using ergodic quantum dynamics

Daniel K. Mark  Center for Theoretical Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Joonhee Choi  California Institute of
Technology, Pasadena, CA 91125, USA  Adam L. Shaw  California Institute of
Technology, Pasadena, CA 91125, USA  Manuel Endres  California Institute of
Technology, Pasadena, CA 91125, USA  Soonwon Choi<EMAIL_ADDRESS>Center for
Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA

###### Abstract

We propose and analyze a sample-efficient protocol to estimate the fidelity
between an experimentally prepared state and an ideal target state, applicable
to a wide class of analog quantum simulators without advanced spatiotemporal
control. Our protocol relies on universal fluctuations emerging from generic
Hamiltonian dynamics, that we discover in the present work. It does not
require fine-tuned control over state preparation, quantum evolution, or
readout capability, while achieving near optimal sample complexity: a percent-
level precision is obtained with $\sim 10^{3}$ measurements, independent of
system size. Furthermore, the accuracy of our fidelity estimation improves
exponentially with increasing system size. We numerically demonstrate our
protocol in a variety of quantum simulator platforms including quantum gas
microscopes, trapped ions, and Rydberg atom arrays. We discuss applications of
our method for tasks such as multi-parameter estimation of quantum states and
processes.

Introduction.—Recent advances in quantum technology have opened new ways to
probe quantum many-body physics, leading to the first observations of novel
phases of matter [1, 2, 3, 4, 5], quantum thermalization [6, 7, 8], and non-
equilibrium phenomena [9, 10, 11, 12, 13, 14]. However, in order to advance to
the stage where quantum devices produce highly accurate data, it is important
to quantify the performance of said devices. One method to do so is quantum
device benchmarking [15]—verifying that a device accurately produces a state
$\rho$ close to the desired state $\ket{\Psi}$, in the presence of
imperfections and noise, measured by the fidelity
$F=\bra{\Psi}\rho\ket{\Psi}$. A high fidelity certifies that any property of
the prepared state is close to that of the target state [16], hence is widely
used in theory to quantify the goodness of state preparation. Experimentally
measuring the fidelity is important for building, characterizing, and
improving increasingly complex and precise systems.

Several methods to benchmark quantum devices have been proposed. A naïve
approach is to perform quantum state tomography [16, 17, 18, 19], in which an
experimental state is fully characterized by measurements in many different
bases. This approach, however, is impractical even for relatively small
systems as it requires prohibitively many measurements. Alternatively, recent
proposals pointed out that one can directly estimate the fidelity with a small
number of measurements in randomly chosen bases [20, 21, 22, 23, 24, 25, 26,
27]. These methods rely on implementing highly engineered quantum gates that
satisfy certain statistical properties and are not readily applicable to
quantum devices with limited controllability. Other existing benchmarking
protocols require sophisticated controls and are challenging to implement [28,
29, 30, 31, 32, 33, 34, 35, 36]. In particular, we emphasize that analog
quantum simulators are typically designed to realize specific forms of many-
body Hamiltonians and lack the ability to implement arbitrary unitary
operations. Thus, it remains an outstanding challenge to develop a general
benchmarking method with minimal requirements on hardware capability.

Figure 1:  (a) Schematic of our benchmarking protocol. The fidelity estimator
$F_{d}$ is evaluated from experimental snapshots $\\{z_{1},\dots,z_{M}\\}$ of
a state $\rho(t)$ obtained after quench dynamics, compared against classical
computation of the ideal state $\ket{\Psi(t)}$ in the absence of error (see
Table 1). (b-e) Numerical demonstrations. $F_{d}$ closely tracks the fidelity
decay over evolution time normalized in units of Rabi frequency or tunnelling
strength (black dashed) between noisy and ideal quench dynamics in a wide
class of analog simulators, including 1D Bose-Hubbard, integrable 1D Fermi-
Hubbard, 1D trapped-ion, and 2D Rydberg array models at finite effective
temperature, see SM [37] for details. Previously proposed benchmarks
$F_{\textrm{XEB}}$ [30] [green dotted, out of scale in (b,d)] or $F_{c}$ [36]
(purple dot-dashed) fail to estimate $F$ for these systems.

In this work, we propose and analyze a benchmarking protocol that requires
minimal experimental control: one prepares an initial state, time-evolves it
under a natural Hamiltonian of the system, and performs measurements in a
fixed basis (Fig. 1). We show that, with appropriate data-processing (enabled
by classical computation), this simple experiment gives an estimate for the
fidelity $F$—encapsulating the combined effects of errors in state
preparation, quench evolution, and readout 111More accurately, the fidelity is
a quantity defined between two quantum states and does not include readout
errors. Our benchmark estimates the state fidelity when readout errors are
negligible.—with a small number of measurements. Most importantly, our method
works for generic quench dynamics far from fine-tuned cases, including: at
finite effective temperatures, in the presence of symmetries, and in non-qubit
based systems such as itinerant particles on optical lattices, making it
suitable for a wide class of existing platforms.

The key behind our approach is our discovery of universal statistical
fluctuations in the measurement outcome distributions $p(z)$ that arise from
generic quantum dynamics (Fig. 2). Previously, such universal fluctuations in
$p(z)$ were only known to occur in ideal, controlled dynamics such as random
unitary circuits (RUCs), where $\\{p(z)\\}$ approximately follows the Porter-
Thomas distribution [39, 31]. Leveraging our discovery and classical
computation, we design a novel statistic: a real number $f(z)$ associated to
every measurement outcome $z$ such that its average over experimentally
obtained samples, $\hat{F}_{d}\equiv\langle f(z)\rangle_{\textrm{exp}}$,
converges quickly to the many-body fidelity $F$. In other words, $\hat{F}_{d}$
is a computationally-assisted, efficiently measurable observable [40, 41] that
estimates the fidelity.

The ability to estimate fidelity serves as a foundation for two tasks: (i)
_target state benchmarking_ , where the overlap between an experimentally
prepared state and a pure target state is measured via a high-fidelity quench
time evolution, and (ii) _quantum process benchmarking_ , in which the
fidelity decay of quench dynamics is monitored over the course of evolution.

Protocol.—We focus on describing and numerically demonstrating our protocol,
before returning to why it works. Our benchmarking method consists of three
steps: experiment, computation, and data processing [Table 1 and Fig. 1(a)].
The initial state of our protocol can either be an easy-to-prepare state or a
more complex state that one wishes to benchmark. After quench evolution for a
fixed time $t$, the experimental state $\rho(t)$ is measured in any fixed
basis $\\{\ket{z}\\}$. Convenient choices of $\\{\ket{z}\\}$ include the set
of bitstrings in two-level (qubit) systems or real-space particle number
configurations in quantum gas microscopes. Repeating the state preparation and
measurement $M$ times, one obtains measured configurations $\\{z_{1},$ $\dots$
$,z_{M}\\}$, with each $z_{i}$ sampled from the distribution $q(z,t)$ $\equiv$
$\bra{z}\rho(t)\ket{z}$. Our protocol estimates the fidelity by using a small
number of samples to compare the empirical distribution $q(z,t)$ against a
theoretical, target distribution $p(z,t)$. By classical computation, we obtain
$p(z,t)$ and its infinite-time average
$p_{\text{avg}}(z)\equiv\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}p(z,t)dt$.
In practice, one may average over a finite duration $T$ as an approximation,
at the expense of slightly larger statistical errors. Then, we evaluate the
rescaled outcome probabilities $\tilde{p}(z,t)\equiv p(z,t)/p_{\text{avg}}(z)$
and the normalization factor
$\mathcal{Z}(t)\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{2}$.

The classical computation determines our statistic $f(z)$, while the
experimental samples determine which outcomes $z$’s to use when evaluating the
statistic. More specifically, we estimate the fidelity with the empirical
average

$\hat{F}_{d}(t)=\langle
f(z)\rangle_{\textrm{exp}}=\frac{1}{M}\sum_{i=1}^{M}2\tilde{p}(z_{i},t)/\mathcal{Z}(t)-1\leavevmode\nobreak\
.$ (1)

This explicitly defines the statistic $f(z)$, that also depends on
$\ket{\Psi_{0}}$, $H$, and $t$. In the limit $M\\!\rightarrow\\!\infty$, this
converges to our benchmark
$F_{d}(t)=2\Big{[}\sum_{z}q(z,t)\tilde{p}(z,t)\Big{]}/\mathcal{Z}(t)-1$ . This
benchmark can be understood as a weighted covariance between the empirical and
ideal distributions. We show that $F_{d}(t)$ approximates the fidelity $F$ for
a wide class of quantum systems, both for uncorrelated, infrequent incoherent
errors, or local and weak global coherent errors [42, 43, 44, 37], and
rigorously prove our statement for isolated single errors and long evolution
times.

Table 1: Proposed benchmarking protocol. Experiment:
---
|  1\. Prepare an initial state $\rho_{0}$, which approximates a pure state
$\outerproduct{\Psi_{0}}{\Psi_{0}}$.
|  2\. Evolve the system under its natural Hamiltonian $H$ for a time $t$.
|  3\. Measure the evolved state $\rho(t)$ in a natural basis, obtaining
configurations $\\{z_{j}\\}_{j=1}^{M}$.
Computation: Classically compute
|  1\.
$p(z,t)\equiv\absolutevalue{\bra{z}\ket{\Psi(t)}}^{2}=\absolutevalue{\bra{z}\exp(-iHt)\ket{\Psi_{0}}}^{2}$,
|  2\.
$p_{\text{avg}}(z)\equiv\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}p(z,t)\leavevmode\nobreak\
dt$,
|  $\leavevmode\nobreak\ \tilde{p}(z,t)\equiv p(z,t)/p_{\text{avg}}(z)$.
|  3\. $\mathcal{Z}(t)\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{2}$.
Data processing: Evaluate
$\hat{F}_{d}(t)\equiv\frac{2}{M}\Big{[}\sum_{i=1}^{M}\tilde{p}(z_{i},t)\Big{]}/\mathcal{Z}(t)-1\approx
F_{d}(t)$ ,
which approximates the fidelity $F=\bra{\Psi(t)}\rho(t)\ket{\Psi(t)}$.

Figure 1(b-e) numerically demonstrates the use of our estimator for process
benchmarking: tracking the decay of fidelity over time in four different
quantum simulation platforms. For each platform, we simulate an initial
product state undergoing natural Hamiltonian dynamics in the presence of
experimentally relevant errors [37]. We confirm that $F_{d}$ successfully
traces the fidelity decay in regimes where previously proposed fidelity
estimators $F_{\textrm{XEB}}$ [30] and $F_{c}$ [36] do not. This is because
$F_{c}$ and $F_{\textrm{XEB}}$ (reviewed in the SM) assume that $p(z,t)$
satisfy statistical properties (discussed below) that are in general not
satisfied by natural Hamiltonian dynamics, e.g. $F_{c}$ requires the system to
evolve at infinite effective temperature. We now turn to the underlying
principles of our protocol: emergent universal statistics, speckle based
benchmarking, and measurement-basis independence.

Figure 2: Emergent universal statistics. (a) During many-body Hamiltonian
evolution of a pure state, the probability $p(z,t)$ of measuring an outcome
$z$ fluctuates around its average value $p_{\text{avg}}(z)$ (dashed, three
distinct $z$’s highlighted). Here, we consider a Bose-Hubbard model and
present $p(z,t)$ rescaled by the Hilbert space dimension $D$. (b) The
histogram of $p(z,t)$ over all $z$ at a fixed $t$ [red dots in (a)], is non-
universal (red dashed), far from the Porter-Thomas (PT) distribution (black
dashed). (c) Rescaling each $p(z,t)$ by $p_{\text{avg}}(z)$ yields
$\tilde{p}(z,t)$, whose histogram follows the PT distribution. (d,e,f) The
histograms of $\tilde{p}(z,t)$ (blue bars) follow the universal PT
distribution in all models considered in this work, whereas those of the bare
$p(z,t)$ rescaled by $D$ are non-universal (red dashed).

Emergent universal statistics.—The statistical properties of $p(z)$ have been
extensively studied in deep random unitary circuits (RUCs). When the output
state of a typical deep RUC is measured, $p(z)$ is not perfectly uniform, but
exhibits a speckle pattern: over different $z$’s, $p(z)$ fluctuates about
$1/D$ due to random interference in coherent quantum dynamics, with $D$ the
Hilbert space dimension. While the details of the fluctuations—which $p(z)$’s
are larger—sensitively depends on the particular choice of RUC, the
statistical properties of $p(z)$ are universal. Specifically, the fraction of
$p(z)$’s in a given interval $p(z)\in[x,x+dx]$ is given by the Porter-Thomas
(PT) distribution: $P[p(z)=x]dx=\mu^{-1}\exp(-x/\mu)dx$ with mean $\mu=1/D$.
This enables the existing benchmarks $F_{\textrm{XEB}}$ and $F_{c}$.
Specifically, they utilize the fact that the PT distribution has a second
moment equal to two [30, 31, 36]. Previously, it was unclear under what
conditions the PT distribution can arise, other than from RUCs and fine-tuned
Hamiltonian dynamics.

In fact, for generic time-independent Hamiltonian dynamics, the raw
distribution $p(z)$ does not follow the PT distribution [Fig. 2(a,b)]. This is
due to the presence of energy conservation or symmetries, which causes
systematic trends in $p(z)$ and distorts its distribution away from PT. For
example, in any state with positive effective temperature, low-energy
configurations are measured more frequently than high-energy ones [37]. While
previous work discovered PT distributions in certain local observables [36],
this work concerns global observables under general conditions such as finite
effective temperature.

Our key insight is that the systematic trends in $p(z)$ can be removed simply
by rescaling $p(z)$ with its time-averaged value, leaving only random relative
fluctuations $\tilde{p}(z)\equiv p(z)/p_{\text{avg}}(z)$ that follow the PT
distribution with mean $\mu=1$ [Fig. 2(c-f)].

######  Theorem.

(informal) Consider an initial state $\ket{\Psi_{0}}$, evolved for time $t$
under a Hamiltonian $H$ satisfying the $k$-th no-resonance condition for a
large integer $k$, and measured in a complete basis $\\{\ket{z}\\}$. For
sufficiently late times $t$, the rescaled probabilities $\tilde{p}(z,t)$
follow the Porter-Thomas distribution, up to a correction bounded by the
inverse effective Hilbert space dimension
$D_{\beta}^{-1}\equiv\sum_{z,E}\absolutevalue{\innerproduct{z}{E}}^{4}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{4}/p_{\text{avg}}(z)$,
where $\\{\ket{E}\\}$ are the eigenstates of $H$.

Rigorous statements of our Theorem and their proofs are presented in the SM
[37]. The only assumption of this Theorem is the $k$-th no-resonance
condition, stating that $\sum_{i=1}^{k}E_{\alpha_{i}}$ $=$
$\sum_{i=1}^{k}E_{\beta_{i}}$ if and only if the $k$ indices $(\alpha_{i})$
are a permutation of $(\beta_{i})$ [45, 46, 47, 48, 49]. That is, the
eigenvalues $\\{E_{j}\\}$ of $H$ possess no resonant structures. This
condition is expected to hold for generic ergodic Hamiltonians [45, 46], and
we find that it even holds in some integrable systems such as the 1D Fermi-
Hubbard model [Fig. 1(c)] [37, 50]. The effective dimension $D_{\beta}$
quantifies the size of the Hilbert space explored during quench evolution
(that can be probed in the $\\{\ket{z}\\}$ basis). $D_{\beta}$ is similar to a
participation ratio of $\ket{\Psi_{0}}$ and $\ket{z}$, when they are
decomposed in the energy eigenbasis, which generically increases exponentially
with increasing system size, leading to a better agreement with the PT
distribution and enabling our protocol to be increasingly accurate.

Our Theorem states that the outcome distribution factorizes into two parts
$p(z,t)=p_{\text{avg}}(z)\times\tilde{p}(z,t)$: systematic values
$p_{\text{avg}}(z)$ and random Porter-Thomas fluctuations $\tilde{p}(z,t)$.
The systematic value is related to thermalization and does not distinguish
between pure and mixed states, usually set by coarse-grained information such
as the total energy $\bra{\Psi_{0}}H\ket{\Psi_{0}}$. Meanwhile, the
fluctuations originate from random interference and average away in a mixed
state [37]. These fluctuations are highly sensitive to details of the initial
state and evolution, serving as a “fingerprint” that enables their
benchmarking.

Speckle-based benchmarking.—We provide an intuitive explanation of our
benchmark $F_{d}$, based on two properties: (i) the second moment of the
rescaled $\tilde{p}(z)$ is
$\mathcal{Z}=\sum_{z}p_{\text{avg}}(z)\tilde{p}(z)^{2}\approx 2$ for $p(z)$
arising from ideal unitary evolution, (ii) the experimental distribution
$q(z)$ in the presence of errors can be expressed as a linear combination
$q(z)=Fp(z)+(1-F)p_{\perp}(z)$, where $p_{\perp}(z)$ is uncorrelated with the
ideal distribution $p(z)$ in the following sense:
$\mathbb{E}_{z}[\tilde{p}(z)\tilde{p}_{\perp}(z)]\approx\mathbb{E}_{z}[\tilde{p}(z)]\mathbb{E}_{z}[\tilde{p}_{\perp}(z)]=1$,
with $\mathbb{E}_{z}[\cdot]\equiv\sum_{z}p_{\text{avg}}(z)[\cdot]$ and
$\tilde{p}_{\perp}\equiv p_{\perp}/p_{\text{avg}}$. The second property relies
on an assumption that the speckle pattern in $\tilde{p}$ significantly changes
under any error, which has been rigorously proven for RUCs [43, 42].

Using these properties, our estimator $F_{d}$ is designed to isolate the
desired “fingerprint,” taking value 1 when $q(z)=p(z)$ and 0 when
$q(z)=p_{\perp}(z)$, which in turn implies $F_{d}\approx F$ [37]. We emphasize
that it is essential to use the rescaled $\tilde{p}(z)$ to estimate the
fidelity; otherwise $p(z)$ and $p_{\perp}(z)$ exhibit large correlations due
to their shared systematic values $p_{\text{avg}}(z)$.

In fact, the second condition can be relaxed. Under local coherent or
incoherent errors, the relation $F_{d}\approx F$ can be shown at late times
without any assumption on $q(z)$, solely based on the Eigenstate
Thermalization Hypothesis and no-resonance conditions [37]. We also argue that
this result extends to multiple stochastic errors, and small coherent errors
in the quench Hamiltonian, as verified by various numerical simulations [37].
Furthermore, we also verify that $F_{d}\approx F$ holds even for relatively
short quenches, well before the PT distribution emerges in $\tilde{p}(z)$,
owing to the time-dependent adjustment factor $\mathcal{Z}(t)$ in $F_{d}$.
This factor is inspired by $F_{c}$ [36] and its effect is illustrated in Fig.
1(b-e) and the SM.

Measurement-basis independence.—A surprising feature of our approach is that
the fidelity is estimated from measurements in a fixed basis, despite the fact
that the fidelity also depends on phase information not accessed from such
measurements. Nevertheless, our protocol works because quench evolution
transforms the effects of physically relevant errors, including phase errors,
into a form detectable by generic local measurements, after a short delay time
[Fig. 3(a)].

In our examples above, the quench dynamics plays two roles simultaneously: it
enables our protocol, but also generates imperfect quantum evolution (due to
errors), whose fidelity decay is measured. If the quench evolution were
perfect, the measured fidelity would only reflect the state preparation error
of the (potentially interesting) target initial state. We present numerical
demonstrations of such target state benchmarking in the SM [37].

Figure 3:  Performance analysis. (a) Numerical simulation of open system
dynamics of a 1D Bose-Hubbard model with $N=9$ particles on $N$ sites
($D=24310$). When a single error occurs (black arrow), $F_{d}$ (solid line,
green) closely approximates the fidelity $F$ (dashed line, green) after a
short delay time $\tau_{r}$. Averaging $F_{d}$ and $F$ over (potentially many)
stochastic errors at different times gives their values for the mixed state
$\rho$ (blue). The inset shows uncertainties and errors in our method. In
particular, $F_{d}$ slightly deviates from $F$, quantified by the difference
$\delta_{\textrm{sys}}$ between $F$ and the time-averaged $F_{d}$ (red arrow)
and by the fluctuations $\delta_{\textrm{temp}}$ of $F_{d}$ over time (green
arrow). Furthermore, a finite number of samples $M$ results in a statistical
uncertainty $\delta_{\textrm{stat}}$ in the unbiased estimator $\hat{F}_{d}$
(orange error bar). (b,c) Both $\delta_{\textrm{sys}}$ and
$\delta_{\textrm{temp}}$ decrease exponentially with system size,
quantitatively agreeing with analytic predictions (dashed lines). (d) The
sample complexity $M\delta^{2}_{\textrm{stat}}$ increases weakly with $N$ at
fixed, early times (dotted line). At late times, it approaches the
$N$-independent value $1+2F-F^{2}$ (dashed line). Error bars in (b-d) indicate
variations over an ensemble of disordered Hamiltonians. See SM [37] for
details.

Performance analysis.—Our Theorem enables us to predict how well $F_{d}$
estimates the fidelity. First, we point out that it suffices to study the
effect of a single error when the error rate is sufficiently small. This is
because incoherent noisy dynamics can be “unravelled” into an ensemble of
stochastic pure state trajectories [51] each corresponding to a fixed
occurrence of errors [Fig. 3(a)]. As long as they are sufficiently infrequent,
the effect of multiple errors can be understood from that of a single error
[42, 43].

We showcase the performance of $F_{d}$ under realistic conditions by
numerically simulating the 1D Bose-Hubbard model. See Figure 3(a). We quantify
the performance of $F_{d}$ along several axes [Fig. 3(a), inset]: the
systematic error $\delta_{\textrm{sys}}$ refers to the difference between the
true fidelity $F$ and the time-averaged $F_{d}$, while the
$\delta_{\textrm{temp}}$ quantifies how $F_{d}$ fluctuates over time. Finally,
the statistical fluctuation (or, sample complexity) $\delta_{\textrm{stat}}$
measures the uncertainty of the estimated $\hat{F}_{d}$ associated with a
finite number of samples $M$, and hence the number of samples required to
determine $F_{d}$ up to a desired precision. See SM [37] for further details.

Our Theorem allows analytical estimation of these quantities in the limit of
long evolution: $\delta_{\textrm{sys}}$ and $\delta_{\textrm{temp}}$ are
respectively $O(D_{\beta}^{-1})$ and $O(D_{\beta}^{-1/2})$ [37]. Hence, both
the accuracy and precision of our benchmark improve exponentially with
increasing system size, explicitly confirmed in our numerical simulations
[Fig. 3(b,c)]. Meanwhile, the sample complexity has optimal scaling. It is
system size independent for long evolutions:
$M\delta^{2}_{\textrm{stat}}\approx 1+2F-F^{2}$ [Fig. 3(d), dashed line]. For
short quench evolution, the sample complexity grows weakly with system size
[Fig. 3(d), dotted line]. Finally, in the presence of incoherent errors, the
finite response time $\tau_{r}$ leads to a slight delay between $F_{d}$ and
the continuously-decaying fidelity [42, 52, 37].

Figure 4:  State and process benchmarking with $F_{d}$. We numerically
simulate the estimation of: (a) the phase $\phi$ of a GHZ-like initial state
in a 2D Rydberg system; (b) the ratio of the interaction and tunneling
strengths $U/\Omega$ in a 1D Bose-Hubbard model; and (c) ten disordered on-
site potentials in a trapped ion model. Parameters are estimated by maximizing
$\hat{F}_{d}$ over simulated parameter values, using measurements after error-
free (blue) or noisy (red) quench evolution. The error bars and shaded regions
indicate the statistical uncertainties in $\hat{F}_{d}$ and the parameter
values with 1000 samples. See SM [37] for details. (a,b) Both $F$ (black
lines) and $\hat{F}_{d}$ (markers) are consistent and simultaneously maximized
at the true parameter value (dotted line). (c) Reconstructed disorder
potential values (markers) are consistent with their true values (lines).

Limitations.—While our protocol is applicable to generic quantum many-body
systems, it may fail in special cases in which our Theorem is not applicable.
Examples include systems with weakly- or non-ergodic dynamics, or the presence
of correlated nonlocal errors. We provide detailed analysis and potential
resolutions of known failure cases in the SM [37].

Applications.— The ability to measure the fidelity enables further
applications. As examples, we show that one can simultaneously estimate
multiple parameters of quantum states or Hamiltonians. The key observation is
that, given the ability to measure the fidelity between a theoretical model
and experiment, one can vary classical simulation model parameters to maximize
the estimated fidelity [31, 36]. This optimization requires classical
computation but no further data acquisition. We numerically demonstrate this
idea in three different examples and verify that the extracted parameter
values are close to the actual ones, even in the presence of noise. See Figure
4.

Acknowledgments.— We would like to thank Wen Wei Ho, Anant Kale, Eun Jong Kim,
Sherry Zhang, and Peter Zoller for useful discussions. In particular, we thank
Sherry Zhang for a careful reading of this manuscript and Anant Kale for
insightful discussion at the early stage of the project. We acknowledge
funding provided by the NSF Physics Frontiers Centers IQIM (NSF PHY-1733907)
and CUA (NSF PHY-1734011), the AFOSR YIP (FA9550-19-1-0044), the DARPA ONISQ
program (W911NF2010021), the Army Research Office MURI program
(W911NF2010136), the NSF QLCI program (2016245), the DOE (DE-SC0021951), NSF
CAREER award 1753386 and NSF CAREER award 2237244. Support is also
acknowledged from the U.S. Department of Energy, Office of Science, National
Quantum Information Science Research Centers, Quantum Systems Accelerator. JC
acknowledges support from the IQIM postdoctoral fellowship. ALS acknowledges
support from the Eddleman Quantum graduate fellowship.

## References

  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature (London) 415, 39 (2002).
  * Aidelsburger  _et al._ [2013] M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Realization of the Hofstadter Hamiltonian with ultracold atoms in optical lattices, Phys. Rev. Lett. 111, 185301 (2013).
  * Jotzu  _et al._ [2014] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Experimental realization of the topological Haldane model with ultracold fermions, Nature (London) 515, 237 (2014).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, _et al._ , Quantum phases of matter on a 256-atom programmable quantum simulator, Nature (London) 595, 227 (2021).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).
  * Tang  _et al._ [2018] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Thermalization near integrability in a dipolar quantum Newton’s cradle, Phys. Rev. X 8, 021030 (2018).
  * Chen  _et al._ [2021] F. Chen, Z.-H. Sun, M. Gong, Q. Zhu, Y.-R. Zhang, Y. Wu, Y. Ye, C. Zha, S. Li, S. Guo, H. Qian, H.-L. Huang, J. Yu, H. Deng, H. Rong, J. Lin, Y. Xu, L. Sun, C. Guo, N. Li, F. Liang, C.-Z. Peng, H. Fan, X. Zhu, and J.-W. Pan, Observation of strong and weak thermalization in a superconducting quantum processor, Phys. Rev. Lett. 127, 020602 (2021).
  * Neyenhuis  _et al._ [2017] B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C. Lee, P. Richerme, Z.-X. Gong, A. V. Gorshkov, and C. Monroe, Observation of prethermalization in long-range interacting spin chains, Science Adv. 3, e1700672 (2017).
  * Zhang  _et al._ [2017] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator, Nature (London) 551, 601 (2017).
  * Choi  _et al._ [2019] J. Choi, H. Zhou, S. Choi, R. Landig, W. W. Ho, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, D. A. Abanin, and M. D. Lukin, Probing quantum thermalization of a disordered dipolar spin ensemble with discrete time-crystalline order, Phys. Rev. Lett. 122, 043603 (2019).
  * Rubio-Abadal  _et al._ [2020] A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Floquet prethermalization in a Bose-Hubbard system, Phys. Rev. X 10, 021044 (2020).
  * Peng  _et al._ [2021] P. Peng, C. Yin, X. Huang, C. Ramanathan, and P. Cappellaro, Floquet prethermalization in dipolar spin chains, Nature Physics 17, 444 (2021).
  * Mi  _et al._ [2021a] X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, _et al._ , Time-crystalline eigenstate order on a quantum processor, Nature (London) 601, 531–536 (2021a).
  * Eisert  _et al._ [2020] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi, Quantum certification and benchmarking, Nature Rev. Phys. 2, 382 (2020).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Cramer  _et al._ [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Efficient quantum state tomography, Nature Comm. 1, 149 (2010).
  * Gross  _et al._ [2010] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, Quantum state tomography via compressed sensing, Phys. Rev. Lett. 105, 150401 (2010).
  * Christandl and Renner [2012] M. Christandl and R. Renner, Reliable quantum state tomography, Phys. Rev. Lett. 109, 120403 (2012).
  * Flammia and Liu [2011] S. T. Flammia and Y.-K. Liu, Direct fidelity estimation from few Pauli measurements, Phys. Rev. Lett. 106, 230501 (2011).
  * da Silva  _et al._ [2011] M. P. da Silva, O. Landon-Cardinal, and D. Poulin, Practical characterization of quantum devices without tomography, Phys. Rev. Lett. 107, 210404 (2011).
  * Aaronson [2019] S. Aaronson, Shadow tomography of quantum states, SIAM Journal on Computing 49, STOC18 (2019).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Elben  _et al._ [2019] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states, Phys. Rev. A 99, 052323 (2019).
  * Elben  _et al._ [2020] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Cross-platform verification of intermediate scale quantum devices, Phys. Rev. Lett. 124, 010504 (2020).
  * Liu  _et al._ [2021] Y. Liu, M. Otten, R. Bassirianjahromi, L. Jiang, and B. Fefferman, Benchmarking near-term quantum computers via random circuit sampling, arXiv:2105.05232 (2021).
  * Elben  _et al._ [2022] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox, Nature Reviews Physics , 1 (2022).
  * Ohliger  _et al._ [2013] M. Ohliger, V. Nesme, and J. Eisert, Efficient and feasible state tomography of quantum many-body systems, N. J. Phys. 15, 015024 (2013).
  * Brydges  _et al._ [2019] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Rényi entanglement entropy via randomized measurements, Science 364, 260 (2019).
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Characterizing quantum supremacy in near-term devices, Nature Physics 14, 595 (2018).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandão, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature (London) 574, 505 (2019).
  * Li  _et al._ [2020] Z. Li, L. Zou, and T. H. Hsieh, Hamiltonian tomography via quantum quench, Phys. Rev. Lett. 124, 160502 (2020).
  * Kokail  _et al._ [2021] C. Kokail, R. van Bijnen, A. Elben, B. Vermersch, and P. Zoller, Entanglement hamiltonian tomography in quantum simulation, Nature Physics 17, 936 (2021).
  * Gluza and Eisert [2021] M. Gluza and J. Eisert, Recovering quantum correlations in optical lattices from interaction quenches, Phys. Rev. Lett. 127, 090503 (2021).
  * Hu  _et al._ [2021] H.-Y. Hu, S. Choi, and Y.-Z. You, Classical shadow tomography with locally scrambled quantum dynamics, arXiv:2107.04817 (2021).
  * Choi  _et al._ [2023] J. Choi, A. L. Shaw, I. S. Madjarov, X. Xie, R. Finkelstein, J. P. Covey, J. S. Cotler, D. K. Mark, H.-Y. Huang, A. Kale, _et al._ , Preparing random states and benchmarking with many-body quantum chaos, Nature (London) 613, 468 (2023).
  * [37] See Supplemental material online for the proof of our theorem, detailed analysis of the performance, numerical demonstration, and limitations of our protocol, which includes Refs. [53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67].
  * Note [1] More accurately, the fidelity is a quantity defined between two quantum states and does not include readout errors. Our benchmark estimates the state fidelity when readout errors are negligible.
  * Porter and Thomas [1956] C. E. Porter and R. G. Thomas, Fluctuations of nuclear reaction widths, Phys. Rev. 104, 483 (1956).
  * Garratt  _et al._ [2022] S. J. Garratt, Z. Weinstein, and E. Altman, Measurements conspire nonlocally to restructure critical quantum states, arXiv:2207.09476 (2022).
  * Lee  _et al._ [2022] J. Y. Lee, W. Ji, Z. Bi, and M. P. A. Fisher, Decoding measurement-prepared quantum phases and transitions: from ising model to gauge theory, and beyond, arXiv:2208.11699 (2022).
  * Gao  _et al._ [2021] X. Gao, M. Kalinowski, C.-N. Chou, M. D. Lukin, B. Barak, and S. Choi, Limitations of linear cross-entropy as a measure for quantum advantage, arXiv:2112.01657 (2021).
  * Dalzell  _et al._ [2021] A. M. Dalzell, N. Hunter-Jones, and F. G. S. L. Brandão, Random quantum circuits transform local noise into global white noise, arXiv:2111.14907 (2021).
  * Noh  _et al._ [2020] K. Noh, L. Jiang, and B. Fefferman, Efficient classical simulation of noisy random quantum circuits in one dimension, Quantum 4, 318 (2020).
  * Goldstein  _et al._ [2006] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghi, On the distribution of the wave function for systems in thermal equilibrium, J. Stat. Phys. 125, 1193 (2006).
  * Reimann [2008] P. Reimann, Foundation of statistical mechanics under experimentally realistic conditions, Phys. Rev. Lett. 101, 190403 (2008).
  * Linden  _et al._ [2009] N. Linden, S. Popescu, A. J. Short, and A. Winter, Quantum mechanical evolution towards thermal equilibrium, Phys. Rev. E 79, 061103 (2009).
  * Kaneko  _et al._ [2020] K. Kaneko, E. Iyoda, and T. Sagawa, Characterizing complexity of many-body quantum dynamics by higher-order eigenstate thermalization, Phys. Rev. A 101, 042126 (2020).
  * Huang [2021] Y. Huang, Extensive entropy from unitary evolution, arXiv:2104.02053 (2021).
  * Essler  _et al._ [2005] F. H. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The one-dimensional Hubbard model_ (Cambridge University Press, 2005).
  * Daley [2014] A. J. Daley, Quantum trajectories and open many-body quantum systems, Adv. Phys. 63, 77 (2014).
  * Mi  _et al._ [2021b] X. Mi, P. Roushan, C. Quintana, S. Mandrà, J. Marshall, C. Neill, F. Arute, K. Arya, J. Atalaya, R. Babbush, _et al._ , Information scrambling in quantum circuits, Science 374, 1479 (2021b).
  * Short [2011] A. J. Short, Equilibration of quantum systems and subsystems, N. J. Phys. 13, 053009 (2011).
  * Cotler  _et al._ [2023] J. S. Cotler, D. K. Mark, H.-Y. Huang, F. Hernández, J. Choi, A. L. Shaw, M. Endres, and S. Choi, Emergent quantum state designs from individual many-body wave functions, PRX Quantum 4, 010311 (2023).
  * Sarkar  _et al._ [2014] S. Sarkar, S. Langer, J. Schachenmayer, and A. J. Daley, Light scattering and dissipative dynamics of many fermionic atoms in an optical lattice, Phys. Rev. A 90, 023618 (2014).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Quantum scarred eigenstates in a rydberg atom chain: Entanglement, breakdown of thermalization, and stability to perturbations, Phys. Rev. B 98, 155134 (2018).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Probing many-body dynamics on a 51-atom quantum simulator, Nature (London) 551, 579 (2017).
  * Vasseur and Moore [2016] R. Vasseur and J. E. Moore, Nonequilibrium quantum dynamics and transport: from integrability to many-body localization, J. Stat. Mech. 2016, 064010 (2016).
  * Dankert  _et al._ [2009] C. Dankert, R. Cleve, J. Emerson, and E. Livine, Exact and approximate unitary 2-designs and their application to fidelity estimation, Phys. Rev. A 80, 012304 (2009).
  * Deutsch [2018] J. M. Deutsch, Eigenstate thermalization hypothesis, Rep. Prog. Phys. 81, 082001 (2018).
  * Collins and Nechita [2016] B. Collins and I. Nechita, Random matrix techniques in quantum information theory, Journal of Mathematical Physics 57, 015215 (2016).
  * Nahum  _et al._ [2018] A. Nahum, S. Vijay, and J. Haah, Operator spreading in random unitary circuits, Phys. Rev. X 8, 021014 (2018).
  * Khemani  _et al._ [2018] V. Khemani, A. Vishwanath, and D. A. Huse, Operator spreading and the emergence of dissipative hydrodynamics under unitary evolution with conservation laws, Phys. Rev. X 8, 031057 (2018).
  * Jaksch  _et al._ [1998] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Cold bosonic atoms in optical lattices, Phys. Rev. Lett. 81, 3108 (1998).
  * Pichler  _et al._ [2010] H. Pichler, A. J. Daley, and P. Zoller, Nonequilibrium dynamics of bosonic atoms in optical lattices: Decoherence of many-body states due to spontaneous emission, Phys. Rev. A 82, 063605 (2010).
  * Bluvstein  _et al._ [2021] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, _et al._ , Controlling quantum many-body dynamics in driven rydberg atom arrays, Science 371, 1355 (2021).
  * Kühner and Monien [1998] T. D. Kühner and H. Monien, Phases of the one-dimensional Bose-Hubbard model, Phys. Rev. B 58, R14741 (1998).
  * Gross  _et al._ [2007] D. Gross, K. Audenaert, and J. Eisert, Evenly distributed unitaries: On the structure of unitary designs, Journal of mathematical physics 48, 052104 (2007).

Supplemental Material for:
Benchmarking quantum simulators using ergodic quantum dynamics

###### Contents

  1. A Previous fidelity benchmarks
  2. B Emergence of universal statistical properties from Hamiltonian dynamics
    1. B.1 $k$-th no-resonance condition
    2. B.2 $k$-th no resonance with degeneracies
    3. B.3 $k$-th Hamiltonian twirling identity
    4. B.4 Porter-Thomas distribution over time
    5. B.5 Porter-Thomas distribution over outcomes
  3. C Detailed performance analysis for long quench dynamics
    1. C.1 Ansatz for how speckle patterns change under errors
    2. C.2 Systematic error
      1. C.2.1 Argument based on Eigenstate Thermalization Hypothesis
      2. C.2.2 Extensions beyond local errors
    3. C.3 Variance over time
    4. C.4 Sample complexity
    5. C.5 Moderate error rate: Effects of delay time $\tau_{d}$
    6. C.6 Large error rate: Critical error rate and operator scrambling time
    7. C.7 Symmetry Resolution
    8. C.8 Alternative fidelity estimator $F_{e}$
  4. D Detailed performance analysis for short quench dynamics
    1. D.1 Overview: random matrix product states
    2. D.2 Remarks
    3. D.3 Derivation of result
  5. E Details of numerical simulations
    1. E.1 Bose-Hubbard model
      1. E.1.1 Disorder Averaging
    2. E.2 Fermi-Hubbard model
    3. E.3 Trapped ion model
    4. E.4 Rydberg blockaded model
  6. F Limitations of benchmarking protocol: Failure cases
    1. F.1 Violating the no-resonance condition: Particle-hole symmetry
      1. F.1.1 Real Porter-Thomas distribution
    2. F.2 Violating the no-resonance condition: Non-interacting models
    3. F.3 Small effective dimension: Low temperature quench
    4. F.4 Slow equilibration: Quantum many-body scar states
    5. F.5 Slow equilibration: Inhomogeneous initial state
    6. F.6 Non-local error on low entanglement states
    7. F.7 Non-ergodicity due to symmetry and readout in symmetry-resolved basis
  7. G Greedy algorithm for multiparameter estimation
    1. G.1 Algorithm
  8. H Target state benchmarking
    1. H.1 Coherent local errors
    2. H.2 Comparing non-locally deformed states

##  Appendix A Previous fidelity benchmarks

The cross-entropy benchmark was introduced in Ref. [30] and was prominently
used in Ref. [31] to demonstrate quantum computational advantage. It is
defined as

$F_{\text{XEB}}\equiv D\sum_{z}p(z)q(z)-1\leavevmode\nobreak\ ,$ (S1)

where as in the main text, $p(z)$ and $q(z)$ are the probabilities of
obtaining the outcome $z$ from the ideal simulated state and experimental
states respectively, and $D$ is the Hilbert space dimension, here $D=2^{N}$
where $N$ is the number of qubits. $F_{\text{XEB}}$ has been argued to
estimate the fidelity in deep random unitary circuits (RUCs) [30, 31, 43, 42].

The estimator $F_{c}$ was introduced in Ref. [36], defined as

$F_{c}\equiv 2\frac{\sum_{z}p(z)q(z)}{\sum_{z}p(z)^{2}}-1\leavevmode\nobreak\
.$ (S2)

$F_{c}$ is a slight modification of $F_{\text{XEB}}$ to reliably estimate the
fidelity in shallow circuits. This is because of the onset of “local
scrambling” even in shallow circuits, discussed in Ref. [54]. In Ref. [36],
$F_{c}$ was experimentally demonstrated to estimate the fidelity in a Rydberg
atom array quantum simulator, at infinite effective temperature and short
quench evolution times. $F_{d}$ can be taken as a further improvement of
$F_{c}$, that works for arbitrary quench initial states and Hamiltonians.

##  Appendix B Emergence of universal statistical properties from Hamiltonian
dynamics

In this section, we present precise statements of our Theorem in the main text
and prove them. The significance of our Theorems 1 and 2 is that a Porter-
Thomas distribution emerges from measurement outcomes
$p(z,t)=\absolutevalue{\innerproduct{z}{\Psi(t)}}^{2}$ after generic
Hamiltonian dynamics, upon appropriate data processing. To prove our theorems,
we first establish a new identity, the $k$-th Hamiltonian twirling identity,
which only requires one assumption: the eigenvalues of the Hamiltonian $H$
satisfy the $k$-th no-resonance condition. We utilize this identity to
establish theorems on the emergence of a Porter-Thomas distribution in the
rescaled measurement outcomes $\tilde{p}(z,t)$.

###  B.1 $k$-th no-resonance condition

Definition ($k$-th no-resonance): Let $E_{j}$ be the energy eigenvalues of a
Hamiltonian $H$. $H$ satisfies the $k$-th no-resonance condition [46, 47, 48,
49] if

$E_{\alpha_{1}}+E_{\alpha_{2}}+\cdots+E_{\alpha_{k}}=E_{\beta_{1}}+E_{\beta_{2}}+\cdots+E_{\beta_{k}}\leavevmode\nobreak\
$ (S3)

implies that the indices $(\alpha_{1},\alpha_{2},\dots,\alpha_{k})$ and
$(\beta_{1},\beta_{2},\dots,\beta_{k})$ are equal up to reordering. In other
words, $\beta_{j}=\alpha_{\sigma(j)}$ for some permutation $\sigma\in S_{k}$,
where $S_{k}$ is the group of permutations of $k$ elements (the symmetric
group). It is easy to see that if the $k$-th no-resonance condition is
satisfied, then so are the $k^{\prime}$-th no-resonance conditions, for
$k^{\prime}<k$. This assumption is considered generically satisfied: it is
believed that if the condition does not hold in a given ergodic system, any
small perturbation will generically break any $k$-th resonances [48].

###  B.2 $k$-th no resonance with degeneracies

The $k$-th no-resonance conditions are, strictly speaking, violated if the
Hamiltonian contains degeneracies: $E_{i}=E_{j},i\neq j$. This violates the
first no-resonance condition and so all higher $k$-th no-resonance conditions.
This seems to preclude systems with symmetry-protected degeneracies, which can
still be ergodic. For our purposes, however, we will consider pure states
$\rho=\outerproduct{\Psi_{0}}{\Psi_{0}}$. In this case, one can replace any
set of degenerate states $\mathcal{E}=\\{\ket{\Psi_{E}}\\}$ with a _single_
representative state, defined as the projection of $\ket{\Psi_{0}}$ onto the
degenerate subspace: $\mathcal{P}_{\mathcal{E}}\ket{\Psi_{0}}$. In this way,
all degeneracies are effectively eliminated. Therefore, while Lemma 1 (below)
requires the full $k$-th no-resonance condition, one can modify this condition
by demanding the $k$-th no-resonance condition only after eliminating all
degenerates as prescribed above. Our main Theorems only require this weakened
version of the $k$-th no-resonance condition. Hence, our Theorems holds
generically even in the presence of degeneracies protected by symmetries.

The $k$-th no-resonance condition is also notably violated for $k\geq 2$ if
there is a particle-hole symmetry, in which the energy spectrum is symmetric
about $E=0$. This case is discussed in Section F.1.

We finally remark that there are other, indirect ways in which degeneracies
can affect the performance of our benchmark. For example, if there is a large
degeneracy of eigenstates, the quench dynamics will be less ergodic and our
protocol will have larger uncertainties. However, this is entirely captured by
our effective dimension parameter $D_{\beta}$, and no further modification is
required.

###  B.3 $k$-th Hamiltonian twirling identity

We first present the $k$-th Hamiltonian twirling identity. This identity is a
rewriting of a $k$-fold quantum channel: we take a state $\rho^{(k)}$ in
$k$-copies of the Hilbert space $\mathcal{H}^{\otimes k}$, and average it over
all possible times of time-evolution $U_{t}^{\otimes k}$ over $k$ copies, with
$U_{t}\equiv\exp(-iHt)$. That is, the channel $C^{(k)}$ is
$C^{(k)}[\rho^{(k)}]=\mathbb{E}_{t}\left[U_{t}^{\otimes
k}\rho^{(k)}U_{-t}^{\otimes k}\right]$ where
$\mathbb{E}_{t}[f(t)]\equiv\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(t)dt$
here and below. We dub this the $k$-th Hamiltonian twirling channel, in
analogy to twirling channels, in which a state is averaged over all possible
unitaries in a group, most commonly the ensemble of Haar unitaries or the
Clifford or Pauli groups [59].

We find that the $k$-th Hamiltonian twirling channel can be written as sum of
(i) a dephasing operator in the energy eigenbasis, followed by a sum of
permutations of $k$ elements and (ii) additional terms associated with the
double counting of certain permutations. We associate these terms with
additional correlations associated with energy conservation, hence we dub
these “energy-conservation terms”.

Such a rewriting is motivated by Haar-unitary twirling channels, which can be
written as a sum of permutations [68]. In our Theorem below, we consider a
specific quantity: $\mathbb{E}_{t}[p^{k}(z,t)]=\bra{z}^{\otimes
k}C^{(k)}[\outerproduct{\Psi_{0}}{\Psi_{0}}^{\otimes k}]\ket{z}^{\otimes k}$.
In this case, the terms in (i) are the dominant contributions and we are able
to bound the magnitude of the energy-conservation terms in (ii).

We explicitly state our identity for $k=2$.

######  Lemma (Second Hamiltonian twirling identity).

If a Hamiltonian $H$ satisfies the second no-resonance condition, the
following identity holds:

$\displaystyle
C^{(2)}\big{[}\rho^{(2)}\big{]}\equiv\mathbb{E}_{t}\left[U_{t}^{\otimes
2}\rho^{(2)}U_{-t}^{\otimes 2}\right]$ (S4)
$\displaystyle=\sum_{E,E^{\prime}}\bigg{[}\outerproduct{E,E^{\prime}}{E,E^{\prime}}\rho^{(2)}\outerproduct{E,E^{\prime}}{E,E^{\prime}}+\outerproduct{E,E^{\prime}}{E,E^{\prime}}\rho^{(2)}\outerproduct{E^{\prime},E}{E^{\prime},E}\bigg{]}-\sum_{E}\outerproduct{E,E}{E,E}\rho^{(2)}\outerproduct{E,E}{E,E}$
$\displaystyle=\mathcal{D}^{\otimes
2}\big{[}\rho^{(2)}\big{]}\mathbb{I}+\mathcal{D}^{\otimes
2}\big{[}\rho^{(2)}\mathcal{S}\big{]}\mathcal{S}-\sum_{E}\outerproduct{E,E}{E,E}\rho^{(2)}\outerproduct{E,E}{E,E}\leavevmode\nobreak\
,$

where $U_{t}\equiv\exp(-iHt)$,
$\mathbb{E}_{t}[f(t)]=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(t)dt$
is the average of $f(t)$ over infinite time. $\rho^{(2)}$ is any state on two
copies of the Hilbert space, and $\ket{E}$ are energy eigenstates of the
Hamiltonian. $\mathbb{I}$ and $\mathcal{S}$ are the identity and swap
operators on $\mathcal{H}^{\otimes 2}$, defined as
$\mathbb{I}(\ket{\psi_{1}}\otimes\ket{\psi_{2}})=\ket{\psi_{1}}\otimes\ket{\psi_{2}}$
and
$\mathcal{S}(\ket{\psi_{1}}\otimes\ket{\psi_{2}})=\ket{\psi_{2}}\otimes\ket{\psi_{1}}$
for any $\ket{\psi_{1}},\ket{\psi_{2}}$. Lastly, $\mathcal{D}^{\otimes 2}$ is
the dephasing channel in the energy eigenbasis: $\mathcal{D}^{\otimes
2}\big{[}\rho^{(2)}\big{]}\equiv\sum_{E,E^{\prime}}\outerproduct{E,E^{\prime}}{E,E^{\prime}}\rho^{(2)}\outerproduct{E,E^{\prime}}{E,E^{\prime}}$.
Finally, we remark that the state $\rho^{(2)}$ on $\mathcal{H}^{\otimes 2}$
can be generalized to any operator $A$ acting on $\mathcal{H}^{\otimes 2}$.

###### Proof.

$\displaystyle\mathbb{E}_{t}\left[U_{t}^{\otimes 2}\rho^{(2)}U_{-t}^{\otimes
2}\right]=\lim_{T\rightarrow\infty}\frac{1}{T}\int^{T}_{0}dt\leavevmode\nobreak\
U_{t}^{\otimes 2}\rho^{(2)}U_{-t}^{\otimes 2}$ (S5)
$\displaystyle=\sum_{E_{1},E_{2},E_{3},E_{4}}\left[\lim_{T\rightarrow\infty}\frac{1}{T}\int^{T}_{0}dt\leavevmode\nobreak\
e^{-i(E_{1}-E_{2}+E_{3}-E_{4})t}\right]\outerproduct{E_{1},E_{3}}{E_{1},E_{3}}\rho^{(2)}\outerproduct{E_{2},E_{4}}{E_{2},E_{4}}$
(S6)
$\displaystyle=\sum_{E_{1},E_{2},E_{3},E_{4}}\delta(E_{1}-E_{2}+E_{3}-E_{4})\leavevmode\nobreak\
\outerproduct{E_{1},E_{3}}{E_{1},E_{3}}\rho^{(2)}\outerproduct{E_{2},E_{4}}{E_{2},E_{4}}$
(S7)
$\displaystyle\overset{\text{$2$-NR}}{=}\sum_{E,E^{\prime}}\bigg{[}\outerproduct{E,E^{\prime}}{E,E^{\prime}}\rho^{(2)}\outerproduct{E,E^{\prime}}{E,E^{\prime}}+\outerproduct{E,E^{\prime}}{E,E^{\prime}}\rho^{(2)}\outerproduct{E^{\prime},E}{E^{\prime},E}\bigg{]}-\sum_{E}\outerproduct{E,E}{E,E}\rho^{(2)}\outerproduct{E,E}{E,E}\leavevmode\nobreak\
,$ (S8)

where “(2-NR)” indicates the use of the second no-resonance condition. The
second no-resonance condition states that the delta function
$\delta(E_{1}-E_{2}+E_{3}-E_{4})$ is non-zero in only two cases:
$E_{1}=E_{2},E_{3}=E_{4}$ or $E_{1}=E_{4},E_{3}=E_{2}$, which give the two
terms above. The final, energy-conservation term in Eq. (S8) corresponds to
$E_{1}=E_{2}=E_{3}=E_{4}$. This special case belongs to the above two terms
simultaneously and is double counted. Therefore, the energy-conservation term
is subtracted to correct for this double counting.∎

A diagrammatic representation of this Lemma is shown in Fig. S1(a,b). The
energy-conservation terms are illustrated in Fig. S1(c). As with the second
Hamiltonian twirling identity, one can define the $k$-th twirling identity.

######  Lemma 1 ($k$-th Hamiltonian twirling identity).

If the Hamiltonian satisfies the $k$-th no-resonance condition, its $k$-th
twirling channel is:

$\displaystyle C^{(k)}[\rho^{(k)}]$
$\displaystyle\equiv\mathbb{E}_{t}\left[U_{t}^{\otimes
k}\rho^{(k)}U_{-t}^{\otimes k}\right]$
$\displaystyle=\sum_{M=(E_{1},\cdots,E_{k})}\bigg{[}\sum_{\sigma\in\pi(M)}\outerproduct{E_{1},\cdots,E_{k}}{E_{1},\cdots,E_{k}}\rho^{(k)}\outerproduct{E_{\sigma(1)},\cdots,E_{\sigma(k)}}{E_{\sigma(1)},\cdots,E_{\sigma(k)}}\bigg{]},$
(S9) $\displaystyle=\mathcal{D}^{\otimes
k}\big{[}\rho^{(k)}\hat{\sigma}^{-1}\big{]}\sum_{\sigma\in
S_{k}}\hat{\sigma}+\mathrm{E.C.}$ (S10)

where
$\mathbb{E}_{t}[f(t)]=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(t)dt$
is averaging over infinite time, $\rho^{(k)}$ is any state on $k$ copies of
the Hilbert space $\mathcal{H}^{\otimes k}$, and $\hat{\sigma}$—with
$\sigma\in S_{k}$—is a permutation operator acting on the $k$-copy Hilbert
space:

$\displaystyle\hat{\sigma}=\sum_{\\{E_{i}\\}}\outerproduct{E_{1},E_{2},\dots,E_{k}}{E_{\sigma(1)},E_{\sigma(2)},\dots,E_{\sigma(k)}}.$
(S11)

In (S9), the sum is taken over all possible lists of $k$ energies
$M=(E_{1},\cdots,E_{k})$. If $M$ has $r$ unique energies with multiplicities
$n_{1},\cdots,n_{r}$ (such that $\sum_{i=1}^{n}n_{r}=k$), its group of
permutations $\pi(M)$ has $k!/\prod_{i=1}^{n}n_{i}!$ elements and is a
subgroup of the symmetric group $S_{k}$. In (S10), this identity is re-written
in terms of permutations, the dephasing channel $\mathcal{D}^{\otimes
k}\big{[}\rho^{(k)}\big{]}\equiv\sum_{E_{1},\cdots,E_{k}}\outerproduct{E_{1},\cdots,E_{k}}{E_{1},\cdots,E_{k}}\rho^{(k)}\outerproduct{E_{1},\cdots,E_{k}}{E_{1},\cdots,E_{k}}$,
and energy-conservation terms E.C. associated with double, or multiple,
counting (similar to the $k=2$ case). All terms can be explicitly derived from
(S9).

The proof of the $k$-th twirling identity is analogous to that of the second
twirling identity: in this case, there are $k!$ permutations relevant for the
$k$-th no resonance theorem, along with more energy-conservation terms than
the simplest one in the $k=2$ case. As an example, the third twirling identity
is illustrated diagrammatically in Fig. S1(d). In addition to the energy-
conservation terms associated with cases where two energies are the same
$E_{i}=E_{j}$, there is a higher order energy-conservation term when all three
energies are the same: $E_{1}=E_{2}=E_{3}$.

Lastly, we note that this identity also holds in Floquet systems, as long the
$k$-th no resonance condition is tightened to
$\sum_{i}E_{\alpha_{i}}\equiv\sum_{i}E_{\beta_{i}}(\text{mod
}2\pi)\iff\\{\alpha_{i}\\}=\\{\beta_{i}\\}$, and the integral over all time is
replaced by a discrete sum over integer Floquet periods $nT$.

Figure S1: Second and third Hamiltonian twirling identities. (a) The second
Hamiltonian twirling identity comprises identity and swap permutations terms,
along with an energy-conservation correction term that in our context will
have an exponentially small contribution. (b) The same as (a), but open legs
representing the channel input/outputs are rearranged in order to simplify the
diagram. We keep this convention for the rest of diagrams in this document.
(c) Diagrammatic notation for the basic building block in our diagrams: the
dephasing operator in the energy eigenbasis:
$D[\rho]=\sum_{E}\outerproduct{E}{E}\rho\outerproduct{E}{E}$. (d) The third
Hamiltonian twirling identity comprises of six permutation and several energy-
conservation terms. For our purposes, these additional terms will have small
contributions. Figure S2: (a) Outcome probabilities $p(z,t)$ over time, for
three outcomes $z_{1},z_{2},z_{3}$ in a Bose-Hubbard model with 8 particles on
8 sites in red, green and blue respectively. For each $z$, $p(z,t)$ fluctuates
about its mean value $p_{\text{avg}}(z)$ (dashed lines). (b) Plotting the
histogram of $p(z,t)$ over time, each histogram follows an exponential
distribution with differing means $p_{\text{avg}}(z)$. (c) $p_{\text{avg}}(z)$
is approximately exponentially dependent on the outcome energy
$\bra{z}H\ket{z}$, analogous to a Boltzmann distribution. (d) The rescaled
outcome probabilities $\tilde{p}(z,t)$ follow (e) a universal exponential
(Porter Thomas) distribution, which follows from Theorem 1. (f) Taking a
histogram at fixed time [purple, (d)] over bitstrings $z$ gives the same
distribution (Theorem 2).

###  B.4 Porter-Thomas distribution over time

We will use the $k$-th Hamiltonian twirling identity, with $\rho^{\otimes
k}=\outerproduct{\Psi_{0}}{\Psi_{0}}^{\otimes k}$, to prove that the quantity
$\tilde{p}(z,t)\equiv p(z,t)/p_{\text{avg}}(z)$ asymptotically follows the
Porter-Thomas, or exponential distribution [39].

The Porter-Thomas (PT) distribution naturally arises from the distribution of
overlaps $p$ of the Haar ensemble — the unique distribution of states on the
Hilbert space that is invariant under any unitary — to a fixed vector
$\ket{\Phi}$: $p=\absolutevalue{\langle\Psi|\Phi\rangle}^{2}$, with
$\ket{\Psi}$ drawn from the Haar ensemble. A given $\ket{\Psi}$ and
$\ket{\Phi}$ have overlap $\absolutevalue{\langle\Psi|\Phi\rangle}^{2}$. The
probability density of obtaining an overlap with value $p$ is $P(p)=\int
d\Psi\delta(\absolutevalue{\langle\Psi|\Phi\rangle}^{2}-p)$. In deep random
circuits, the wavefunction $\ket{\Psi}$ is approximately Haar-random and we
can evaluate $P(p)$ for any Hilbert space dimension $D$, which converges to
the Porter-Thomas distribution in the limit $D\rightarrow\infty$.

$\displaystyle
P(p)dp=(D-1)(1-p)^{D-2}dp\xrightarrow{D\rightarrow\infty}\exp(-Dp)d(Dp)\leavevmode\nobreak\
,$ (S12)

It is natural to rescale the quantity $p$ by $D$. In this case, the
exponential distribution has a $k$-th moment
$\int_{0}^{\infty}(Dp)^{k}P(Dp)d(Dp)=k!$.

We formally state the Theorem from the main text in Theorems 1 and 2 below.
Informally, our theorems show that the rescaled distribution
$\tilde{p}(z,t)\equiv p(z,t)/p_{\text{avg}}(z)$ asymptotically follows the PT
distribution. As a reminder,
$p(z,t)=\absolutevalue{\innerproduct{z}{\Psi(t)}}^{2}$ is the probability of
measuring outcome $z$, and
$p_{\text{avg}}(z)=\mathbb{E}_{t}\left[p(z,t)\right]$ is its time-averaged
value. By showing that the $k$-th moments of $\tilde{p}(z,t)$ are close to
$k!$, we show that $\tilde{p}(z,t)$ asymptotically approaches the Porter-
Thomas distribution. We can choose to fix the bitstring $z$, and consider the
ensemble over times $t$ [Theorem 1, Fig. S2(d,e)], or fix a large but
arbitrary time $t$, and consider the ensemble of $\tilde{p}(z,t)$ over
outcomes $z$ [Theorem 2, Fig. S2(d,f)]. For brevity, we will occasionally
denote these ensembles as $\\{\tilde{p}(z,t)\\}_{t}$ and
$\\{\tilde{p}(z,t)\\}_{z}$ respectively, with the subscript denoting the
variable over which the ensemble is formed (the other variable is kept fixed).
Formally stated, our first Theorem is the following:

######  Theorem 1 (Temporal Porter-Thomas distribution).

For an initial state $\ket{\Psi_{0}}$ and a Hamiltonian $H$ satisfying the
$k$-th no-resonance condition, the temporal $k$-th moment of $\tilde{p}(z,t)$
(with fixed $z$), approaches a universal value
$\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]=k!+O(D_{\beta}^{-1}(z))$.
$D_{\beta}(z)$ is the effective Hilbert space dimension of the outcome $z$
defined as
$D_{\beta}^{-1}(z)\equiv\left(\sum_{E}\absolutevalue{\innerproduct{z}{E}}^{4}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{4}\right)/p_{\text{avg}}(z)^{2}$.

###### Proof.

Our proof is a direct application of our $k$-th Hamiltonian twirling identity
on the temporal $k$-th moment
$\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]=k!+O(D_{\beta}^{-1})$: To
compute $\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]$, we apply Lemma 1 on
$\rho^{(k)}=\outerproduct{\Psi_{0}}{\Psi_{0}}^{\otimes k}$:

$\displaystyle\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]=\frac{\bra{z}^{\otimes
k}\mathbb{E}_{t}\left[U_{t}^{\otimes
k}\outerproduct{\Psi_{0}}{\Psi_{0}}^{\otimes k}U_{-t}^{\otimes
k}\right]\ket{z}^{\otimes k}}{p_{\text{avg}}(z)^{k}}$ (S13)
$\displaystyle=\frac{\sum_{M=(E_{1},\cdots,E_{k})}\sum_{\sigma\in\pi(M)}\leavevmode\nobreak\
\left[\prod_{i=1}^{k}\absolutevalue{\innerproduct{z}{E_{i}}\innerproduct{E_{\sigma_{i}}}{\Psi_{0}}}^{2}\right]}{p_{\text{avg}}(z)^{k}}=\frac{\sum_{M=(E_{1},\cdots,E_{k})}\absolutevalue{\pi(M)}\leavevmode\nobreak\
\left[\prod_{i=1}^{k}\absolutevalue{\innerproduct{z}{E_{i}}\innerproduct{E_{i}}{\Psi_{0}}}^{2}\right]}{p_{\text{avg}}(z)^{k}}$
(S14)

where we have used Eq. (S9). Now using Eq. (S10), where the contribution from
the first term is equivalent to replacing $\absolutevalue{\pi(M)}$ with $k!$.
This gives:

$\displaystyle\frac{\bra{z}^{\otimes k}\mathcal{D}^{\otimes
k}\left[\ket{\Psi_{0}}\bra{\Psi_{0}}^{\otimes k}\right]\left[\sum_{\sigma\in
S_{k}}\hat{\sigma}\right]\ket{z}^{\otimes
k}}{p_{\text{avg}}(z)^{k}}=k!\frac{\sum_{\\{E_{j}\\}}\prod_{j=1}^{k}\absolutevalue{\innerproduct{z}{E_{j}}\\!\innerproduct{E_{j}}{\Psi_{0}}}^{2}}{p_{\text{avg}}(z)^{k}}=k!\leavevmode\nobreak\
,$ (S15)

where in the first equality we used the fact that $\hat{\sigma}$ operator acts
trivially on the product state $\ket{z}^{\otimes k}$ and in the second
equality, we used the relation

$p_{\text{avg}}(z)\equiv\mathbb{E}_{t}\left[p(z,t)\right]=\sum_{E}\absolutevalue{\innerproduct{z}{E}\\!\innerproduct{E}{\Psi_{0}}}^{2}\leavevmode\nobreak\
.$ (S16)

This relation is used several times, motivating our definition of the quantity

$q_{z}(E)\equiv\frac{\absolutevalue{\innerproduct{z}{E}\\!\innerproduct{E}{\Psi_{0}}}^{2}}{p_{\text{avg}}(z)}$
(S17)

In fact, $q_{z}(E)$ has the properties of a probability distribution over the
energies $E$: $0\leq q_{z}(E)\leq 1$ and $\sum_{E}q_{z}(E)=1$.

Now we wish to bound the deviation
$\absolutevalue{\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]-k!}$, that is
the contribution from all energy-conservation (E.C.) terms. To do so, we note
the E.C. terms must come from lists of energies $M=(E_{1},\cdots,E_{k})$ with
at least one repeated energy. That is:

$\displaystyle\absolutevalue{\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]-k!}\equiv\absolutevalue{\frac{\bra{z}^{\otimes
k}(\text{E.C.})\ket{z}^{\otimes k}}{p_{\text{avg}}(z)^{k}}}$
$\displaystyle=\sum_{\begin{subarray}{c}M=(E_{1},\cdots,E_{k})\\\ \text{not
all $E_{j}$
distinct}\end{subarray}}\absolutevalue{\pi(M)-k!}\prod_{j=1}^{k}q_{z}(E_{j})$
(S18) $\displaystyle\leq
k!\left[\sum_{E_{1}=E_{2}}\sum_{E_{3},\cdots,E_{k}}\prod_{j=1}^{k}q_{z}(E_{j})+\sum_{E_{1}=E_{3}}\sum_{E_{2},E_{4},\cdots,E_{k}}\prod_{j=1}^{k}q_{z}(E_{j})+\cdots\right]$

In the last line, we have used two separate inequalities: we bound the
difference $\absolutevalue{\pi(M)-k!}\leq k!$. This enables us to contain all
possible lists $M$ with at least one repeated energy in the sums above, where
we have only explicitly demanded that one pair of energies coincide:
$E_{a}=E_{b}$. There are ${k\choose 2}$ choices of $a$ and $b$. Each choice
has the same value. Therefore without loss of generality, we set $a=1$ and
$b=2$. Then we obtain:

$\displaystyle\absolutevalue{\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]-k!}\leq
k!{k\choose
2}\sum_{E_{1}}q_{z}(E_{1})^{2}\prod_{j=3}^{k}\sum_{E_{j}}q_{z}(E_{j})=k!{k\choose
2}\sum_{E_{1}}q_{z}(E_{1})^{2}\equiv k!{k\choose
2}D_{\beta}^{-1}(z)\leavevmode\nobreak\ .$ (S19)

Therefore, we obtain the desired bound

$\absolutevalue{\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]-k!}\leq
k!{k\choose 2}D_{\beta}^{-1}(z)=O(D_{\beta}^{-1}(z))\leavevmode\nobreak\ .$
(S20)

∎

In fact, we expect this bound to be close to tight: the contribution from the
simplest energy-conservation terms — those with exactly one repeated energy —
is close to our bound. Other energy-conservation terms are expected to be
exponentially smaller than these simplest terms. Here, the factor
$\pi(M)=k!/2$, and we obtain:

$\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]\approx
k!\left(1-\frac{k(k-1)}{2}\frac{D_{\beta}^{-1}(z)}{2}\right)\leavevmode\nobreak\
.$ (S21)

Eq. (S21) agrees with the leading order (in $D^{-1}$) correction to the $k$-th
moment of the Porter-Thomas distribution:

$\displaystyle\int_{0}^{1}dp\leavevmode\nobreak\
p^{k}(D-1)(1-p)^{D-2}=k!\frac{D^{k}}{D\cdots(D+k-1)}=k!\left[1-\frac{k(k-1)}{2}D^{-1}+O(D^{-2})\right]\leavevmode\nobreak\
.$ (S22)

Therefore, it is apt to call $D_{\beta}(z)$ as the effective dimension: the
deviation of the $k$-th moment from $k!$ for the temporal distribution is the
same as that of the Porter-Thomas distribution for a system with finite
dimension $2D_{\beta}$. $\\{\tilde{p}(z,t)\\}_{t}$ approximately follows a
Porter-Thomas distribution over a vector space of dimension $2D_{\beta}(z)$.
Note, however, that this agreement does not extend to corrections of higher
orders in $D^{-1}$.

For a fixed $z$, the inverse effective dimension $D_{\beta}^{-1}(z)$ can be
understood as the collision probability associated with $q_{z}(E)$. More
explicitly, we can define $D_{\beta}^{-1}(z)=\sum_{E}q_{z}(E)^{2}$ such that
$D_{\beta}(z)$ describes the participation ratio (or equivalently the size of
effective Hilbert space dimension) corresponding to a particular bitstring
$z$. Then the overall inverse participation ratio (inverse effective
dimension) $D_{\beta}^{-1}$ is the weighted average:
$D_{\beta}^{-1}=\sum_{z}p_{\text{avg}}(z)D_{\beta}^{-1}(z)$.

###  B.5 Porter-Thomas distribution over outcomes

Theorem 2 concerns the distribution of $\tilde{p}(z,t)$ for fixed $t$, treated
as samples labeled by $z$.

######  Theorem 2 (Porter-Thomas distribution over outcomes).

Consider an initial state $\ket{\Psi_{0}}$ and a Hamiltonian $H$ satisfying
the $2k$-th no-resonance condition. At a typical time $t\in[0,\infty)$, the
weighted $k$-th moment of $\tilde{p}(z,t)$ approaches a universal value:
$P^{(k)}(t)\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{k}=k!+O(D_{\beta}^{-1/2})$,
where
$D_{\beta}^{-1}\equiv\sum_{z,E}\absolutevalue{\innerproduct{z}{E}}^{4}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{4}/p_{\text{avg}}(z)$
is the inverse effective Hilbert space dimension associated with the quench
evolution, and $\\{\ket{E}\\}$ are the eigenstates of $H$.

###### Proof.

Our proof strategy is the following. We first compute the time-averaged value
$\mathbb{E}_{t}\left[P^{(k)}(t)\right]$. From Theorem B, this quantity is
$k!+O(D_{\beta}^{-1})$. We next compute the fluctuations of $P^{(k)}(t)$ over
time, and we will find that
$\text{std}_{t}\left[P^{(k)}(t)\right]=O(D_{\beta}^{-1/2})$. Combining these
two results leads to our desired conclusion:
$P^{(k)}(t)=k!+O(D_{\beta}^{-1/2})$ with high probability.

In Fig. S3(a), we illustrate this computation of
$\mathbb{E}_{t}\left[P^{(k)}(t)\right]$ for the case $k=2$. In general, we
have:

$\absolutevalue{\mathbb{E}_{t}\left[P^{(k)}(t)\right]-k!}\leq\sum_{z}p_{\text{avg}}(z)\absolutevalue{\mathbb{E}_{t}\left[\tilde{p}(z,t)^{k}\right]-k!}\leq{k\choose
2}k!\sum_{z}p_{\text{avg}}(z)D_{\beta}^{-1}(z)=O(D_{\beta}^{-1})$ (S23)

Up to this point, we have only used the $k$-th Hamiltonian twirling identity.
To estimate the temporal fluctuations of $P^{(k)}(t)$, we need the $2k$-th
twirling identity.

$\displaystyle\text{var}_{t}[P^{(k)}(t)]$
$\displaystyle=\mathbb{E}_{t}[P^{(k)}(t)^{2}]-\mathbb{E}_{t}[P^{(k)}(t)]^{2}=\sum_{z,z^{\prime}}p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})\frac{\mathbb{E}_{t}[p(z,t)^{k}p(z^{\prime},t)^{k}]-\mathbb{E}_{t}[p(z,t)^{k}]\mathbb{E}_{t}[p(z^{\prime},t)^{k}]}{p_{\text{avg}}(z)^{k}p_{\text{avg}}({z^{\prime}})^{k}}$
(S24)

By using the Hamiltonian $k$-th and $2k$-th twirling identities, any term in
$\mathbb{E}_{t}[P^{(k)}(t)^{2}]$ that does not cancel with a term in
$\mathbb{E}_{t}[P^{(k)}(t)]^{2}$ must involve a permutation between $z$ and
$z^{\prime}$ [illustrated for $k=2$ in Fig. S3(b)]. As in our proof of Theorem
1, without loss of generality, we let the energies corresponding to this
permutation be $E_{1}$ and $E_{2}$, and multiply our overall expression by a
factor of $2{2k\choose 2}$ (here, $E_{1}$ and $E_{2}$ are inequivalent). Then
their sum is of the form

$\displaystyle\text{var}_{t}[P^{(k)}(t)]=\sum_{z,z^{\prime}}p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})2{2k\choose
2}\frac{\sum_{E_{1},E_{2}}\absolutevalue{\innerproduct{E_{1}}{\Psi_{0}}}^{2}\innerproduct{z}{E_{1}}\\!\innerproduct{E_{1}}{z^{\prime}}\absolutevalue{\innerproduct{E_{2}}{\Psi_{0}}}^{2}\innerproduct{z^{\prime}}{E_{2}}\\!\innerproduct{E_{2}}{z}}{p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})}(\cdots)$
(S25)
$\displaystyle\equiv\sum_{z,z^{\prime}}p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})2{2k\choose
2}\frac{\bra{z}\rho_{d}\outerproduct{z^{\prime}}{z^{\prime}}\rho_{d}\ket{z}}{p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})}(\cdots)\leavevmode\nobreak\
,$ (S26)
$\displaystyle(\cdots)=\sum_{\text{perms.}}\prod_{j=3}^{2k}\frac{\sum_{E_{j}}\absolutevalue{\innerproduct{E_{j}}{\Psi_{0}}}^{2}\innerproduct*{z_{j}^{L}}{E_{j}}\\!\innerproduct*{E_{j}}{z_{j}^{R}}}{\sqrt{p_{\text{avg}}({z_{j}^{L}})p_{\text{avg}}({z_{j}^{R}})}}\equiv\sum_{\text{perms.}}\prod_{j=3}^{2k}\frac{\bra*{z_{j}^{L}}\rho_{d}\ket*{z_{j}^{R}}}{\sqrt{p_{\text{avg}}({z_{j}^{L}})p_{\text{avg}}({z_{j}^{R}})}}\leq\sum_{\text{perms.}}1\leq(2k-2)!$
(S27)

where $z_{j}^{L}$ and $z_{j}^{R}$ can either be $z$ or $z^{\prime}$, such that
there exactly $k-1$ of the $\\{z_{j}^{L}\\}$ are $z$ and the other $k-1$ are
$z^{\prime}$, and likewise for the $\\{z_{j}^{R}\\}$. For brevity, we have
also rewritten our expressions with
$\rho_{d}\equiv\sum_{E}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{2}\outerproduct{E}{E}$,
the commonly studied diagonal ensemble [47, 53]. Note that
$p_{\text{avg}}(z)=\bra{z}\rho_{d}\ket{z}$.

The inequalities in (S27) are obtained as follows. The term involving
$\rho_{d}$ is bounded by a simple application of the Cauchy-Schwarz identity,
gives $\absolutevalue{\bra{z}\rho_{d}\ket{z^{\prime}}}^{2}\leq
p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})$. Since
$\absolutevalue{\bra{z}\rho_{d}\ket{z}}^{2}\equiv
p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})$, it follows that
$\absolutevalue{\bra*{z_{j}^{L}}\rho_{d}\ket*{z_{j}^{R}}}/\sqrt{p_{\text{avg}}({z_{j}^{L}})p_{\text{avg}}({z_{j}^{R}})}\leq
1$ regardless of the values of $z_{j}^{L}$ and $z_{j}^{R}$. Finally, the sum
$\sum_{\text{perms.}}$ is a sum over allowed permutations of a given list of
energies $(E_{1},\cdots,E_{2k})$; we bound this by $(2k-2)!$. (Formally,
$\sum_{\text{perms.}}$ should be inside the summation over energies, with
different permutations giving rise to different realizations of
$\\{z_{j}^{L},z_{j}^{R}\\}$. However, we use the notation above because this
sum only gives a constant factor.)

Combining (S26) and (S27), we obtain

$\displaystyle\text{var}_{t}[P^{(k)}(t)]=(2k)!\sum_{z,z^{\prime}}p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})\frac{\bra{z}\rho_{d}\outerproduct{z^{\prime}}{z^{\prime}}\rho_{d}\ket{z}}{p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})}=(2k)!\tr(\rho_{d}^{2})\leq(2k)!D_{\beta}^{-1}$
(S28)

In our last inequality, we used $\tr(\rho_{d}^{2})\leq D_{\beta}^{-1}$. This
is derived from Sedrakyan’s inequality, which states that for positive $u_{i}$
and $v_{i}$,
$\left(\sum_{i}u_{i}\right)^{2}/\sum_{i}v_{i}\leq\sum_{i}\left(u_{i}^{2}/v_{i}\right)$.
In our context,

$\tr(\rho_{d}^{2})\equiv\sum_{E}\absolutevalue{\innerproduct{\Psi_{0}}{E}}^{4}=\sum_{E}\absolutevalue{\innerproduct{\Psi_{0}}{E}}^{4}\frac{\left(\sum_{z}\absolutevalue{\innerproduct{z}{E}}^{2}\right)^{2}}{\sum_{z}p_{\text{avg}}(z)}\leq\sum_{z,E}\frac{\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{4}}{p_{\text{avg}}(z)}\equiv
D_{\beta}^{-1}\leavevmode\nobreak\ .$ (S29)

This completes our proof that for a typical $t$,
$P^{(k)}(t)=k!+O(D_{\beta}^{-1/2})$. ∎

Our theorem states that when $D_{\beta}$ is sufficiently large, the $k$-th
moment of $\\{p(z,t)\\}$ is close to $k!$. We can then conclude that
$\\{p(z,t)\\}$ follows the exponential distribution when $D_{\beta}$ is
sufficiently large. Our numerical results support this: fixing a late enough
time $t$ (we numerically find that a time scaling linearly in system size is
sufficient for our purposes) and plotting the histogram of the rescaled
probabilities $\\{\tilde{p}(z,t)\\}_{z}$ (weighted by $p_{\text{avg}}(z)$)
shows a good agreement with the exponential distribution, as seen in Fig.
S2(f).

Theorems 1 and 2 together suggest an analogue to the ergodic theorem in
classical dynamical systems. Consider an ensemble of particles in a classical
dynamical system, such as a billiard stadium. The ergodic theorem states that
the spatial distribution of this ensemble at a fixed time is the same as the
distribution of a fixed particle over time. In our context, the distributions
of $\tilde{p}(z,t)$ are asymptotically the same in two settings: with fixed
$t$ over configurations $z$ and with fixed $z$ over $t$. Furthermore, both
distributions limit to the universal Porter-Thomas distribution.

##  Appendix C Detailed performance analysis for long quench dynamics

In this section we apply the above results to characterize the accuracy of our
fidelity estimator $F_{d}$. Our analysis is valid when the state has reached
global equilibrium, quantified by the quantity
$P^{(2)}(t)\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{2}$ reaching its
equilibrium value of 2. In Theorem 2, we have bounded the deviation of the
$k$-th moments $P^{(k)}(t)$ from its ideal value of $k!$. Under an
approximation [Eq. (S35)], our bound for $P^{(2)}(t)$ is sufficient to claim
that the $F_{d}$ formula approximates $F$ with accuracy that is exponentially
good in system size. By keeping the simplest energy-conservation terms that
are the leading-order corrections in Theorem 2, we make quantitative
predictions for the systematic error $\delta_{\textrm{sys}}$ and temporal
fluctuations $\delta_{\textrm{temp}}$. We explicitly confirm our theoretical
analysis via numerical simulations presented in Fig. 3 of our main text and
Fig. S4, which shows excellent agreement between theory prediction and
numerical results. This suggests that our formalism presented here is able to
capture the performance of the $F_{d}$ formula to first order in
$D_{\beta}^{-1}$, at least for the system sizes investigated.

We define the following relevant quantities:

$\displaystyle F_{\text{id,d}}(t)$
$\displaystyle\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{2}-1=P^{(2)}(t)-1\leavevmode\nobreak\
,$ (S30) $\displaystyle F_{\text{XEB,d}}(t)$
$\displaystyle\equiv\sum_{z}q(z,t)\tilde{p}(z,t)-1\leavevmode\nobreak\ .$
(S31)

$F_{d}$ is given by

$F_{d}=2\frac{\sum_{z}q(z)\tilde{p}(z)}{\sum_{z}p(z)\tilde{p}(z)}-1=2\frac{F_{\text{XEB,d}}+1}{F_{\text{id,d}}+1}-1$
(S32)

$F_{\text{id,d}}$ is a measure of whether the state has reached global
equilibrium: after equilibration, $F_{\text{id,d}}\approx 1$, while at very
short times, $F_{\text{id,d}}\sim O(D)$ is exponentially large when the
initial state is a product state or a lowly-entangled state.
$F_{\text{XEB,d}}$ is analogous to the linear cross-entropy benchmark
$F_{\text{XEB}}$.

We can bound $F_{\text{id,d}}\geq 0$ and $F_{\text{XEB,d}}\geq-1$, and we
obtain (weak) bounds on $F_{d}$:

$\displaystyle-1\leq$ $\displaystyle F_{d}\leq
2\sqrt{\frac{\sum_{z}p_{\text{avg}}(z)\tilde{q}(z)^{2}}{\sum_{z}p_{\text{avg}}(z)\tilde{p}(z)^{2}}}-1\leavevmode\nobreak\
,$ (S33)

We compute the average values over time of $F_{d}$ as well as its variance
over time. We find that the typical variance is exponentially small in system
size, therefore sampling at a single late time will give values close to its
mean.

We use “exponentially small in system size” as a shorthand to mean that a
quantity is proportional to a measure of the effective dimension:
$D_{\beta}^{-1}$ or $\tr(\rho_{d}^{2})$. These are related by the following
inequalities:

$\displaystyle
D^{-1}\leq\sum_{z}p_{\text{avg}}(z)^{2}\leq\tr(\rho_{d}^{2})\leq
D_{\beta}^{-1}\leq 1\leavevmode\nobreak\ .$ (S34)

The first inequality is a straightfoward inequality between the arithmetic and
quadratic means. The second inequality is obtained immediately by noting that
$\tr(\rho_{d})^{2}=\sum_{z}p_{\text{avg}}(z)^{2}+\sum_{z\neq
z^{\prime}}\absolutevalue{\bra{z}\rho_{d}\ket{z^{\prime}}}^{2}\geq\sum_{z}p_{\text{avg}}(z)^{2}$
, and the final inequality was established in Eq. (S29). $D_{\beta}^{-1}$ is
in general only upper bounded by $1$, since choosing
$\ket{z}=\ket{E}=\ket{\Psi_{0}}$, gives $D_{\beta}^{-1}=1$. However, for
generic initial states, we expect $D_{\beta}^{-1}$ to be exponentially small
in system size.

###  C.1 Ansatz for how speckle patterns change under errors

In the long-time limit, it is fruitful to approximate $q(z)$ as:

$q(z)\approx Fp(z)+(1-F)p_{\perp}(z)\leavevmode\nobreak\ ,$ (S35)

where $p_{\perp}(z)$ is a distribution uncorrelated with $p(z)$. In general,
in the case of incoherent errors, $\tilde{p}_{\perp}(z)\equiv
p_{\perp}(z)/p_{\text{avg}}(z)$ does not follow a PT distribution.

We assume that $p(z)$ and $p_{\perp}(z)$ are uncorrelated. Specifically, we
assume that:

$\displaystyle\sum_{z}p_{\perp}(z)\tilde{p}(z)^{k}$
$\displaystyle\approx\mathbb{E}_{z}[p_{\perp}(z)]\sum_{z^{\prime}}\tilde{p}(z^{\prime})^{k}\approx
k!\leavevmode\nobreak\ ,$ (S36)

This illustrates that we do not require $\tilde{p}_{\perp}(z)$ to follow the
Porter-Thomas distribution — the energy density of the perturbed state can be
quite different from the target state, as is the case when many errors occur.
As an example, applying this on Eq. (S31) gives the desired
$F_{\text{XEB,d}}\approx F$.

A related approximation applicable to higher order correction terms, is to
assume a coherent error $V$, and write the experimental state as
$\ket{\psi^{\prime}(t)}=V(t)\ket{\psi_{0}(t)}$. $V(t)$ is unitary, and
therefore $V^{\dagger}(t)V(t)=\mathbb{I}$. However, we replace instances of
$V(t)$ or $V^{\dagger}(t)$ that do not cancel with their Hermitian conjugates
with

$V(t),V^{\dagger}(t)\rightarrow\sqrt{F}\leavevmode\nobreak\ .$ (S37)

While crude, this approximation is useful to make quick estimates.
Furthermore, its predictions quantitatively agree with our numerical results,
indicating that the dominant sources of error are due to the structure of the
$F_{d}$ formula, which we discuss below.

###  C.2 Systematic error

Here, we characterize the systematic error
$\delta_{\text{sys.}}\equiv\mathbb{E}_{t}\left[F_{d}(t)\right]-F$: the
difference between the time-averaged $F_{d}$ and the true fidelity $F$ after
an isolated error. To do so, we first study the time average of the simpler
quantities
$F_{\text{id,d}}(t)\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{2}-1$ and
$F_{\text{XEB,d}}(t)\equiv\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)\tilde{q}(z,t)-1$.

The $k$-th moment of the PT distribution (Theorems 1,2) immediately allows us
to estimate $F_{\text{id,d}}$:

$\displaystyle\mathbb{E}_{t}\left[F_{\text{id,d}}(t)\right]$
$\displaystyle=\mathbb{E}_{t}\left[\sum_{z}p_{\text{avg}}(z)\tilde{p}^{2}(z,t)-1\right]=1-D_{\beta}^{-1}\leavevmode\nobreak\
.$ (S38)

Using (S37) gives

$\displaystyle F_{\text{XEB,d}}\approx F(1-D_{\beta}^{-1})\leavevmode\nobreak\
.$ (S39)

We would like to use these results to obtain the time-averaged value of
$F_{d}$. The time average of the numerators and denominators
$F_{\text{XEB,d}}$ and $F_{\text{id,d}}$ should be done together. We can
estimate its effects with the following Taylor expansion:

$\displaystyle\mathbb{E}\left[\frac{X}{Y}\right]\approx\frac{\mathbb{E}[X]}{\mathbb{E}[Y]}-\frac{\text{cov}[X,Y]}{\mathbb{E}[Y]^{2}}+\frac{\text{var}[Y]\mathbb{E}[X]}{\mathbb{E}[Y]^{3}}\leavevmode\nobreak\
,$ (S40)

which gives the estimate for the systematic error in the presence of a single
error:

$\displaystyle\delta_{\textrm{sys}}=\mathbb{E}_{t}\left[F_{d}\right]-F$
$\displaystyle\approx(1-F)\left[\frac{D_{\beta}^{-1}}{2}+\frac{\text{var}_{t}\left[F_{\text{id,d}}\right]}{8}\right]\leavevmode\nobreak\
.\leavevmode\nobreak\ $ (S41)

This prediction agrees well with numerical simulations in Fig. S4(a) and Fig.
3(b) in the main text.

####  C.2.1 Argument based on Eigenstate Thermalization Hypothesis

In the limit of low error rates and long-time evolution, we can show that
$F_{d}\approx F$, independent of our approximation for $\tilde{q}$ (S35). To
do so, we consider the effect of a single error $V$ that occurs at a variable
time $t$. In the limit of long evolution, $F_{d}\approx
F_{\text{XEB,d}}\equiv\sum_{z}p_{\text{avg}}(z)\tilde{q}(z)\tilde{p}(z)-1$
since the denominator in $F_{d}$ approaches 2. Therefore, we will show that
$F_{\text{XEB,d}}\approx F$ for a single error averaged over the error
occurrence time. In the presence of multiple errors, the above relation will
hold, as long as the errors are sufficiently far apart in time. We average
over two quantities: (i) the time $t$ at which the error $V$ occurs, and (ii)
the time $\tau$ after the error happens when the measurement is made. At this
time, the ideal state is $\ket{\Psi(t+\tau)}$, while the perturbed state is
$V(\tau)\ket{\Psi(t+\tau)}$, with $V(\tau)\equiv\exp(-iH\tau)V\exp(iH\tau)$
the Heisenberg evolution of $V$. Here, the fidelity
$F(t)=\absolutevalue{\bra{\Psi(t+\tau)}V(\tau)\ket{\Psi(t+\tau)}}^{2}=\absolutevalue{\bra{\Psi(t)}V\ket{\Psi(t)}}^{2}$
is independent of $\tau$ and only depends on $t$ (the time when the error is
applied), while the XEB depends on both $t$ and $\tau$. We will average both
$t$ and $\tau$. We first average over $t$ to compute the average value of
$F_{\text{XEB,d}}$ at a time $\tau$ after an error is applied:

$\displaystyle\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]=\sum_{z}\frac{\absolutevalue{\bra{z}V(\tau)\rho_{d}\ket{z}}^{2}}{p_{\text{avg}}(z)}-\sum_{z}\frac{\bra{z,z}\left(V^{\dagger}(\tau)\otimes\mathbb{I}\right)\rho_{d}^{(2)}\left(V(\tau)\otimes\mathbb{I}\right)\ket{z,z}}{p_{\text{avg}}(z)}\leavevmode\nobreak\
,$ (S42)

where we have defined the diagonal ensemble $\rho_{d}$ and its analog
$\rho_{d}^{(2)}$ over $\mathcal{H}^{\otimes 2}$ as follows:

$\displaystyle\rho_{d}$
$\displaystyle\equiv\sum_{E}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{2}\outerproduct{E}{E}=\mathbb{E}_{t}\left[\outerproduct{\Psi(t)}{\Psi(t)}\right]$
(S43) $\displaystyle\rho_{d}^{(2)}$
$\displaystyle\equiv\sum_{E}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{4}\outerproduct{E,E}{E,E}$
(S44)

$\rho_{d}^{(2)}$ is an analog of $\rho_{d}$ defined over two-copies of the
Hilbert space. The first term in Eq. (S42) is the dominant contribution, while
the second term comes from the energy-conservation term in Eq. (S4).

We next average over the time $\tau$ after the error, in order to claim that
$\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]\approx\mathbb{E}_{t}\left[F(t)\right]$.
The averaging over $\tau$ gives:

$\displaystyle\mathbb{E}_{\tau}\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]$
$\displaystyle=\sum_{z}\frac{\bra{z}\rho_{d}V_{d}\outerproduct{z}{z}\rho_{d}V^{\dagger}_{d}\ket{z}}{p_{\text{avg}}(z)}-\sum_{z}\frac{\bra{z,z}\rho^{(2)}_{d}\left(V_{d}\otimes
V_{d}^{\dagger}\right)\ket{z,z}}{p_{\text{avg}}(z)}\leavevmode\nobreak\ ,$
(S45)

where $V_{d}\equiv\sum_{E}\outerproduct{E}{E}V\outerproduct{E}{E}$ is the
error operator dephased in the energy basis, or the time-averaged operator
$\mathbb{E}_{\tau}V(\tau)$. We interpret the first term — the dominant
contribution to (S51) — as the average fidelity of states
$\ket{\tilde{z}}\propto\sqrt{\rho_{d}}\ket{z}$ after the error operator
$V_{d}$:

$\sum_{z}\frac{\bra{z}\rho_{d}V_{d}\outerproduct{z}{z}\rho_{d}V^{\dagger}_{d}\ket{z}}{p_{\text{avg}}(z)}=\sum_{z}\frac{\bra{z}\sqrt{\rho_{d}}V_{d}\sqrt{\rho_{d}}\outerproduct{z}{z}\sqrt{\rho_{d}}V^{\dagger}_{d}\sqrt{\rho_{d}}\ket{z}}{p_{\text{avg}}(z)}=\sum_{z}p_{\text{avg}}(z)\absolutevalue{\bra{\tilde{z}}V_{d}\ket{\tilde{z}}}^{2}\leavevmode\nobreak\
,$ (S46)

$\ket{\tilde{z}}\equiv\sqrt{\rho_{d}}\ket{z}/\sqrt{p_{\text{avg}}(z)}$ is a
normalized wavefunction, after a transformation of a computational basis
state. We interpret this transformation as an imaginary time evolution
generated by the modular Hamiltonian of $\rho$. We now interpret $F_{d}$ as an
average fidelity $\absolutevalue{\bra{\tilde{z}}V_{d}\ket{\tilde{z}}}^{2}$ of
an error operator $V_{d}$ is applied to $\ket{\tilde{z}}$.

We want to relate
$\mathbb{E}_{\tau}\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]$ to the
fidelity $\mathbb{E}_{t}\left[F(t)\right]$ averaged over the time $t$ that the
error occurred. This is given by:

$\displaystyle\mathbb{E}_{t}\left[F(t)\right]$
$\displaystyle=\mathbb{E}_{t}\left[\absolutevalue{\bra{\Psi(t)}V\ket{\Psi(t)}}^{2}\right]=\absolutevalue{\tr(\rho_{d}V)}^{2}+\tr(\rho_{d}V\rho_{d}V^{\dagger})-\tr[\rho_{d}^{(2)}(V\otimes
V^{\dagger})]\leavevmode\nobreak\ ,$ (S47)

with the terms ordered by (expected) magnitude.

In order to relate Eq. (S42) and (S47), we make two observations: (i) for a
product initial state $\ket{\Psi_{0}}$ and local Hamiltonian $H$, $\rho_{d}$
approaches a microcanonical state since its energy variance vanishes relative
to the size of the spectrum. This is because $H$ is local, implying that
$\tr(H^{2})\propto N^{2}$. Therefore, the eigenvalues of $H$ are spread over
an energies scaling with total system size $N$. Meanwhile, (ii) by the
locality of $H$ the uncertainty in energy of a product state is of order
$\sqrt{N}$: $\Delta
E^{2}\equiv\bra{\Psi_{0}}H^{2}\ket{\Psi_{0}}-\bra{\Psi_{0}}H\ket{\Psi_{0}}^{2}\propto
N$. We model the distribution of energies as a Gaussian distribution with mean
$\bar{E}\equiv\bra{\Psi_{0}}H\ket{\Psi_{0}}$ and variance $\Delta E^{2}$:

$\sum_{E=\bar{E}-dE/2}^{\bar{E}+dE/2}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{2}\approx\rho\left(\bar{E}\right)\overline{\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{2}}dE\approx\frac{1}{\Delta
E\sqrt{2\pi}}\exp(-\frac{(E-\bar{E})^{2}}{2(\Delta
E)^{2}})dE\leavevmode\nobreak\ ,$ (S48)

where $\rho(\bar{E})$ is the density of states around $\bar{E}$. We define
this such that $\int\rho(E)dE=D$, such that $\rho(\bar{E})=\exp(cN)$ for some
$c>0$. Eq. (S48) also gives the estimate
$\sum_{E=\bar{E}-dE/2}^{\bar{E}+dE/2}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{4}\propto
1/[\rho(\bar{E})\Delta E]\sim\exp(-cN)$ (omitting polynomial factors of $N$).

(iii) Lastly, the Eigenstate Thermalization Hypothesis (ETH) states that the
diagonal elements $\bra{E}V\ket{E}$ of a local operator $V$ are well described
by a smooth function of $E$, $\bra{E}V\ket{E}\approx v(E)$. The off-diagonal
elements are exponentially smaller:
$\absolutevalue{\bra{E}V\ket{E^{\prime}}}^{2}\approx f(E-E^{\prime})\langle
V^{2}\rangle_{\text{micro}}\exp(-S(\bar{E})/2)$, for some function $f$ which
is $O(1)$ when the separation $E-E^{\prime}$ is small, and
$\bar{E}=(E+E^{\prime})/2$. $\langle V^{2}\rangle_{\text{micro}}$ is the
microcanonical value of $V^{2}$ near $\bar{E}$, and
$S(\bar{E})=-\ln(\rho(\bar{E}))$ is the thermodynamic entropy [60]. Then the
expression (S47) for the fidelity can be approximated as follows:

$\displaystyle\mathbb{E}_{t}\left[F\right]$ $\displaystyle\leavevmode\nobreak\
=\absolutevalue{\tr(\rho_{d}V)}^{2}+\tr(\rho_{d}V\rho_{d}V^{\dagger})-\tr(\rho_{d}^{(2)}(V\otimes
V^{\dagger}))$ (S49) $\displaystyle\leavevmode\nobreak\
=\big{|}\sum_{E}\bra{E}V\ket{E}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{2}\big{|}^{2}+\sum_{E,E^{\prime}}\bra{E}V\outerproduct{E^{\prime}}{E^{\prime}}V^{\dagger}\ket{E}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{2}\absolutevalue{\innerproduct{E^{\prime}}{\Psi_{0}}}^{2}-\sum_{E}\absolutevalue{\bra{E}V\ket{E}}^{2}\absolutevalue{\innerproduct{E}{\Psi_{0}}}^{4}$
$\displaystyle\overset{\text{ETH}}{\approx}\absolutevalue{v(\bar{E})}^{2}+\langle
V^{2}\rangle_{\text{micro}}\exp(-c^{\prime}N)\leavevmode\nobreak\ ,$ (S50)

for some constant $c^{\prime}>0$.

Meanwhile, to analyze
$\mathbb{E}_{\tau}\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]$, we
approximate
$\absolutevalue{\bra{z}\rho_{d}V_{d}\ket{z}}=\sum_{E}\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{2}\absolutevalue{\bra{E}V\ket{E}}\overset{\text{ETH}}{\approx}\overline{\bra{E}V\ket{E}}_{\bar{E}}\sum_{E}\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{2}=|v(\bar{E})|p_{\text{avg}}(z)$.
Then

$\displaystyle\mathbb{E}_{\tau}\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]$
$\displaystyle\leavevmode\nobreak\
=\sum_{z}\frac{\bra{z}\rho_{d}V_{d}\outerproduct{z}{z}\rho_{d}V^{\dagger}_{d}\ket{z}}{p_{\text{avg}}(z)}-\sum_{z}\frac{\bra{z,z}\rho^{(2)}_{d}\left(V_{d}\otimes
V_{d}^{\dagger}\right)\ket{z,z}}{p_{\text{avg}}(z)}\leavevmode\nobreak\ ,$
(S51) $\displaystyle\leavevmode\nobreak\
=\sum_{z}\frac{1}{p_{\text{avg}}(z)}\Big{[}|\sum_{E}\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{2}\bra{E}V\ket{E}|^{2}-\sum_{E}\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{4}\absolutevalue{\bra{E}V\ket{E}}^{2}\Big{]}$
(S52)
$\displaystyle\overset{\text{ETH}}{\approx}\sum_{z}\frac{\absolutevalue{v(\bar{E})}^{2}}{p_{\text{avg}}(z)}\Big{[}\Big{(}\sum_{E}\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{2}\Big{)}^{2}-\sum_{E}\absolutevalue{\innerproduct{z}{E}\innerproduct{E}{\Psi_{0}}}^{4}\Big{]}$
(S53) $\displaystyle\leavevmode\nobreak\
=\absolutevalue{v(\bar{E})}^{2}(1-D_{\beta}^{-1})$ (S54)

Since $D_{\beta}^{-1}$ is also typically exponentially small in system size,
we conclude that
$\absolutevalue{\mathbb{E}_{\tau}\mathbb{E}_{t}\left[F_{\text{XEB,d}}(t,\tau)\right]-\mathbb{E}_{t}\left[F(t)\right]}=O(\exp[-dN])$,
for some $d>0$. This is our desired result, obtained here through the ETH.

####  C.2.2 Extensions beyond local errors

The above argument shows that in the presence of local coherent errors,
$F_{d}$ faithfully estimates $F$ in the long time limit. This can be
immediately extended to incoherent errors by averaging over coherent errors.
By unraveling an error channel into quantum trajectories, these results also
apply to stochastic, local error channels, as long as the error rate is
sufficiently low (numerically investigated in Fig. S4).

Finally, we argue that these results also apply to coherent Hamiltonian
errors, i.e. quench evolution by different Hamiltonian parameter, utilized in
our process benchmarking examples. Time evolution under a slightly perturbed
local Hamiltonian can be expanded into a coherent sum of paths representing
evolution by the original Hamiltonian, weakly interspersed by the perturbation
term. This is similar to a stochastic error channel, except we can have
coherent interference between different paths. We argue that this interference
effect from off-diagonal terms is negligible: different paths have vanishing
overlap and essentially uncorrelated speckle patterns, and their sum is
further suppressed by their random phases. An analogous argument applies to
target state benchmarking within a ground state phase diagram [Fig. S9], which
can be understood as applying imaginary-time evolution (cooling) between
states. This effect is perturbative and can also be decomposed into paths of
local excitations, as long as the two states are sufficiently nearby.

###  C.3 Variance over time

The variance over time utilizes the fourth Hamiltonian twirling identity. The
proof of Theorem 2 estimates this as $O(\tr[\rho_{d}^{2}])$. We explicitly
compute the prefactors of the expected leading order contribution arising from
the simplest energy-conserving terms, illustrated in Fig. S3(b):

$\displaystyle\text{var}_{t}\left[F_{\text{id,d}}\right]$
$\displaystyle=\mathbb{E}_{t}\left[F_{\text{id,d}}(t)^{2}\right]-\mathbb{E}_{t}\left[F_{\text{id,d}}(t)\right]^{2}=16\left[\cancelto{0}{\tr(\rho_{d}^{2})-\tr(\rho_{d}^{(2)})}\right]+4\sum_{z,z^{\prime}}\frac{\absolutevalue{\bra{z}\rho_{d}\ket{z^{\prime}}}^{4}}{p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}})}+\cdots\leavevmode\nobreak\
.$ (S55)

Figure S3: Diagrammatic derivation of terms that contribute to (a)
$\mathbb{E}_{t}\left[F_{\text{id,d}}(t)\right]$ and (b)
$\text{var}_{t}\left[F_{\text{id,d}}(t)\right]$. In (a), the two permutations
and collision term give $2-D_{\beta}^{-1}$. In (b), we list the (expected)
leading order contributions. 16 out of the 24 possible permutations each give
$\tr(\rho_{d}^{2})$, these are cancelled out by corresponding collision terms:
representative permutations and collision terms are illustrated here. Another
four permutations give the leading order contribution to
$\text{var}_{t}\left[F_{\text{id,d}}(t)\right]=4\sum_{z,z^{\prime}}\absolutevalue{\bra{z}\rho_{d}\ket{z^{\prime}}}^{4}/(p_{\text{avg}}(z)p_{\text{avg}}({z^{\prime}}))\leq
4\tr(\rho_{d}^{2})$. Figure S4: Performance metrics for $F_{d}$ under single
and continuous errors. In all cases, our numerical results quantitatively
agree with our theoretical predictions (black dashed lines), with all
parameters numerically extracted from independent numerical data. (a,b) After
a single error, the (a) systematic error $\delta_{\text{sys.}}$ and (b)
temporal standard deviation $\delta_{\text{temp.}}$ quantitatively agrees with
the ansatz $\delta_{\text{sys.}}=(1-F)D_{\beta}^{-1}$ and
$\delta_{\text{temp.}}=\sqrt{5/4}(1-F)\text{std}_{t}\left[F_{\text{id,d}}\right]$,
across all four systems studied: a Bose-Hubbard model (red), Fermi-Hubbard
model (blue), trapped-ion spin chain Hamiltonian (purple), and 2D Rydberg
model (orange). (c,d) In the presence of continuous errors, $F$ decays with
time from 1 to 0 (inset). Therefore, we can test our predictions for the
dependence of our quantities of interest—the sample complexity (c), and the
temporal fluctuations (d)—on the fidelity $F$ in a single time series of open
system dynamics. We present such a time series for a Bose-Hubbard model with 8
bosons on 8 sites, where we plot the (c) sample complexity at each point in
time and the (d) temporal fluctuations in a narrow time window against the
fidelity $F$ at that point in time. After a short thermalization time, the
sample complexity follows our universal prediction in Eq. (S60):
$M\delta_{\text{stat.}}^{2}=(1+2F-F^{2})$, and the temporal fluctuations
satisfy Eq. (S57):
$\delta_{\text{temp.}}=(1-F)\text{std}_{t}\left[F_{\text{id,d}}\right]/2$. (e)
Under continuous errors, $F_{d}(t)$ will appear to lag behind $F(t)$ with a
delay time. In particular, the rescaled curves of $F_{d}/F$ for a range of
error rates $\gamma=0.0005$ to $\gamma=0.005$ collapse, and are close to our
prediction in Eq. (S63), obtained from the single-error response curve
(inset). (f) At low error rates, $F_{d}$ estimates $F$ well, but they deviate
significantly at high error rates. Inset: Across all platforms, a normalized
cost function $C(\gamma)/C(\infty)$ (S66) shows a crossover when the average
time between errors is less than an error autocorrelation time $t^{*}$ (S65).

While this expression allows us to theoretically understand the expected the
magnitude of $\text{var}_{t}\left[F_{\text{id,d}}\right]$, in practice we
numerically evaluate $\text{var}_{t}\left[F_{\text{id,d}}\right]$ and
investigate its relation to $\text{var}_{t}\left[F_{d}\right]$. In the
presence of errors, the variance of $F_{\text{XEB,d}}$ can be calculated using
Eq. (S36). This gives $\text{var}_{t}\left[F_{\text{XEB,d}}\right]\approx
F^{2}\text{var}_{t}\left[F_{\text{id,d}}\right]$ for incoherent errors. In the
case of a single, coherent error, we must include the coherence of
$\tilde{p}_{\perp}^{2}$: $\langle\tilde{p}_{\perp}^{2}\rangle\approx 2$ and we
obtain
$\text{var}_{t}\left[F_{\text{XEB,d}}\right]\approx\left[F^{2}+(1-F)^{2}\right]\text{var}_{t}\left[F_{\text{id,d}}\right]$
instead. In both cases, we also have
$\text{cov}_{t}\left[F_{\text{XEB,d}},F_{\text{id,d}}\right]\approx
F\text{var}_{t}\left[F_{\text{id,d}}\right]$.

Performing the same Taylor expansion as in Eq. (S40),

$\displaystyle\text{var}\left[\frac{X}{Y}\right]\approx\frac{\mathbb{E}[X]^{2}}{\mathbb{E}[Y]^{2}}\left[\frac{\text{var}[X]}{\mathbb{E}[X]^{2}}+\frac{\text{var}[Y]}{\mathbb{E}[Y]^{2}}-2\frac{\text{cov}[X,Y]}{\mathbb{E}[X]\mathbb{E}[Y]}\right]\,$
(S56)

gives the estimate

$\displaystyle\delta^{2}_{\textrm{temp}}=\text{var}_{t}\left[F_{d}\right]$
$\displaystyle\approx\frac{1}{4}(1-F)^{2}\text{var}_{t}\left[F_{\text{id,d}}\right]\leavevmode\nobreak\
\text{(open system dynamics)}$ (S57)
$\displaystyle\approx\frac{5}{4}(1-F)^{2}\text{var}_{t}\left[F_{\text{id,d}}\right]\leavevmode\nobreak\
\text{(coherent errors)}$ (S58)

Here, the larger fluctuations for a single coherent error are associated with
the fact that it is pure. We numerically evaluate
$\text{var}_{t}\left[F_{\text{id,d}}\right]^{-1}$, and compare our predictions
for $\text{var}_{t}\left[F_{d}\right]$ in terms of
$\text{var}_{t}\left[F_{\text{id,d}}\right]$ and $F$ in Fig. 3(c) of the main
text and Fig. S4(b,d) (respectively single error and open system dynamics).

###  C.4 Sample complexity

The sample complexity of $F_{d}$ is related to the variance of the estimator
$\hat{F}_{d}$ with a finite number of samples $M$. The variance
$\delta_{\text{stat}}^{2}$ is given by

$\delta_{\text{stat}}^{2}\equiv\text{var}\big{[}\hat{F}_{d}\big{]}=\frac{4\left(\sum_{z}q(z)\tilde{p}(z)^{2}-\left[\sum_{z}q(z)\tilde{p}(z)\right]^{2}\right)}{M\left[\sum_{z}p_{\text{avg}}(z)\tilde{p}(z)^{2}\right]^{2}}\overset{\text{large
}t}{\approx}\frac{1}{M}\left(\sum_{z}q(z)\tilde{p}(z)^{2}-\left[\sum_{z}q(z)\tilde{p}(z)\right]^{2}\right)$
(S59)

We estimate this by computing computing its time-averaged value, with the
third Hamiltonian twirling identity, which gives
$\mathbb{E}_{t}\left[\sum_{z}p_{\text{avg}}(z)\tilde{p}(z,t)^{3}\right]=3!-(2!)^{2}+\cdots\approx
2$. Our ansatz for $\tilde{q}(z,t)$ (S36) gives:

$\displaystyle\mathbb{E}_{t}\left[M\delta^{2}_{\text{stat}}\right]=\mathbb{E}_{t}\left[\sum_{z}q(z,t)\tilde{p}(z,t)^{2}-\left(\sum_{z}q(z,t)\tilde{p}(z,t)\right)^{2}\right]\approx
1+2F-F^{2}\leavevmode\nobreak\ .$ (S60)

This non-trivial prediction is verified in Fig. S4(c) and Fig. 3(d) of the
main text.

###  C.5 Moderate error rate: Effects of delay time $\tau_{d}$

Here we argue that under open system dynamics with a moderate error rate, the
ratio $F_{d}/F$ approaches a constant $F_{d}/F\approx\exp(\gamma\tau_{d})$,
where $\gamma$ is the total error-rate and $\tau_{d}$ is a response time
independent of system size. When the fidelity $F$ shows a simple exponential
decay, the delay time leads to $F_{d}$ lagging behind $F$ by a constant time
$\tau_{d}$. This is in contrast to the constant difference $F_{d}-F$ seen
under a single error, which is exponentially small in system size. We quantify
this delay time effect by introducing the notion of a response function of
$F_{d}$ after a single error, $F_{d}^{(1)}$ , and making the assumption that
the effects of multiple errors on $F_{d}$ are independent and factorize.

We assume that an error operator $X$ can occur at any time during the open
system dynamics. Without loss of generality, we normalize $X$ so it does not
change the norm of a typical state:
$\bra{\psi(t)}X^{\dagger}X\ket{\psi(t)}\approx 1$. The fidelity $F$ as a
function of time is:

$\displaystyle F(t)$ $\displaystyle\approx
p_{ne}(t)+\int_{0}^{t}dt_{1}p_{e}(t_{1})\absolutevalue{\langle
X(t_{1})\rangle}^{2}+\int_{0}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}p_{e}(t_{1},t_{2})\absolutevalue{\langle
X(t_{2})X(t_{1})\rangle}^{2}+\cdots$ $\displaystyle\approx
p_{ne}+\sum_{n=1}^{\infty}\frac{1}{n!}\left[\int_{0}^{t}dt_{i}p_{e}(t_{i})\absolutevalue{\langle
X(t_{i})\rangle}^{2}\right]^{n}=\exp[\left(F_{\beta}-1\right)\gamma
t]\leavevmode\nobreak\ ,$ (S61)

where we denote $p_{ne}(t)\propto\exp(-\gamma t)$ as the probability of not
having an error from time $0$ to $t$ and
$p_{e}(t_{1},\cdots,t_{n})\propto\gamma^{n}\exp(-\gamma t)$ as the probability
density of errors occurring at times $t_{1},\cdots,t_{n}$. If $n$ errors occur
during a trajectory, the pure state asspcoated with this trajectory will have
a fidelity $\absolutevalue{\bra{\psi(t_{n})}X\cdots
Xe^{-iH(t_{2}-t_{1})}X\ket{\psi(t_{1})}}^{2}$. As long as the error rate is
low enough that errors are typically far apart in time, we assume that this
expression factorizes: $F=\absolutevalue{\langle X(t_{n})\cdots
X(t_{1})\rangle}^{2}\approx\absolutevalue{\langle X\rangle}^{2n}_{\beta}\equiv
F_{\beta}^{n}$. Here the quantity $F_{\beta}\equiv\absolutevalue{\langle
X\rangle}^{2}_{\beta}$ should be interpreted as the fidelity after a single
error $X$ on a typical thermal state.

Meanwhile, we can estimate how $F_{d}$ behaves under open system dynamics in
terms of its behavior in response to a single error. We denote
$F^{(1)}_{d}(\tau)$ as the average $F_{d}$ value at time $\tau$ after a single
error.

$\displaystyle F_{d}(t)$ $\displaystyle\approx
p_{ne}(t)+\int_{0}^{t}dt_{1}p_{e}(t_{1})F_{d}(t-t_{1})+\int_{0}^{t}\int_{t_{1}}^{t}dt_{1}dt_{2}p_{e}(t_{1},t_{2})F_{d}(t-t_{2},t-t_{1})+\cdots$
$\displaystyle\approx
p_{ne}+\sum_{n=1}^{\infty}\frac{1}{n!}\left[\int_{0}^{t}dt_{i}p_{e}(t_{i})F_{d}^{(1)}(t-t_{i})\right]^{n}=\exp[\gamma\int_{0}^{t}d\tau\left(F_{d}^{(1)}(\tau)-1\right)]$
(S62)

where we again assume that the value of $F_{d}$ after $n$ errors factorizes
as: $F_{d}(t_{1},\cdots,t_{n})\approx\prod_{i}F_{d}^{(1)}(t-t_{i})$.

Combining Eqs. (S61), (S62), we have a simple expression for the ratios:

$\frac{F_{d}(t)}{F(t)}\approx\exp[\gamma\int_{0}^{t}d\tau\left(F_{d}^{(1)}(\tau)-F_{\beta}\right)]\approx\exp[\gamma\tau_{d}]\leavevmode\nobreak\
,$ (S63)

where the quantity
$\tau_{d}\equiv\int_{0}^{\infty}d\tau\left(F_{d}^{(1)}(\tau)-F_{\beta}\right)dt$
has a natural interpretation as the total area under the response curve
$F_{d}^{(1)}(\tau)-F_{\beta}$ (Fig. S4(e) inset) — if we ignore the
(exponentially small) systematic error $\bar{F}_{d}^{(1)}-F$, this total area
is a constant. This prediction is supported in Fig. S4(e), where we observe
the ratio $F_{d}(t)/F$ saturating at our predicted value for
$\exp[\gamma\tau_{d}]$, with $\tau_{d}$ extracted from the corresponding
single error curve $F_{d}^{(1)}(t)$.

###  C.6 Large error rate: Critical error rate and operator scrambling time

As the error rate increases, there is a critical rate at which $F_{d}$ shows
significant deviations from $F$ [Fig. S4(f) inset]. Naively, we would expect
that this error rate is solely determined by the delay time discussed above.
However, this is not necessarily the case. The delay time is a property of the
relaxation dynamics of a local perturbation occurring after the state has
reached thermal equilibrium. When the error rate is substantially larger than
the local relaxation time, the fidelity becomes very small even before the
state has a chance to reach thermal equilibrium. Therefore, the delay time
might not be the appropriate timescale for the critical error rate. Instead,
we seek to characterize a different timescale relevant in this non-equilibrium
dynamics. We hypothesize that a suitable timescale is the error operator
autocorrelation time — intuitively an operator scrambling time — $t^{*}$ with
respect to the nonequilibrium initial state. Our numerical results support our
hypothesis: we find that the critical error rate is related to $t^{*}$.

We define the autocorrelation function of the error operator $X_{j}$, for a
given initial state $\ket{\Psi_{0}}$ and quench Hamiltonian $H$, as:

$A_{j}(t)=\frac{\absolutevalue{\bra{\Psi_{0}}e^{iHt}X_{j}e^{-iHt}X_{j}\ket{\Psi_{0}}}^{2}}{\bra{\Psi_{0}}e^{iHt}X_{j}^{2}e^{-iHt}\outerproduct{\Psi_{0}}{\Psi_{0}}X_{j}^{2}\ket{\Psi_{0}}}\leavevmode\nobreak\
,$ (S64)

where we have chosen this normalization such that $0\leq A(t)\leq 1$. The
states $X_{j}\ket{\Psi(t)}$ are generally not normalized because the error
operators $X_{j}$ may not be unitary (for example, in the Bose-Hubbard model,
the errors are associated with photon scattering, and we take $X_{j}$ to be
$n_{j}$). We define the autocorrelation time as:

$t^{*}_{j}\equiv\frac{\int_{0}^{\infty}dt\absolutevalue{A_{j}(t)-A_{j}(\infty)}}{1-A_{j}(\infty)}\leavevmode\nobreak\
,$ (S65)

where $A_{j}(\infty)$ is the infinite time value
$A_{j}(\infty)=\lim_{t\rightarrow\infty}A_{j}(t)$. This can be analytically
evaluated from the second Hamiltonian twirling identity [Eq. (S4)]. We take
the overall operator scrambling time as the average of $t^{*}_{j}$ over error
operators $X_{j}$. This timescale is independent of system size, unlike other
timescales such as the time for the entanglement entropy to saturate, which
scales with system size.

Importantly, we find that the critical error rate is associated with the non-
equilibrium dynamics of the initial state, as opposed to the dynamics after
the state has reached thermal equilibrium. To numerically demonstrate this, we
initialize the noisy dynamics on a time-evolved state $e^{-iHt}\ket{\Psi_{0}}$
for sufficiently large $t$, after the state has thermalized. We find that in
the regime of error rates studied, there is no critical change in the accuracy
of $F_{d}$, merely a constant delay time between $F_{d}$ and $F$, as predicted
in Section C.5, which translates to a smooth increase in the relative error
$F_{d}/F$. This is analogous to Ref. [43, 42], in which the error rate in
random unitary circuits was found to have the largest effect at the start and
the end of the evolution.

We define a cost function

$C(\gamma)\equiv\int_{0}^{T}dt\absolutevalue{F_{d}(t)-F(t)}\leavevmode\nobreak\
,$ (S66)

that captures the total deviation between $F_{d}$ and $F$ over time. In Fig.
S4(f), we find that the rescaled cost function $C(\gamma)/C(\infty)$ shows a
crossover from 0 to 1 at $\gamma t^{*}\sim 1$, indicating that $t^{*}$ is the
relevant timescale for our setting. We choose such a normalization such that
the curves from different simulator platforms can be plotted on the same
graph: they demonstrate a cross over in a common regime $\gamma t^{*}\sim 1$.

###  C.7 Symmetry Resolution

In the main text, we described how $F_{d}$ naturally incorporates symmetries
of the quench Hamiltonian. Here, we discuss this in more detail and provide an
explicit example. If the initial state belongs to a single symmetry sector
$s$, the quenched state will remain in that sector. Assume for simplicity that
the state is Haar random in $s$. The simpler formulae $F_{\text{XEB}}$ and
$F_{c}$ would be appropriate if our measurements $z$ were also states in $s$.
However, this is generically not the case — if the outcome $z$ does not have a
well-defined symmetry quantum number, it should be weighted by its overlap
with $s$. The factor $p_{\text{avg}}(z)$ automatically weights $z$
appropriately, without explicit knowledge of the symmetry.

As an illustration, consider spatial inversion symmetry $\mathcal{I}$ on a
lattice. The even parity sector has basis states
$\ket{e_{z}}=\mathcal{N}(\ket{z}+\ket{\bar{z}})$, where $\bar{z}$ is the
mirrored configuration of $z$. If $\bar{z}=z$, $\mathcal{N}=1/2$, otherwise
$\mathcal{N}=1/\sqrt{2}$. Consider the optimal case where $\ket{\Psi_{0}}$ is
a random even-parity state and the dynamics is a random unitary that conserves
parity. The universal statistics would be apparent if we could measure in the
$\\{\ket{e_{z}}\\}$ basis — a challenging task. Instead, outcome
configurations $z$ have unequal average probabilities:
$p_{\text{avg}}(z)=1/D_{s}$ if $z=\bar{z}$ and $p_{\text{avg}}(z)=1/(2D_{s})$
otherwise, where $D_{s}$ is the dimension of the symmetry sector. While $p(z)$
does not follow the Porter-Thomas distribution, $\tilde{p}(z)$ does.
Therefore, $F_{d}$ correctly estimates $F$, while $F_{c}$ unequally weights
some basis states. This argument generalizes to other symmetries such as
translational invariance. Since $p_{\text{avg}}(z)$ automatically accounts for
symmetries, our conclusions hold even in the presence of extensively many
symmetries, as in the case of integrable systems. Finally, in a special case
where the readout basis states have well-defined symmetry quantum number,
$F_{d}$ does not estimate $F$. This is discussed in Section F.7.

###  C.8 Alternative fidelity estimator $F_{e}$

In the following section, we use the random MPS model to derive an alternative
fidelity estimator

$F_{e}\equiv\frac{\sum_{z}q(z)\tilde{p}(z)-1}{\sum_{z}p(z)\tilde{p}(z)-1}=\frac{F_{\text{XEB,d}}}{F_{\text{id,d}}}\leavevmode\nobreak\
.$ (S67)

Here, we summarize some of its properties. While the temporal fluctuations and
sample complexity of $F_{e}$ are similar to those of $F_{d}$, the systematic
error (after a long time) is smaller.

(i) Systematic error. The analysis in Section C.2 indicates that the
systematic error after a single error
$\delta_{\text{sys}}=\mathbb{E}_{t}F_{d}(t)-F=O(D_{\beta}^{-1})$. Repeating
the same analysis for $F_{e}$ indicates that the systematic error vanishes for
$F_{e}$, to first order. This indicates that $F_{e}$ is less sensitive to
small values of $D_{\beta}$, and hence may perform better in quench dynamics
with low effective temperature. This is supported in Fig. S6(e,f). More
fundamentally, $F_{e}$ is more robust against deviations of the denominator
$\sum_{z}p_{\text{avg}}(z)\tilde{p}(z)^{2}$ from 2. For example, in the
presence of particle hole symmetry, this denominator can be 3. $F_{d}$ can be
appropriately modified to take this into account, but $F_{e}$ estimates the
fidelity without any modification [Fig. S6(a,b)].

(ii) Non-local errors. Our analysis based on random matrix product states
indicates that $F_{e}$ will estimate $F$ for highly non-local errors. This is
numerically supported in Fig. S7(f).

(iii) Short-time response. A tradeoff of $F_{e}$ is that it typically has a
longer response time than $F_{d}$ to learn the fidelity $F$. This is
demonstrated in Fig. S7(e) and Fig. S8(b).

##  Appendix D Detailed performance analysis for short quench dynamics

###  D.1 Overview: random matrix product states

Figure S5: Random matrix product state model. (a) We consider two states
evolved by a short quench. In one copy, an error has occurred, which spreads
to a small subsystem $A$, and no error has occured in the other copy. (b) A
toy model for the non-equilibrium state is the random MPS — a matrix product
state with a random $\chi d\times\chi d$ matrix at every site, with $\chi$ the
MPS bond dimension, and $d$ the local Hilbert space. The error $V$ is also
modeled with a random unitary $W$ and an operator $\exp(i\theta X)$. (c) The
expression for $F_{d}$, averaged over random unitaries is mapped to the
partition function of an effective ferromagnetic Ising spin chain, with Ising
degrees of freedom $\sigma_{j},\sigma_{W}\in\\{\mathbb{I},\mathcal{S}\\}$.
This partition function is given in Eq. (S69). (d) We predict regimes of
validity for our benchmarking formulae $F_{d},F_{e}$, and $F_{\text{XEB}}$,
depending on the error operator size $N_{A}$ and MPS bond dimension $\chi$. A
typical quench is always in the regime of validity of $F_{c}$ (red arrow).

The infinite-time averaging in Section C tells us that our fidelity estimators
work at late times. Our formulae are also designed to work at early times,
before the state has reached equilibrium. To capture this early time physics,
we introduce a toy model, the random matrix product state.

We model the state $\ket{\Psi_{0}(t)}$ as a matrix product state (MPS) with
bond dimension $\chi$ and local Hilbert space dimension $d$. The MPS is
constructed by choosing a $\chi d\times\chi d$ random unitary $U_{j}$ on each
site $j$, as illustrated in Fig. S5(b). At every site, a $d$ degree of freedom
is eliminated by contracting to the state $\ket{0}$. At the boundary sites, an
additional $\chi$ degree of freedom is removed by contracting to a state
$\ket{0}_{\chi}$ on the bond-dimension space. This provides a toy model for a
state with entanglement $S\sim\log(\chi)$.

We model the error operator $V$ as $V=W\exp(i\theta X)W^{\dagger}$, where $V$
is taken to be supported on some region $A$, and $W$ is a random unitary on
$\mathcal{H}_{A}$. $\exp(i\theta X)$ is an operator whose precise form will be
irrelevant. However, the parameter $\theta$ measures the strength of the
coherent error, and can be used to tune the fidelity between the perturbed and
initial states. The random unitary $W$ models the scrambling nature of
$V(\tau)$ while making the calculation tractable using techniques of averaging
over random unitaries [61].

This toy model captures what we believe to be the essential physics:
entanglement and operator growth. In particular, it is well known that both
entanglement entropy and operator size grow linearly with time. This
corresponds to:

$\displaystyle\chi\sim\exp(v_{S}t)\leavevmode\nobreak\ ,\leavevmode\nobreak\
N_{A}\sim vt\leavevmode\nobreak\ ,$ (S68)

where $v_{S}$ is an entanglement velocity, and $v$ is an operator spreading
velocity [62, 63].

Our calculations give us expressions for how $F_{\text{XEB}}\equiv
D^{-1}\sum_{z}q(z)p(z)-1$ and $F_{\text{id}}\equiv D^{-1}\sum_{z}p(z)^{2}-1$,
and the fidelity $F$ depend on $N_{A}$ and $\chi$. Because of the random
unitary nature, there are no systematic patterns in $p(z,t)$:
$p_{\text{avg}}(z)=1/D$. Therefore, in this case $F_{d}=F_{c}$,
$F_{\text{XEB,d}}=F_{\text{XEB}}$ and $F_{\text{id,d}}=F_{\text{id}}$.
Furthermore, we will also consider $F_{e}=F_{\text{XEB}}/F_{\text{id}}$. Our
main result is:

$\displaystyle\mathbb{E}_{W,\mathbf{U}}\left[F_{\text{XEB}}+1\right]\approx\mathbb{E}_{\mathbf{U}}\left[F_{\text{id}}+1\right]\left(F+\exp[-\frac{d-1}{d}\frac{N_{A}}{\chi}]\frac{1-F}{2}\right)\leavevmode\nobreak\
,$ (S69)

where $\mathbb{E}_{W,\mathbf{U}}$ denotes averaging over the random unitaries
$W$ and $\mathbf{U}=\\{U_{j}\\}$, and $\mathbb{E}_{\mathbf{U}}$ denotes
averaging over $\mathbf{U}$.

The second term describes a crossover between two regimes: $\chi\ll N_{A}$ and
$\chi\gg N_{A}$. In the former regime,
$\exp[-\frac{d-1}{d}\frac{N_{A}}{\chi}]\approx 0$, and
$F_{\text{XEB}}+1\approx\left(F_{\text{id}}+1\right)F$. Furthermore since
$\chi$ is small, both $F_{\text{id}}$ and $F_{\text{XEB}}$ are exponentially
large and

$F\approx\frac{F_{\text{XEB}}+1}{F_{\text{id}}+1}\approx\frac{F_{\text{XEB}}}{F_{\text{id}}}=F_{e}\leavevmode\nobreak\
.$ (S70)

In the latter regime, $\exp[-\frac{d-1}{d}\frac{N_{A}}{\chi}]\approx 1$, and
$F_{\text{XEB}}+1\approx\left(F_{\text{id}}+1\right)\frac{1+F}{2}$, and

$F\approx
2\frac{F_{\text{XEB}}+1}{F_{\text{id}}+1}-1=F_{d}\leavevmode\nobreak\ .$ (S71)

Both $F_{d}$ and
$F_{e}\equiv\left[\sum_{z}p_{\text{avg}}(z)\tilde{q}(z)\tilde{p}(z)-1\right]/\left[\sum_{z}p_{\text{avg}}(z)\tilde{p}(z)^{2}-1\right]$
accurately estimate $F$ at late times. $F_{d}$ is valid since $\chi\approx
D\gg N_{A}$. However $F_{e}$ is also valid since $F_{\text{id}}\approx 1$ and
so $F_{\text{XEB}}\approx F$. These results give the regime of validity of
$F_{d}$ and $F_{e}$ sketched in Fig. S5(d). Further properties of $F_{e}$ are
discussed in Sec. C.8.

###  D.2 Remarks

The random MPS model qualitatively captures the short time physics. However,
it is not quantitatively accurate, at least to first order in $D^{-1}$. The
simplest example of this is that the state norm has non-zero variance:
$\text{var}_{\mathbf{U}}\left[\absolutevalue{\innerproduct{\Psi_{\mathbf{U}}}{\Psi_{\mathbf{U}}}}^{2}\right]\neq
0$. As another example, $F_{\text{id}}$ is, on average, bigger than 1, while
$\mathbb{E}_{t}\left[F_{\text{id,d}}\right]\leq 1$ for Hamiltonian dynamics
[Eq. (S38)]. Furthermore, the randomness in the MPS does not capture the
systematic trends in finite temperature states that lead to non-trivial
$p_{\text{avg}}(z)$ factors. Therefore, the random MPS model cannot be a
qualitative model for non-equilibrium quench-evolved states at either finite
or infinite temperature. Nevertheless, the qualitative conclusions of the
random MPS model — in particular the regimes of validity of $F_{d}$ and
$F_{e}$ — are consistent with numerical results at finite temperature.

###  D.3 Derivation of result

We model the error operator $V(\tau)=We^{i\theta X}W^{\dagger}$, where $W$ is
a Haar-random unitary supported on $A$ [Fig. S5(b)]. This models the effect of
scrambling and makes the calculation tractable. We first integrate over $W$ (a
unitary two-design calculation [68, 61]), and obtain two terms corresponding
to the permutations of two elements: $\sigma_{W}=\mathbb{I}$ and
$\sigma_{W}=\mathcal{S}$. We then integrate over $U$ to obtain effective Ising
chains for each value of $\sigma_{W}$: the integration over unitaries $U_{j}$
leaves a degree of freedom on each sites with two possible states, again
corresponding to (local) $\mathbb{I}$ and $\mathcal{S}$ permutations [Fig.
S5(c)]. These degrees of freedom “interact” with their nearest neighbors
ferromagnetically: the lowest energy configurations are those with all sites
in the $\mathbb{I}$ state or all sites in the $\mathcal{S}$ state. This
ferromagnetic interaction is encoded in the transfer matrix (of the partition
function) $T_{0}$ (S75).

The presence of the error operator $V(\tau)$ modifies this effective Ising
chain when $\sigma_{W}=\mathcal{S}$: it introduces an effective magnetic field
on the region $A$ supporting the error.

Our quantities of interest can be expressed in terms of the partition
functions of these two Ising chains. We call these partition functions
$f_{X,\mathbb{I}}$ and $f_{X,\mathcal{S}}$, corresponding to the value of
$\sigma_{W}$. For example, the average value of $F_{\text{id}}+1$ over all
possible choices of $U_{j}$ is exactly the partition function
$f_{X,\mathbb{I}}$. The presence of an error changes the value of
$F_{\text{XEB}}+1$. $F_{\text{XEB}}+1$, again averaged over unitaries $U_{j}$
and $W$, can be expressed in terms of $f_{X,\mathbb{I}}$, $f_{X,\mathcal{S}}$,
and the (averaged) fidelity $F$.

Our calculation gives the expressions:

$\displaystyle\mathbb{E}_{W,\mathbf{U}}\left[F\right]=\frac{\absolutevalue{\tr
V}^{2}-1}{D_{A}^{2}-1}\mathbb{E}_{\mathbf{U}}\left[\absolutevalue{\innerproduct{\Psi_{\mathbf{U}}}{\Psi_{\mathbf{U}}}}^{2}\right]+\frac{\absolutevalue{D_{A}^{2}-\tr
V}^{2}}{D_{A}\left(D_{A}^{2}-1\right)}\mathbb{E}_{\mathbf{U}}\left[\tr(\rho_{A}^{2})\right]\approx\frac{\absolutevalue{\tr
V}^{2}-1}{D_{A}^{2}-1}\equiv\bar{F}\leavevmode\nobreak\ ,$ (S72)
$\displaystyle\mathbb{E}_{W,\mathbf{U}}\left[F_{\text{XEB}}+1\right]=f_{X,\mathbb{I}}\bar{F}+\frac{f_{X,\mathcal{S}}}{D_{A}}\left(1-\bar{F}\right)\leavevmode\nobreak\
,$ (S73)
$\displaystyle\mathbb{E}_{\mathbf{U}}\left[F_{\text{id}}+1\right]=f_{X,\mathbb{I}}\leavevmode\nobreak\
,$ (S74)

where the coefficients $f_{X,\mathbb{I}}$ and $f_{X,\mathcal{S}}$ are the
partition functions of the Ising spin models, defined by the transfer
matrices:

$\displaystyle T_{0}$
$\displaystyle=\frac{1}{\chi^{2}d^{2}-1}\begin{pmatrix}\chi^{2}d^{2}-d&\chi
d(d-1)\\\ \chi d(d-1)&\chi^{2}d^{2}-d\end{pmatrix}\leavevmode\nobreak\
,\leavevmode\nobreak\ T_{1}=\begin{pmatrix}1&0\\\
0&\sqrt{d}\end{pmatrix}T_{0}\begin{pmatrix}1&0\\\
0&\sqrt{d}\end{pmatrix}\leavevmode\nobreak\ .$ (S75)

These transfer matrices are obtained from the Weingarten coefficients that
arise from the integration over random unitaries [61]. The partition functions
are

$\displaystyle f_{X,\mathbb{I}}$
$\displaystyle\equiv\frac{\chi}{\chi+1}\frac{\chi d+d}{\chi
d+1}\begin{pmatrix}1&1\end{pmatrix}T_{0}^{N-1}\begin{pmatrix}1\\\
1\end{pmatrix}=2\frac{\chi}{\chi+1}\left(\frac{\chi d+d}{\chi
d+1}\right)^{N}\leavevmode\nobreak\ ,$ (S76) $\displaystyle f_{X,\mathcal{S}}$
$\displaystyle\equiv\frac{\chi}{\chi+1}\frac{\chi d+d}{\chi
d+1}\begin{pmatrix}1&1\end{pmatrix}T_{0}^{k-1}\begin{pmatrix}1&0\\\
0&\sqrt{d}\end{pmatrix}T_{1}^{N_{A}}\begin{pmatrix}1&0\\\
0&\sqrt{d}\end{pmatrix}T_{0}^{N-N_{A}-k}\begin{pmatrix}1\\\
1\end{pmatrix}\leavevmode\nobreak\ .$ (S77)

As a reminder, $f_{X,\mathbb{I}}$ is the partition function of the
ferromagnetic effective Ising model, while $f_{X,\mathcal{S}}$ is the
partition function of the same Ising spin chain with an additional (very
strong) effective magnetic field applied on region $A$, that effectively locks
all the Ising spins in $A$ to the swap state $\mathcal{S}$.

We can simplify $f_{X,\mathbb{I}}$ and $f_{X,\mathcal{S}}$ by noting that the
dominant contributions to $f_{X,\mathbb{I}}$ are the configurations with all
sites in the $\mathbb{I}$ state, and all sites in the $\mathcal{S}$ state.
Meanwhile, the dominant contribution to $f_{X,\mathcal{S}}$ is the
configuration with all sites in the $\mathcal{S}$ state. Furthermore, $T_{0}$
has eigenvalues $\frac{\chi d-d}{\chi d-1}$ and $\frac{\chi d+d}{\chi d+1}$,
while $T_{1}$ has eigenvalues $\frac{(\chi d)^{2}-d^{2}}{(\chi d)^{2}-1}$ and
$d$. Therefore, for large $N_{A}$,
$\frac{f_{X,\mathcal{S}}}{D_{A}f_{X,\mathbb{I}}}$ is proportional to the ratio
of their largest eigenvalues:

$\frac{1}{D_{A}}\frac{f_{X,\mathcal{S}}}{f_{X,\mathbb{I}}}\approx\frac{1}{2}\left(\frac{\chi
d+1}{\chi d+d}\right)^{N_{A}}\\\
\approx\frac{1}{2}\exp(-\frac{d-1}{d}\frac{N_{A}}{\chi})\leavevmode\nobreak\
.$ (S78)

where the factor of a half comes from the breaking of the global Ising
$\mathbb{I}\leftrightarrow\mathcal{S}$ symmetry.

This gives our main result Eq. (S69), which predicts the relationship between
the cross-entropy benchmark $F_{\text{XEB}}$, fidelity, and the parameters in
our model: the operator size $N_{A}$ and MPS bond dimension $\chi$.

##  Appendix E Details of numerical simulations

Below we summarize the details of the models simulated: their Hamiltonians,
initial states, error models, and Hamiltonian parameter values.

###  E.1 Bose-Hubbard model

We numerically simulate quench evolution under the one-dimensional (1D) Bose-
Hubbard Hamiltonian [64] in the presence of noise, which is often dominated by
the photon scattering from optical lattices. In the quantum trajectory
formalism [51], such noisy dynamics is described by stochastically applying
the particle number operator $\\{n_{j}\\}$ at a rate $\gamma$ for each site.
An ensemble of such pure state evolution represents a mixed state $\rho(t)$ in
open dynamics. The 1D Bose-Hubbard Hamiltonian is:

$H_{\text{BH}}=-\Omega\sum_{j=1}^{N-1}\left(b^{\dagger}_{j}b^{\mathstrut}_{j+1}+\text{h.c.}\right)+\frac{U}{2}\sum_{j=1}^{N}n_{j}(n_{j}-1)\leavevmode\nobreak\
,$ (S79)

where $b^{\dagger}_{j}$ is the bosonic creation operator at site $j$, and
$n_{j}=b^{\dagger}_{j}b^{\mathstrut}_{j}$ is the bosonic number operator at
site $j$. In Fig. 1 and 2, we use the parameters $(\Omega,U)=(1,2.87)$. To
increase the ergodicity of our quench, in Fig. 3,4, we use the parameters
$(\Omega,U)=(1,0.5)$, whose ground state is in the superfluid phase. We take
as initial state the “unity filling state” $\ket{1,1,\dots,1}$, which is the
ground state deep in the Mott insulating phase $U\gg\Omega$.

Our error operators are $\\{n_{j}\\}$, which model the effects of scattering
from optical lattice photons [65]. In Fig. 1(b), we use an overall error rate
of $\gamma=0.05$.

In Fig. 2(a), we show data for a system of 6 bosons on 6 sites, and in Fig.
1(b) and Fig. 2(b,c), we show data for a system with 10 bosons on 10 sites
(with total Hilbert space dimension 92378).

####  E.1.1 Disorder Averaging

In Fig. 3(b-d), we present disorder averaged results for the response curve
and the sample complexity of $F_{d}$ against time. Disorder averaging is done
to remove the large fluctuations present at short times in order to clearly
observe systematic behavior.

For simplicity, we choose the following ensemble of Bose-Hubbard models with
disordered hopping strengths:

$H_{\text{BH}}=-\sum_{j=1}^{N-1}\Omega_{j}\left(b^{\dagger}_{j}b^{\mathstrut}_{j+1}+\text{h.c.}\right)+\frac{U}{2}\sum_{j=1}^{N}n_{j}(n_{j}-1)\leavevmode\nobreak\
,$ (S80)

with $\Omega_{j}=0.2m_{j}$, for $m_{j}$ a random integer from 1 to 10. As in
Eq. (S80), we choose $U=0.5$. We chose to disorder the hopping strengths in
order to observe different dynamics at short terms. If, for example, we had
chosen a disordered potential $\Delta_{j}n_{j}$ instead, due to the
homegeneous initial state $\ket{1,\dots,1}$, the resultant states
$\ket{\psi(t)}$ would not be significantly different from each other across
disorder realizations. Lastly, an error at $n_{\lfloor L/2\rfloor}$ is applied
at a fixed, early time $t=5$. This is because the error operator $n_{j}$ acts
trivially on the initial state $\ket{1,\dots,1}$.

###  E.2 Fermi-Hubbard model

The 1D Fermi-Hubbard model Hamiltonian is:

$H_{\text{FH}}=\sum_{j=1}^{N}\left[-\Omega\sum_{\sigma=\uparrow,\downarrow}\left(c^{\dagger}_{j,\sigma}c^{\mathstrut}_{j+1,\sigma}+\text{h.c.}\right)+Un_{j,\uparrow}n_{j,\downarrow}\right]\leavevmode\nobreak\
,$ (S81)

where $c^{\dagger}_{j,\sigma}$ is the creation operator of a fermion at site
$j$, with spin $\sigma\in\\{\uparrow,\downarrow\\}$, and
$n_{j,\sigma}=c^{\dagger}_{j,\sigma}c^{\mathstrut}_{j,\sigma}\in\\{0,1\\}$ is
the fermion occupation number. In Fig. 2 we use the parameter value
$(\Omega,U)=(1,3)$, while in other figures we use$(\Omega,U)=(1,1)$. For any
value of $U/\Omega$, this model is integrable by Bethe Ansatz [50], but we
numerically observe that it still passes the $2$nd no-resonance condition,
after disregarding degeneracies [Eq. (S3)]. As discussed in Section B,
degeneracies can be disregarded for the purposes of Theorem 2.

We take as initial state the anti-ferromagnetic state at half filling,
$\ket{\uparrow,\downarrow,\dots}$, and as error operators the total fermion
occupation $\\{n_{j,\uparrow}+n_{j,\downarrow}\\}$ at each site. As with the
Bose-Hubbard model, this models scattering from optical lattice photons [55].
In Fig. 1(c), we use a total error rate of $\gamma=0.2$, in a system with 5 up
and 5 down fermions on 10 sites (with total Hilbert space dimension 63504).

###  E.3 Trapped ion model

Our trapped ion model has Hamiltonian:

$H_{\text{TI}}=\sum_{j=1}^{N-1}\sum_{i>j}\left(\Omega/\absolutevalue{i-j}^{\alpha}\right)\sigma^{x}_{i}\sigma^{x}_{j}+h_{z}\sum_{j=1}^{N}\sigma^{z}_{j}\leavevmode\nobreak\
,$ (S82)

with $\sigma^{x}_{j}$ and $\sigma^{z}_{j}$ the Pauli $x$\- and $z$-operators
on site $j$. $H_{\text{TI}}$ is adapted from the experiment of Ref. [10]. Due
to the long-ranged interactions, following Ref. [10] we fix the “rescaled”
parameter value
$\tilde{h}_{z}=h_{z}/(N^{-1}\sum_{j=1}^{N-1}\sum_{i>j}1/\absolutevalue{i-j}^{\alpha})$
to ensure consistent behavior across system sizes, where $\alpha$ is the
exponent of the power-law interactions. As in Ref. [10], our initial state is
the all-downs computational basis state: $\ket{\downarrow,\dots,\downarrow}$.
We use the rescaled parameters $(\Omega,\tilde{h}_{z},\alpha)=(1,0.7,1)$. In
Fig. 1(d), we use on-site spin flip and dephasing errors:
$\\{\sigma^{x}_{j},\sigma^{z}_{j}\\}$, with a total error rate of
$\gamma=0.05$ in a system of 16 qubits (with a total Hilbert space dimension
of $2^{16}=65536$).

In Fig. 4(c) we learn the disorder strengths $\Delta_{j}$ of a disorder term
$H_{\text{dis}}=\sum_{j}\Delta_{j}\sigma^{z}_{j}$.

###  E.4 Rydberg blockaded model

We simulate a Rydberg blockaded model in two dimensions (2D), as recently
experimentally realized in Ref. [66]. Here each site can be either in the
ground $\ket{\downarrow}$ or Rydberg $\ket{\uparrow}$ state, subject to the
constraint that there are no two adjacent Rydberg excitations
$\ket{\uparrow}$. Its Hamiltonian is:

$H_{\text{PXP}}=\Omega\sum_{j=1}^{N}PX_{j}+\delta\sum_{j}n_{j}\leavevmode\nobreak\
,$ (S83)

with $PX_{j}\equiv\left[\bigotimes_{\langle
i,j\rangle}\outerproduct{\downarrow}{\downarrow}_{i}\right]\sigma^{x}_{j}$ the
projected spin flip operator — a spin flip on site $j$ is allowed only if all
the nearest neighbors $i$ of site $j$ are in state $\ket{\downarrow}$.

We use the all-downs state $\ket{\downarrow,\dots,\downarrow}$ as our initial
state. At zero detuning $\delta=0$, this would be an infinite temperature
state ($\langle H\rangle=\tr(H)/D$). In Fig. 2 we use the parameter values
$(\Omega,\delta)=(1,1.5)$ and in all other figures we use the parameters
$(\Omega,\delta)=(1,0.7)$ to tune $\ket{\downarrow,\dots,\downarrow}$ away