††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# Probing many-body quantum chaos with quantum simulators

Lata Kh Joshi  Center for Quantum Physics, University of Innsbruck, Innsbruck
A-6020, Austria  Institute for Quantum Optics and Quantum Information of the
Austrian Academy of Sciences, Innsbruck A-6020, Austria  Andreas Elben  Center
for Quantum Physics, University of Innsbruck, Innsbruck A-6020, Austria
Institute for Quantum Optics and Quantum Information of the Austrian Academy
of Sciences, Innsbruck A-6020, Austria  Institute for Quantum Information and
Matter and Walter Burke Institute for Theoretical Physics, California
Institute of Technology, Pasadena, CA 91125, USA  Amit Vikram  Joint Quantum
Institute, University of Maryland, College Park, MD 20742, USA  Condensed
Matter Theory Center, Department of Physics, University of Maryland, College
Park, MD 20742, USA  Benoît Vermersch  Center for Quantum Physics, University
of Innsbruck, Innsbruck A-6020, Austria  Institute for Quantum Optics and
Quantum Information of the Austrian Academy of Sciences, Innsbruck A-6020,
Austria  Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France  Victor
Galitski  Joint Quantum Institute, University of Maryland, College Park, MD
20742, USA  Condensed Matter Theory Center, Department of Physics, University
of Maryland, College Park, MD 20742, USA  Peter Zoller  Center for Quantum
Physics, University of Innsbruck, Innsbruck A-6020, Austria  Institute for
Quantum Optics and Quantum Information of the Austrian Academy of Sciences,
Innsbruck A-6020, Austria

###### Abstract

The spectral form factor (SFF), characterizing statistics of energy
eigenvalues, is a key diagnostic of many-body quantum chaos. In addition,
partial spectral form factors (PSFFs) can be defined which refer to subsystems
of the many-body system. They provide unique insights into energy eigenstate
statistics of many-body systems, as we show in an analysis on the basis of
random matrix theory and of the eigenstate thermalization hypothesis. We
propose a protocol that allows the measurement of the SFF and PSFFs in quantum
many-body spin models, within the framework of randomized measurements. Aimed
to probe dynamical properties of quantum many-body systems, our scheme employs
statistical correlations of local random operations which are applied at
different times in a single experiment. Our protocol provides a unified
testbed to probe many-body quantum chaotic behavior, thermalization and many-
body localization in closed quantum systems which we illustrate with numerical
simulations for Hamiltonian and Floquet many-body spin-systems.

##  I Synopsis

The ongoing development of quantum simulators provides us with unique
opportunities to study quantum chaos in many-body systems, and its connections
to random matrix theory (RMT) [1] and Eigenstate Thermalization Hypothesis
(ETH) [2, 3] in highly controlled laboratory settings. This refers to not only
the experimental realization of engineered Hamiltonian dynamics of isolated
quantum systems, which can be tuned from integrable to non-integrable, but
also the ability to measure novel observables beyond standard low-order
correlation functions [4, 5, 6, 7, 8]. It includes recent measurements of the
growth of entanglement entropies in quantum many-body systems [9, 10, 11, 12]
as well as of the decay of out-of-time-ordered correlation functions [13, 14,
15, 16, 17, 18, 19]. In this work, our interests lie in developing
experimentally feasible probes of universal RMT predictions for the statistics
of energy eigenvalues [20, 21, 22, 23, 24, 1] and predictions of the ETH for
the statistics of energy eigenstates [2, 3, 25, 26, 27, 28, 29] of quantum
chaotic many-body systems. Using these probes, we are further interested in
distinguishing many-body localized (MBL) systems [30, 31] from the chaotic
ones, where in the former the eigenvalue statistics are described by the
Poisson distribution [32, 33, 34] and the ETH is violated.

In this paper, we identify the spectral form factor (SFF), and its
generalization to partial SFF (PSFF), as observables of interest to reveal
energy level and eigenstate statistics. The SFF is defined in terms of the
time evolution operator of the quantum many-body system of interest and
provides us with statistics of energy levels [1]. The PSFF will be defined in
terms of the time evolution operator restricted to a subsystem of the many-
body system, and contains information on both, the statistics of energy
eigenvalues and energy eigenstates. We derive analytic expressions for the
PSFF in Wigner-Dyson random matrix ensembles. More generally, in chaotic
quantum many-body systems, the ETH imposes constraints on the statistics of
eigenstates, which are however typically violated in localized systems.
Therefore, the PSFF provides a direct probe of eigenstate thermalization and
localization.

The goal of the present work is to develop measurement protocols for the SFF
and PSFF in quantum spin models of arbitrary dimension, as realized for
instance with trapped ions [4, 6], Rydberg atoms [7] and superconducting
qubits [8]. We extend the randomized measurement toolbox [35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53] to infer the SFF and
PSFF from statistical correlations of local random operations applied at
different times in a single experiment. In contrast to the previous works
utilizing randomized measurements to infer properties of many-body quantum
states [35, 36, 37, 11, 39, 40, 41, 44, 42, 43, 51, 45, 52, 49, 48, 46, 50,
54, 55, 56, 57, 53] and (out-of-time-ordered) correlation functions of
Heisenberg operators [38, 18], the present protocol yields, with the SFF and
PSFF, genuine properties of the time evolution operator. We emphasize that the
present protocol is ancilla-free. This is in contrast to Ref. [58] where a
measurement scheme for the SFF was proposed requiring time evolution of an
extended system comprising of the quantum simulator and an auxiliary spin.

Our protocol to measure the PSFF and SFF in a quantum simulation experiment
can be readily implemented in existing experimental platforms. It requires
only to implement local (single-spin) random unitaries and projective
measurements, which have been previously demonstrated with high fidelity [11,
17, 18, 56]. Interestingly, in our protocol we obtain the SFF and PSFF from
the same experimental dataset. This enables an efficient scheme to test
universal RMT predictions for the energy eigenvalue spectrum and, at the same
time, to probe properties of the energy eigenstates and thermalization via
ETH.

We now turn to an overview of the main results of the paper. We start by
recalling the standard definition of the SFF, define the PSFF and describe
their estimation using randomized measurement protocol. We then illustrate the
key features of the (P)SFF and demonstrate our measurement protocol using an
example of a chaotic, periodically kicked spin$-1/2$ model. We will argue on
the basis of this example and show in later sections with detailed analytical
and numerical calculations that the SFF and PSFF provide unique insights into
the eigenvalue and eigenstate statistics of quantum many-body systems.

###  I.1 Spectral form factor

The SFF in a many-body quantum system with time-independent Hamiltonian $H$
and energy spectrum $\\{E_{j}\\}$ is defined as the Fourier transform of the
two-point correlator of the energy level density [1]. It can be expressed as

$\displaystyle\\!K(t)$
$\displaystyle\equiv\frac{1}{D^{2}}\,\overline{\sum_{i,j}e^{i(E_{i}-E_{j})t}}=\frac{1}{D^{2}}\,\overline{\textrm{Tr}\left[T(t)\right]\textrm{Tr}\left[T^{\dagger}(t)\right]}~{}.\\!\\!\\!$
(1)

Here, we normalize $K(t)$ such that
$K(0)=D^{-2}\textrm{Tr}\left[\mathbb{1}\right]^{2}=1$, with $D$ the Hilbert
space dimension and have defined the unitary time-evolution operator
$T(t)\equiv\exp(-iHt)$. The overline denotes a possible disorder or ensemble
average over an ensemble of $T(t)$, which is needed due to non-self-averaging
behavior of the SFF [59]. Replacing the energies $E_{i}$ with quasi-energies,
this definition carries over to Floquet models with time-periodic evolution
operator $T(t=n\tau)=V^{n}$ ($n\in\mathbb{N}$) and $V$ the Floquet time
evolution operator for a single period $\tau$ 111Denoting the set of
eigenvalues of the Floquet operator $V$ with $\\{\exp(-iE_{i}\tau)\\}$, the
quasi-energy eigenvalues $\\{E_{i}\\}$ are only defined up to multiples of the
driving frequencies $\omega=2\pi\tau^{-1}$. We fix them to lie in the interval
$[0,\omega]$. .

The SFF is a probe of the universal properties of the statistics of energy
eigenvalues in chaotic and localized systems. Lately, it has played a key role
in a variety of different fields, interconnecting quantum chaos [1], quantum
dynamics of black holes [61, 62, 63, 64], condensed matter systems [65, 66,
67, 68, 69, 70, 71, 72, 73, 74, 75], and the dynamics of thermalization [76].
In Fig. 1(a), we illustrate its behavior in the context of a periodically
kicked spin-$1/2$ system. The time evolution operator $T$ at integer multiples
$n\in\mathbb{N}$ of driving period $\tau$ is given by $T(t=n\tau)=V_{3}^{n}$
with,

$V_{3}=e^{-iH^{(x)}\tau/3}e^{-iH^{(y)}\tau/3}e^{-iH^{(z)}\tau/3}~{}.$ (2)

Here, the Hamiltonians $H^{(x,y,z)}$ contain nearest-neighbor interactions
with strength $J=3\tau^{-1}$ and disordered transverse fields with strength
$h_{i}^{(x,y,z)}\in[-J,J]$,

$\displaystyle
H^{(x,y,z)}=J\sum_{i=1}^{N-1}\sigma^{(x,y,z)}_{i}\sigma^{(x,y,z)}_{i+1}+\sum_{i=1}^{N}h_{i}^{(y,z,x)}\sigma^{(y,z,x)}_{i}~{},$

and $\sigma^{a}$ [$a\in(x,y,z)$] denote the Pauli matrices. We have denoted
the number of spins with $N$ such that $D=2^{N}$. An ensemble average is
naturally performed by averaging over many instances of
$T(t=n\tau)=V_{3}^{n}$, each with local disorder potentials $h_{i}^{(a)}$
sampled independently from the uniform distribution on $[-J,J]$.

Figure 1:  Illustration of the characteristic properties of the SFF and PSFF
using the chaotic spin-$1/2$ Floquet model $V_{3}$. (a) We display the SFF
$K(t)$ for the Floquet model $V_{3}$ with $N=6$ qubits as a function of time
$t$. We observe characteristic features such as the ramp between $t\sim\tau$
to $t=t_{H}=2^{N}\tau$ and a plateau for $t>t_{H}$. (b) For the PSFF
$K_{A}(t)$ we observe ramp, plateau and, in particular, a constant, additive
shift of the PSFF compared to the SFF, which depends on the subsystem size
$N_{A}$ of the subsystem $A$. We have chosen subsystems A from the middle of
the total system. In both, the colored lines show the numerically calculated
SFF and PSFFs, averaged over $8000$ disorder realizations. In addition, we
illustrate our measurement protocol (see Sec. I.3) by simulating $M=2\times
10^{5}$ experimental runs (single-shot randomized measurements) at each time
and display the estimated SFF and PSFF as black dots with associated error
bars. The dashed green line in panel (a) sketches the form of the SFF
generically expected in a many-body localized model.

As shown in Fig. 1(a), the SFF $K(t)$ for this model and choice of parameters
exhibits a period of linear growth, before transitioning to a constant at time
$t/\tau\approx D=2^{N}$. This ramp-plateau structure of the SFF is a
characteristic feature of quantum chaotic systems [77, 78, 1], originating
from (quasi-)energy level repulsion and spectral rigidity [77], and is
predicted by RMT [1, 24]. In particular, as we briefly review in App. A, RMT
for time evolution operators $T(t=\tau n)=V^{n}$, with $V$ from the circular
unitary ensemble (CUE), yields

$\displaystyle K(t)=\frac{1}{D}\begin{cases}{t}/{t_{H}},&0<t\leq t_{H}\\\
1,&t>t_{H}~{}.\end{cases}$ (3)

Here, the slope of the ramp and the onset of the plateau is determined by the
Heisenberg (or plateau) time $t_{\text{H}}$ which is connected to the mean
inverse spacing of adjacent (quasi-) energies. It typically scales with the
Hilbert space dimension $t_{\text{H}}/\tau\sim D$ — for $V$ from CUE,
$t_{H}/\tau=D$ [1, 64]. Thus, the SFF is expected to drop with increasing
Hilbert space dimension $D=2^{N}$, as $D^{-2}$ at times $1\lesssim t/\tau\ll
D$ and as $D^{-1}$ at times $t/\tau\gtrsim D$. Fig. 1(a) shows that the SFF
$K(t)$ for the $V_{3}$ model closely follows the CUE prediction after the
initial few time steps. This time after which the many-body model shows the
same SFF as the one in RMT is known as the Thouless time $t_{\rm{Th}}$ [65].
For the model $V_{3}$ we note that $t_{\rm{Th}}\approx 5\tau$ (see also Sec.
III). Therefore, the quasi-energy eigenvalues of the Floquet operator $V_{3}$
exhibit Wigner-Dyson statistics (see also Ref. [58]).

In contrast to the example of a chaotic system $V_{3}$ presented above, the
energy eigenvalues of integrable and localized models are known to exhibit
Poissonian statistics [32, 30, 33, 34]. This corresponds to a flat SFF without
a ramp which is, after an initial transient regime, constant in time [1],
$K(t\gg 0)=1/D$. This is schematically shown in Fig. 1(a) with green dashes.
These distinct features of the SFF have been pivotal in characterizing many-
body chaotic and MBL phases [69, 70, 58].

###  I.2 Partial Spectral Form Factor

The SFF reveals information on the statistics of (quasi-) energy eigenvalues.
It is however by definition insensitive to properties of the (quasi-) energy
eigenstates. In this subsection, we define the PSFF and illustrate its
essential properties connected to properties of eigenvalues and eigenstates.

For a fixed subsystem $A\subseteq\mathcal{S}$ of the total system
$\mathcal{S}$ with complement $B$ ($A\cup B=\mathcal{S}$) and Hilbert space
dimensions $D_{A}$ and $D_{B}$ respectively ($D=D_{A}D_{B}$), we define the
PSFF as

$\displaystyle K_{\mathrm{A}}(t)$
$\displaystyle\equiv\frac{1}{DD_{A}}\,\overline{\sum_{i,j}e^{i(E_{i}-E_{j})t}\textrm{Tr}_{B}\left[\rho_{B}(E_{i})\rho_{B}(E_{j})\right]}$
$\displaystyle=\frac{1}{DD_{A}}\,\overline{\textrm{Tr}_{B}\left[\textrm{Tr}_{A}\left[T(t)\right]\textrm{Tr}_{A}\left[T^{\dagger}(t)\right]\right]}~{},$
(4)

where $\rho_{B}(E_{i})=\textrm{Tr}_{A}\left[\ket{E_{i}}\bra{E_{i}}\right]$
denotes the reduced density matrix obtained after partial trace of the
eigenstate $\ket{E_{i}}$ of the Hamiltonian $H$ (the Floquet time evolution
operator $V$) with energy (quasi-energy) $E_{i}$. Here, the normalization of
$K_{A}(t)$ is chosen such that
$K_{A}(0)=\textrm{Tr}_{B}\left[\textrm{Tr}_{A}\left[\mathbb{1}\right]^{2}\right]/(DD_{A})=1$.
Hence, the SFF and PSFF coincide when $A=\mathcal{S}$, i.e.
$K_{\mathrm{A}=\mathcal{S}}(t)=K(t)$. We emphasize that for
$A\subset\mathcal{S}$, the PSFF $K_{A}(t)$ contains non-trivial contributions
from the eigenstates $\ket{E_{i}}$: We obtain terms of the form
$\mathrm{Tr}(\rho_{B}(E_{i})^{2})$ and
$\mathrm{Tr}(\rho_{B}(E_{i})\rho_{B}(E_{j}))$ ($i\neq j$) which correspond to
the purity and overlap of reduced eigenstates. As shown below, a measurement
of the PSFF allows to extract these purities and overlaps, averaged over
spectrum and ensemble, i.e. allows to characterize (second-order moments of)
the statistics of eigenstates.

We remark that $K_{A}(t)$ has been previously discussed as a topological
invariant in the classification of symmetry-protected matrix product unitaries
in Ref. [79]. Its limiting cases for special subsystems ($A$ or $B$ consisting
of a single site, in the limit of a large local Hilbert space dimension) have
been used to study matrix elements of local operators in the energy eigenbasis
in 1D Floquet circuits, with comparisons to random matrix predictions for
eigenstate statistics in these subsystems (as a special case of ETH) [80].

In this work, we identify a general shift-ramp-plateau structure of the PSFF,
which reveals a direct connection to ETH contained in the subsystem dependence
of the PSFF. In Fig. 1(b), we display the PSFF for the Floquet model (2) for
various subsystems $A$, where $N_{A}$ denotes number of qubits in the
subsystem such that $D_{A}=2^{N_{A}}$. We first note that the PSFF also has a
ramp and plateau, similar to the full SFF. The slope of the ramp is nearly
identical for the displayed subsystem sizes $N_{A}\gtrsim N/2$ $=3$, which
holds more generally for $D_{A}\gg 1$ in the CUE model, and the onset of the
plateau in the PSFF takes place at the Heisenberg time $t_{H}$. Crucially, we
find that, at late times comparable to the onset of the ramp, there is a
subsystem dependent additive shift of the PSFF $K_{A}(t)$ compared to the full
SFF $K(t)$.

Similar to the case of the full SFF, we can compare the behavior of the PSFF
to predictions of RMT. As detailed in Sec. II, we find that RMT yields for
time evolution operators $T(t=\tau n)=V^{n}$, with $V$ from the CUE, and
sufficiently large subsystems $A,B$, ($D_{A},{{D}_{B}}\gg 1$),

$\displaystyle
K_{A}(t)=\frac{1}{D_{A}^{2}}+\frac{1}{D}\begin{cases}{t}/{t_{H}},&0<t\leq
t_{H}\\\ 1,&t>t_{H}~{}.\end{cases}$ (5)

As shown in Fig. 1(b), and analyzed in detail by further numerical studies in
Sec. III, the PSFF (and SFF) for the $V_{3}$ model follows closely the RMT
predictions. This indicates that both (quasi-) energy eigenvalues and
eigenstates of $V_{3}$ exhibits the Wigner-Dyson statistics of the CUE. We
remark that this is consistent with previous works demonstrating that
(sub-)systems of chaotic Floquet systems thermalize to infinite temperature
states as per RMT [81, 82, 83, 33, 84, 80].

Partial spectral form factor and eigenstate thermalization hypothesis – Using
the example of a chaotic Floquet model, we have illustrated above the
essential features of the PSFF in chaotic quantum systems. In Sec. II, we
analyze its behavior in detail invoking subsystem ETH [29] for the reduced
eigenstates, which is a conjecture regarding the distribution of eigenstates
responsible for the thermal behavior (in the standard sense of ETH) of few-
body observables in chaotic systems.

By separating out the components of the reduced density matrix into maximally
mixed, smooth and fluctuating parts as a function of energy, a generic late
time expression for PSFF can be obtained. From here, we later conclude that
the features of the ramp, plateau and shift are generic features of the PSFF
in chaotic quantum many-body systems. These features are directly connected to
the spectrum and ensemble averages of the subsystem purities
$\mathrm{Tr}_{B}(\rho_{B}(E)^{2})$ and of the overlaps of reduced eigenstates
$\mathrm{Tr}_{B}(\rho_{B}(E_{i})\rho_{B}(E_{j}))$. Furthermore, the magnitudes
of these features in the chaotic systems follow specific constraints when the
eigenstates satisfy subsystem ETH, see Sec. II.2.2. In particular, we show
that this shift, connected to the average overlaps, enables the detection of
thermalization of eigenstates in the framework of subsystem ETH.

Let us take for instance the shift seen in the Fig. 1, defined precisely in
terms of the fluctuating part of the density matrix later in Sec II.2. For
chaotic models, the shift can be identified as the time independent constant
during the linear ramp phase, and for $D_{A}\ll D$ it is approximated by
$K_{A}(t_{0})-K(t_{0})$ where $t_{Th}<t_{0}\ll t_{H}$. If the eigenstates
follow ETH, it is expected that,

$K_{A}(t_{0})-K(t_{0})\approx O\left(\frac{1}{D_{A}^{2}}\right)~{}.$ (6)

This can be noted for the CUE in the Eqs. (3) and (5) as well as for the
$V_{3}$ model in Fig. 1, where the shift above SFF is seen to be increasing as
the $N_{A}$ decreases and is found to follow Eq. (6) (see Sec. III for more
numerical details). On the other hand, for eigenstates which do not
thermalize, the time independent shift above SFF is generically much larger
than $O(1/D_{A}^{2})$.

As illustrated above, the SFF and PSFF of a quantum many-body system provide
crucial insights into the statistics of energy eigenvalues and eigenstates,
which results in a joint observation of chaos and validity of ETH. The
question arises of how to probe the SFF and PSFF in today’s quantum devices.
In the next subsection, we present our measurement protocol which can be
directly implemented in state-of-the-art quantum simulation platforms
realizing lattice spin models. It builds on the toolbox of randomized
measurements.

###  I.3 Randomized measurements of spectral form factors

Initially, randomized measurements have been proposed and experimentally
implemented to characterize many-body quantum states [35, 36, 37, 11, 39, 40,
41, 44, 42, 43, 51, 45, 52, 49, 48, 46, 50, 54, 55, 56, 57, 53] and (out-of-
time-ordered) correlation functions of Heisenberg operators [38, 18].
Randomized measurements on quantum states exploit statistical correlations
obtained between measurements obtained from different random bases. However,
for measuring an object like the SFF, we need to access the full trace of the
time evolution operator $T(t)$, summing contributions from all its
eigenstates. Therefore, we need to devise a protocol that can measure how
various initial states are propagated via $T(t)$, in a way that allows to
extract the SFF from standard projective measurements. This subsection
provides this protocol and the estimation formulas to achieve this. We also
comment on statistical errors arising from a finite number of experimental
runs which are elaborated in detail in Sec. V.

####  I.3.1 Description of the protocol

Before describing the experimental sequence in detail, we first outline the
key idea of our protocol: As visualized in Fig. 2, we consider a system
$\mathcal{S}$ of $N$ qubits. The first step of our protocol is to prepare a
random product state of these qubits. Next, this state is evolved with $T(t)$.
Finally, a local measurement in the conjugate random product basis is
performed, in order to probe how the time-evolved state compares to the
initial random product state. This is repeated for many random product states
in order to sample the complete trace $\textrm{Tr}\left[T(t)\right]$ of the
time evolution operator and its adjoint uniformly. For instance, in the
trivial case $T(t=0)=\mathbb{1}$, we obtain that the ‘time-evolved’ state
always matches to the initial random state corresponding to
$D^{-1}\textrm{Tr}\left[T(0)\right]=1$. At later times $t$, we obtain in
general a more complex statistics of measurement results from which we can
extract the SFF and PSFF.

In our protocol, we note that the ensemble average over time evolution
operators in the definition of SFF and PSFF can be favorably combined with the
averaging over random product states and measurement bases. As detailed in the
prescription of the protocol in the next paragraph, each time evolution
operator can thus in practice be applied only to a single random initial
product state and measured only once in the corresponding randomized basis,
i.e., only a single-shot measurement for each time evolution operator is
sufficient in our protocol.

Figure 2: Probing SFF and PSFF using randomized measurements. We present our
protocol for the measurement of the SFF and PSFF using statistical
correlations of local random unitaries applied at different times in a single
experiment. We begin with a product state $\rho_{0}=\ket{0}\bra{0}^{\otimes
N}$. Before and after the time evolution $T(t)$, we apply random local
rotations $U=\bigotimes_{i}u_{i}$ and $U^{\dagger}$, respectively, where local
unitaries $u_{i}$ are sampled from a unitary $2-$design. Here, $T(t)$ can be
generated as Hamiltonian evolution, $T(t)={\rm{exp}}(-iHt)$, or Floquet
dynamics, $T(t=n\tau)=V^{n},~{}n\in\mathbb{N}$, where $V$ is Floquet evolution
operator for time period $\tau$. In the last step, a single-shot measurement
is performed in the $z-$basis to collect a bitstring of the form
$\mathbf{s}=(s_{1},s_{2},...,s_{N})$ with $s_{i}\in\\{0,1\\}$. This procedure
is repeated $M$ times and $M$ bitstrings are collected to estimate the SFF and
PSFF using Eqs. (7) and (8). The gray shaded region shows one possible choice
of the subsystem $A$.

In detail, the experimental recipe reads as follows: (i) We begin with a
product state $\rho_{0}=\ket{\mathbf{0}}\bra{\mathbf{0}}$ with
$\ket{\mathbf{0}}\equiv\ket{0}^{\otimes N}$. (ii) On this initial state, we
apply local random unitaries $U=\bigotimes_{i=1}^{N}u_{i}$ where $u_{i}$ are
the local unitaries independently sampled from a unitary 2-design [85, 86] on
the local Hilbert space $\mathbb{C}^{2}$. Here, unitary 2-designs are
ensembles of random unitaries whose first and second moments match the moments
of the Haar measure on the unitary group (defining the CUE) [85, 86]. Examples
of unitary 2-designs on $\mathbb{C}^{2}$ include the (discrete) single-qubit
Clifford group as well as uniformly distributed unitary $2\times 2$ matrices
which can be sampled for instance via the algorithm presented in Ref. [87].
(iii) We evolve the system in time, i.e. apply a time evolution operator
$T(t)$, which is generated by a Hamiltonian $H$ (or Floquet operator $V$) with
randomly sampled disorder potentials. (iv) We apply the adjoint local random
unitary $U^{\dagger}$ resulting in the final state
$\rho_{f}(t)=U^{\dagger}T(t)U\rho_{0}U^{\dagger}T^{\dagger}(t)U$. (v) Lastly,
we perform a single-shot measurement in the computational basis with outcome
bitstring $\mathbf{s}=(s_{1},\dots,s_{N})$ with $s_{i}\in\\{0,1\\}$ for
$i=1,\dots,N$. This concludes a single experimental run of our protocol. Steps
(i)-(v) are now repeated $M$ times with new disorder realizations and new
local random unitaries such that a set of outcome bitstrings
$\mathbf{s}^{(r)}$ with $r=1,\dots M$ is collected.

####  I.3.2 Estimation formulas and illustrations

The statistics of the measured bitstrings $\mathbf{s}^{(r)}$, $r=1,\dots M$,
depends on the applied time evolution operators $T(t)$. Using the theory of
unitary $2$\- designs, we can express the SFF as a function of this data. We
define

$\widehat{K(t)}=\frac{1}{M}\sum_{r=1}^{M}\;(-2)^{-|\mathbf{s}^{(r)}|}~{},$ (7)

where $|\mathbf{s}|\equiv\sum_{i}s_{i}$. As we show in Sec. IV,
$\widehat{K(t)}$ yields an (unbiased) estimate of the SFF for a finite number
$M$ of experimental runs and converges to $K(t)$ when $M\to\infty$.

Remarkably, from the same measurement data $\mathbf{s}^{(r)}$, we have also
access to the PSFF $K_{A}(t)$ for arbitrary subsystems $A\subseteq\mathcal{S}$
via post-processing. To this end, we simply project the measured bitrings on
the subsystem $A$ of interest, i.e., define $\mathbf{s}_{A}=(s_{i})_{i\in A}$,
and use

$\widehat{K_{A}(t)}=\frac{1}{M}\sum_{r=1}^{M}\;(-2)^{-|\mathbf{s}_{A}^{(r)}|}~{},$
(8)

which gives an (unbiased) estimate for $K_{A}(t)$ for finite $M$ and converges
to $K_{A}(t)$ when $M\to\infty$ (see Sec. IV).

In Fig. 1(a-b), we illustrate our measurement protocol in the context of the
periodically kicked spin-$1/2$ model $V_{3}$, Eq. (2). We consider a total
system size of $N=6$ qubits and present the simulated experimental results
(black dots and error bars) for $K(t)$ and $K_{A}(t)$ using $M=2\times 10^{5}$
experimental runs for the single-shot sequence shown in Fig. 2 at each time
$t$. We observe that the simulated experiment agrees with the exact numerical
calculations at all times $t$ within error bars. Here, error bars, indicating
the standard error of the mean, quantify statistical errors arising from the
finite measurement budget (i.e. the finite number $M$ of simulated single-shot
measurements), see next subsection.

####  I.3.3 Statistical errors and remarks

The SFF and PSFF can be accessed from the same set of measurement data via the
estimators defined in Eqs. (7) and (8). Statistical errors arise in practice
from a finite number $M$ of experimental runs, and are governed by the
variance of these estimators. We discuss statistical errors in detail via
numerical and analytical calculations in V, and find a typical scaling of
$M\sim 10^{N_{A}}\approx 2^{3.32N_{A}}$ to access the (P)SFF of a (sub-)system
of size $N_{A}$ up to a fixed relative error. Such exponential scaling of the
measurement effort reflects the exponential decrease of the SFF with system
size [see remarks below Eq. (3)]. We emphasize however that this scaling of
the experimental effort is substantially better than for quantum process
tomography which requires at least $\sim 2^{5N_{A}}$ experiments to
reconstruct the full time evolution operator $T(t)$ [88]. Importantly, and in
contrast to quantum process tomography, the initial state and the measurement
basis coincide in our protocol.

As detailed in Sec. V, we can further decrease the required number of
experimental runs to observe the ramp and plateau of the (P)SFF, by
considering an averaged PSFF. Here, an average over PSFFs of all subsystems
with a fixed size is performed. This results in a further improved signal-to-
noise ratio.

Lastly, we remark that our protocol shares some similarities with randomized
benchmarking [89, 90, 91, 92, 93], where however global random unitaries and
their inverses are applied sequentially. In the case of randomized
benchmarking the goal is to characterize noise and decoherence acting during
the implementation of these global random unitaries. In contrast, with our
protocol, the aim is to characterize a unitary time evolution operator $T(t)$
using local random unitaries $U=\bigotimes_{i}u_{i}$ applied before and after
$T(t)$, which can be prepared with high fidelity [11, 40].

##### Organization of the paper:

In the remainder of the manuscript, we elaborate on the contents of the above
synopsis with technical details, derivations, and examples. In Sec. II, we
provide an in-depth theoretical analysis of the PSFF in RMT and in generic
many-body models in relation to ETH. The analytic results are compared with
numerics in Sec. III where we consider many-body models undergoing Floquet and
Hamiltonian evolution. For the latter, we discuss both, chaotic and MBL
phases. Sec. IV contains the necessary background and proof of our protocol to
measure the SFF. In Sec. V, we discuss statistical errors, arising in our
measurement protocol from a finite number of experimental runs, and the
influence of experimental imperfections. Lastly, we summarize in Sec. VI with
some concluding remarks and future directions.

##  II Partial Spectral Form Factor: Analytic Results

In this section, we analyze the origin of the main features observed in the
PSFF, namely the ramp, plateau and shift, based on analytical calculations. We
provide arguments to show that the PSFF generically is a reliable probe of
eigenvalue correlations characterizing chaotic and localized phases, signified
by the presence and absence of a late time ramp-plateau structure
respectively. In addition, we show that the specific features observed in the
PSFF are related to the ensemble and spectrum averaged second-moments of
reduced density matrices of eigenstates at different energies, and therefore
provide a useful measure of eigenstate properties.

This section is organized as follows. In Sec. II.1, we analyze the PSFF in
standard Wigner-Dyson random matrix ensembles (see App. A for a brief
discussion), which are mathematically idealized models of quantum chaotic
systems in which the PSFF can be obtained exactly. These ensembles display the
essential features of the PSFF and present a clear example of the roles of
eigenvalue and eigenstate statistics in these features. This is followed by a
discussion of more general chaotic systems in Sec. II.2, where we show that
the PSFF detects thermalization in the sense of ETH [2, 3, 25, 26, 27, 28, 29]
in addition to level statistics (see also Ref. [80], that compares ETH for
Floquet circuits to random matrix ensembles using the PSFF for specific
subsystem sizes). We then discuss the PSFF in localized systems in Sec. II.3,
and summarize our main conclusions for all cases in Sec. II.4.

Common to all these cases is the fact that the time-independent part of the
PSFF in Eq. (I.2) is given by the plateau value, which depends only on the
eigenstate purities (assuming no degeneracies) i.e.
$K_{\mathrm{A}}(t\rightarrow\infty)={\mathcal{P}}_{B}/D_{A}$, where

$\mathcal{{P}}_{B}=\frac{1}{D}\overline{\sum_{i}\textrm{Tr}_{B}\left[\rho_{B}^{2}(E_{i})\right]}$
(9)

is the (spectrum- and ensemble-)averaged purity of the reduced energy
eigenstates. For later reference, we separate out this time-independent
plateau value,

$\displaystyle K_{\mathrm{A}}(t)=$
$\displaystyle\frac{\mathcal{P}_{B}}{D_{A}}+\frac{1}{DD_{A}}\sum_{i\neq
j}\overline{e^{i(E_{i}-E_{j})t}\textrm{Tr}_{B}\left[\rho_{B}(E_{i})\rho_{B}(E_{j})\right]},$
(10)

and note that the time-dependent second term only involves overlaps of
distinct energy levels.

###  II.1 Random matrix ensembles

To understand the essential features of the PSFF we first analyze it in RMT,
allowing for an exact determination of the PSFF. We choose Hamiltonians $H$
(Floquet operators $V$) from the canonical Wigner-Dyson random matrix
ensembles [20, 21, 24, 1], yielding time evolution operators $T(t)=\exp(-iHt)$
[$T(t=\tau n)=V^{n}$]. To evaluate the ensemble average in Eq. (10), we can
utilize that for these RMT ensembles the eigenvalues and eigenstates of $H$
($V$) are uncorrelated. Thus, their ensemble average factorizes and can be
performed independently. We find

$\displaystyle
K_{A}(t)=\frac{\mathcal{P}_{B}-\mathcal{Q}_{B}}{D_{A}}+{{D}_{B}}\mathcal{Q}_{B}K(t)~{},$
(11)

where $\mathcal{Q}_{B}=(D(D-1))^{-1}\overline{\sum_{i\neq
j}\textrm{Tr}_{B}\left[\rho_{B}(E_{i})\rho_{B}(E_{j})\right]}$ and
$\mathcal{P}_{B}$ are the averaged overlap and purities of the reduced
eigenstates, respectively. We note that here the PSFF is the full SFF with a
scaling factor ${{D}_{B}}\mathcal{Q}_{B}$ and a constant subsystem dependent
shift $(\mathcal{P}_{B}-\mathcal{Q}_{B})/D_{A}$ such that the entire time
dependence of the PSFF is captured in the SFF. Therefore, the PSFF in these
models preserves the characteristic ramp-plateau structure and the relevant
time scales of the SFF.

As shown in App. B, we can evaluate $\mathcal{P}_{B}$ and $\mathcal{Q}_{B}$
explicitly using Wigner-Dyson RMT for the eigenstates of $H$ ($V)$. They are
functions of only the Hilbert space dimensions of subsystems $A$ and $B$, i.e.
$\mathcal{P}_{B}\equiv\mathcal{P}_{B}(D_{A},{{D}_{B}})$ and
$\mathcal{Q}_{B}\equiv\mathcal{Q}_{B}(D_{A},{{D}_{B}})$. The precise
functional form of $\mathcal{P}_{B}$ and $\mathcal{Q}_{B}$ depends on the
symmetry class of the Hamiltonian $H$ (Floquet operator $V$). For the case of
the unitary Wigner-Dyson ensembles, for example $H$ from the Gaussian unitary
ensemble or $V$ from CUE, we find

$\displaystyle\mathcal{P}_{B}=\frac{D_{A}+{{D}_{B}}}{D_{A}{{D}_{B}}+1}\quad;\quad\mathcal{Q}_{B}=\frac{{{D}_{B}}\left({D_{A}}^{2}-1\right)}{{D_{A}}^{2}{{{D}_{B}}}^{2}-1}~{}.$
(12)

The analogous expressions for orthogonal Wigner-Dyson ensembles can be found
in App. B. In both symmetry classes at $D_{A},D_{B}\gg 1$, we find that,
$\mathcal{P}_{B}-\mathcal{Q}_{B}\approx 1/D_{A}$ and $\mathcal{Q}_{B}\approx
1/{{D}_{B}}$. Thus, in this limit, the PSFF has a constant shift of
$1/D_{A}^{2}$ added to the SFF and the slope of the ramp is the same as the
slope of the ramp in the SFF, i.e. $K_{A}(t)\approx K(t)+1/D_{A}^{2}$ [see
also Eq. (5)].

###  II.2 General chaotic systems

In the case of more general chaotic systems, we begin by separating out the
reduced density matrices of the energy eigenstates into smooth and fluctuating
functions of energy,

$\rho_{B}(E)=\frac{\mathbb{1}}{D_{B}}+\Delta\rho_{B}(E)+\delta\rho_{B}(E)~{}.$
(13)

Here, the first term is a constant corresponding to a maximally mixed reduced
density matrix; $\Delta\rho_{B}(E)$ is traceless and a smooth function of $E$,
while $\delta\rho_{B}(E)$ is again traceless but required to fluctuate rapidly
with $E$. For our present purposes, it is useful to define the smooth and
fluctuating parts in terms of their Fourier transforms with respect to a
continuous energy variable as follows: for some cutoff time $t_{\rho}\ll
O(D)$, we take their respective Fourier transforms to satisfy
$(\Delta\tilde{\rho}_{B}(t))_{jk}=0$ for $\lvert t\rvert>t_{\rho}$, and
$(\delta\tilde{\rho}_{B}(t))_{jk}=0$ for $\lvert t\rvert<t_{\rho}$ (with some
additional details in App. C). The essence of the definition is that as a
function of energy, the smooth part varies only over scales much larger than
some energy window of size $t_{\rho}^{-1}$ containing several levels, while
the fluctuating part varies only over scales much smaller than
$t_{\rho}^{-1}$.

We will further assume that $\delta\rho_{B}(E)$ behaves as if it is
‘randomized’ within these energy windows over the ensemble i.e. it is
uncorrelated with the smooth part and satisfies
$\overline{\textrm{Tr}_{B}\left[\delta\rho_{B}(E_{i})\delta\rho_{B}(E_{j})\right]}=\delta_{ij}\overline{\textrm{Tr}_{B}\left[\delta\rho_{B}^{2}(E_{i})\right]}$
for $E_{i},E_{j}$ closer than $\sim t_{\rho}^{-1}$, fluctuating around an
average of zero (we do not require this behavior to persist over larger energy
scales $\lvert E_{i}-E_{j}\rvert\gtrsim t_{\rho}^{-1})$. We note that this
assumption is consistent with the general picture of random behavior over
small energy windows in chaotic systems [27], and we can justify it more
generally (irrespective of whether the system/ensemble is chaotic) as follows.
In evaluating the SFF $K(t)$, the ensemble is usually chosen to have
sufficiently large disorder so that the energy levels are randomly distributed
over some large energy window, across different ensemble realizations. This is
necessary to eliminate the erratic fluctuations of the SFF at large $t$ that
depend on the precise positions of levels, and obtain a smooth ensemble-
averaged behavior (see e.g. Refs. [59, 68] for further discussion of this
point). Our assumption is essentially that, this random redistribution of
levels over different ensemble realizations extends to an energy window of
$\sim t_{\rho}^{-1}$, effectively randomizing the fluctuations
$\delta\rho_{B}(E)$ faster than this scale, while $\Delta\rho_{B}(E)$ which
varies over scales larger than this energy window is not randomized in this
manner. We also note that the eigenstates of a given ensemble realization
themselves may additionally be random superpositions of those of a different
realization, e.g. generally randomly mixing all eigenstates of the latter
within the energy window in fully chaotic systems (i.e. systems with no
‘physical’ conserved quantities other than energy) [2, 94, 95, 96], which
gives further weight to this assumption.

####  II.2.1 Shift-ramp-plateau structure of the PSFF

Using the form in Eq. (13), the overlaps occurring in the definition of the
PSFF in Eq. (I.2) separate out into independent contributions from each part
of the reduced density matrix - the cross terms vanish, due to tracelessness
for terms involving overlaps with the maximally mixed part, or due to the
randomization of $\delta\rho_{B}(E)$ for terms involving overlaps of the
smooth and fluctuating part for $t\gg t_{\rho}$. We can write this as,

$K_{A}(t\gg t_{\rho})=K(t)+\Delta K_{A}(t)+\delta K_{A}(t)~{},$ (14)

where $\Delta K_{A}(t)$ involves only overlaps of the form
$\textrm{Tr}_{B}\left[\Delta\rho_{B}(E_{i})\Delta\rho_{B}(E_{j})\right]$ and
similarly, $\delta K_{A}(t)$ involves only those of the form
$\textrm{Tr}_{B}\left[\delta\rho_{B}(E_{i})\delta\rho_{B}(E_{j})\right]$. On
decomposing $\delta K_{A}(t)$ in a manner analogous to Eq. (10), it follows
that its time dependent part for $t\gg t_{\rho}$ (which sees contributions
only from variations of the overlaps of fluctuating parts within energy
windows smaller than $\sim t_{\rho}^{-1}$) vanishes on ensemble averaging, an
important consequence of the randomization of $\delta\rho_{B}(E)$. This leaves
only a constant contribution from the purity of the fluctuating part, $\delta
K_{A}(t\gg t_{\rho})=\delta\mathcal{P}_{B}/D_{A}$, where
$\delta\mathcal{P}_{B}\equiv
D^{-1}\overline{\sum_{i}\textrm{Tr}_{B}\left[\delta\rho_{B}^{2}(E_{i})\right]}$
(here we use ‘purity’ to generally mean $\textrm{Tr}\left[x^{2}\right]$ for a
Hermitian operator $x$). We see that this constant late-time shift is a
generic feature of the PSFF, independent of the specific form of the full SFF
$K(t)$. It merges into the plateau of the PSFF when $K(t)$ and $\Delta
K_{A}(t)$ show only a plateau behavior - and therefore, the shift is an
independent observable only if the other two terms show non-trivial time
dependence at late times $t\gg t_{\rho}$.

We note that $\Delta K_{A}(t)$ is modulated only by a smooth function of two
energy variables varying over scales larger than $t_{\rho}^{-1}$. For $t\gg
t_{\rho}$, it should then essentially see the contribution to $K(t)$ from each
part of the spectrum but modulated by the value of the function for nearly
equal energies in that part. In App. C, we show this by direct calculation for
a fully chaotic system with Wigner-Dyson level statistics, obtaining a
modulated linear ramp and plateau in addition to the late-time shift, for
$t\gg t_{\rm Th},t_{\rho}$,

$\displaystyle K_{A}(t\gg t_{\rm
Th},t_{\rho})=\frac{\delta\mathcal{P}_{B}}{D_{A}}\qquad\qquad\qquad\qquad\qquad$
$\displaystyle+$ $\displaystyle\frac{1}{D}\begin{cases}(\beta\pi D)^{-1}\gamma
t\left(1+D_{B}\widetilde{\Delta\mathcal{P}}_{B}\right)&\text{for }t\ll
t_{H}~{},\\\ 1+D_{B}\Delta\mathcal{P}_{B}&\text{for }t\gg
t_{H}~{}.\end{cases}$ (15)

Here, $\beta=1,2$ respectively for the orthogonal and unitary classes, while
$\gamma=\sum_{i}\Omega^{-1}(E_{i})$ is the range of energies in the spectrum
with $\Omega(E)$ representing the (smoothened) local density of states, in
agreement with known results for the full SFF (see e.g. Refs. [64, 97]). To
keep the expressions simple, we are ignoring corrections that are prominent
near $t\sim t_{H}$ [see, for instance, the exact form of the GOE SFF in Eq.
(40)]; we focus instead on the $t\ll t_{H}$ regime where the ramp appears
linear for all values of $\beta$ and profiles of $\Omega(E)$, and the $t\gg
t_{H}$ regime with a constant plateau. However, both expressions are exact
throughout the range of times when $\beta=2$ with constant density of states
$\Omega(E)=t_{H}/(2\pi)$. We have also defined two ensemble-averaged
quantities corresponding to slightly different spectrum averages of the purity
of the smooth part,
$\Delta\mathcal{{P}}_{B}=D^{-1}\overline{\sum_{i}\textrm{Tr}_{B}\left[\Delta\rho_{B}^{2}(E_{i})\right]}$
and
$\widetilde{\Delta\mathcal{{P}}_{B}}=\gamma^{-1}\overline{\sum_{i}\Omega^{-1}(E_{i})\textrm{Tr}_{B}\left[\Delta\rho_{B}^{2}(E_{i})\right]}$,
the latter including the contribution to the coefficient of the linear ramp
from each part of the spectrum. We note that the purities of the smooth and
fluctuating parts are (exactly) related to the overall average purity by
$\mathcal{P}_{B}=D_{B}^{-1}+\Delta\mathcal{P}_{B}+\delta\mathcal{P}_{B}$,
giving the expected plateau value of $\mathcal{P}_{B}/D_{A}$ in Eq. (15).
There are also two competing time scales for the onset of the ramp, $t_{\rm
Th}$ and $t_{\rho}$ \- the former entirely determines the behavior of $K(t)$
but the latter appears in $\Delta K_{A}(t)$ and $\delta K_{A}(t)$.

For direct comparison with numerics, it is useful to define the ensemble
averaged overlap of adjacent states,
$Q_{B}=(D-1)^{-1}\sum_{i}\overline{\textrm{Tr}_{B}\left[\rho_{B}(E_{i})\rho_{B}(E_{i+1})\right]}$.
Using Eq. (13), we note that,

$\displaystyle\mathcal{Q}_{B}$
$\displaystyle=\frac{1}{{{D}_{B}}}+\Delta\mathcal{P}_{B}~{},$
$\displaystyle{\mathcal{P}_{B}}-{\mathcal{Q}_{B}}$
$\displaystyle=\delta\mathcal{P}_{B}~{},$ (16)

which follow from the assumption of uncorrelated $\delta\rho_{B}(E)$ in the
ensemble, and taking $\Delta\rho_{B}(E_{i})\simeq\Delta\rho_{B}(E_{i+1})$. We
note that this definition of $Q_{B}$ is equivalent to that in Sec. II.1 for
random matrix ensembles, where the ensemble averaged overlaps between distinct
states are independent of their energies. Sec. III will directly use
$\mathcal{P}_{B}$ and $\mathcal{Q}_{B}$, with the implicit assumption that
$\widetilde{\Delta P}_{B}$ is of similar order of magnitude to
$\Delta\mathcal{P}_{B}$ (due to $\Omega(E)$ being of a similar order of
magnitude throughout the spectrum) and is therefore similarly well represented
by $\mathcal{Q}_{B}$.

####  II.2.2 Constraints from eigenstate thermalization

We have seen that at late times, the PSFF preserves the characteristic
features of the SFF, such as the ramp and the Heisenberg time (as in Eq. (15)
for fully chaotic systems). However, there are non-negative subsystem-
dependent parameters $\mathcal{P}_{B}$, $\delta\mathcal{P}_{B}$ and
$\Delta\mathcal{P}_{B}$ ($\sim\widetilde{\Delta\mathcal{P}}_{B}$) that
respectively influence the plateau value, the magnitude of the shift and the
magnitude i.e., slope of the ramp. The purity $\mathcal{P}_{B}$ measures the
extent of delocalization of eigenstates in a physical basis (e.g. a product
basis of qubits), while we will see that $\delta\mathcal{P}_{B}$ and
$\Delta\mathcal{P}_{B}$ are complementary probes of thermalization of these
eigenstates. Specifically, we mean thermalization in the sense of ETH - that
eigenstates corresponding to sufficiently close energies show nearly identical
behavior in the dynamics of few-body observables [2, 3, 25, 26, 27, 28].

For our purposes, it is convenient to use subsystem ETH [29], which amounts to
imposing ETH on an entire subsystem i.e. for all observables in the subsystem,
and is directly expressed in terms of reduced density matrices. It can be
interpreted as the requirement of a small fluctuating part for the reduced
density matrices of thermal eigenstates, as opposed to large fluctuations for
non-thermal eigenstates. We can therefore apply it directly to the
decomposition of reduced density matrices in Eq. (13). An important advantage
of this version of ETH is that the dependence on subsystem size is made more
explicit, whereas more conventional statements of ETH restrict themselves to
few body operators, corresponding to extremely small subsystems and therefore
negligible subsystem dependence. This subsystem size dependence will turn out
to be the primary non-trivial indicator of the properties of eigenstates in
the PSFF.

In App. D, we discuss the general constraints from (an extension of) subsystem
ETH for eigenstates with an arbitrary extent of delocalization in a physical
basis. Here, we present the results for a system with fully delocalized
eigenstates, characterized by subsystem purities that follow the volume law of
entanglement [31],

$\mathcal{P}_{B}=D_{B}^{-1}+O(D_{B}^{-1})+O(D_{A}^{-1}),$ (17)

which cannot be less than $D_{B}^{-1}$ as well as $D_{A}^{-1}$. This is the
case relevant for the numerical examples of Sec. III. If these eigenstates are
thermal, subsystem ETH requires the smooth and fluctuating parts to satisfy,

$\Delta\mathcal{P}_{B}=O(D_{B}^{-1}),\ \delta\mathcal{P}_{B}=O(D_{A}^{-1}).$
(18)

Non-thermal eigenstates are characterized by much larger fluctuations,
$\delta\mathcal{P}_{B}\gg O(D_{A}^{-1})$, with $\Delta\mathcal{P}_{B}$ being
correspondingly smaller so as to satisfy the constraint
$\mathcal{P}_{B}=D_{B}^{-1}+\Delta\mathcal{P}_{B}+\delta\mathcal{P}_{B}$. A
narrower class of such chaotic systems (e.g. Floquet systems) have uniformly
random eigenstates that are distributed in close agreement with the standard
random matrix ensembles (Sec. II.1); the leading forms of the corresponding
exact results in Eq. (12) are seen to be consistent with Eqs. (17),(18), on
relating the two using Eq. (16). In this context, we note that Ref. [80] has
observed subleading corrections to the random matrix prediction for
eigenstates in 1D Floquet quantum circuits.

###  II.3 Localized systems

Now, we consider localized systems, which show Poisson level statistics (i.e.
uncorrelated neighboring levels) with localized non-thermal eigenstates, for
strong disorder [30, 31]. Here, $K(t)$ shows only a plateau at late times,
allowing us to access only the purity $\mathcal{P}_{B}$ through the PSFF.
Fully localized states are essentially nearly pure states with
$\mathcal{P}_{B}\sim O(1)$ (more precisely, following an area law of
entanglement [31]), and additionally have large fluctuations
$\delta\mathcal{P}_{B}\sim O(1)\leq 1-D_{B}^{-1}$. In other words, fully
localized states cannot thermalize, as they would have to be distributed over
different physical basis states due to orthogonality. An $O(1)$ plateau value
is therefore all we need to characterize the eigenstates of such systems.

On the other hand, when the eigenstates become more delocalized in the
approach to a chaotic phase, thermalization becomes a possibility. The moment
any non-trivial correlations between nearby energy eigenvalues emerge in the
spectrum, leading to a time dependence of $K(t)$ for $t>t_{\rho}$,
$\delta\mathcal{P}_{B}$ becomes a meaningful observable in the PSFF according
to the discussion following Eq. (14). Here, the PSFF can be used to study the
extent of thermalization in addition to the delocalization of the eigenstates.

###  II.4 Summary

Let us summarize the main conclusions of this section from a unified
perspective, before moving on to illustrate them with numerical examples in
the next section. The PSFF in a subsystem $A$ combines energy level
statistics, as reflected in the SFF, with the purities and overlaps of the
reduced energy eigenstates in the complementary subsystem B. The plateau value
of the PSFF encodes the (spectrum and ensemble averaged) purity, which is
$\sim O(1)$ in a fully localized phase, and small for fully delocalized states
in accordance with the volume law of entanglement, Eq. (17). Something more
interesting happens at late times if the SFF has a ramp or other time-
dependent feature due to the existence of local level correlations. The PSFF
inherits the ramp, but the ramp couples only to the smooth, slowly varying
part of the reduced energy eigenstates. The rapidly fluctuating part is left
over as a nearly time-independent shift [Eq. (15)].

Eigenstate thermalization is primarily encoded in the size of the fluctuating
part as measured by the shift - namely, an exponential suppression of the
latter with subsystem size $N_{A}$ is indicative of thermalization [Eq. (18)],
while the lack of such a suppression translates to a failure of the
eigenstates to thermalize. The smooth part is correspondingly large for
thermal eigenstates and small for non-thermal eigenstates, so as to preserve
the overall purity (i.e. extent of delocalization). Finally, there are special
systems for which much more precise predictions for the PSFF can be
theoretically derived/motivated and tested, such as chaotic Floquet systems
with their random matrix-like eigenstates [Eqs. (11) and (12)].

Thus, the PSFF complements the SFF in analyzing late-time quantum chaos by
being able to probe if the eigenstates satisfy ETH, in addition to (and
because of) capturing information about level correlations as contained in the
ramp of the SFF. In particular, we expect that it could potentially be useful
in studying the joint emergence or loss of Wigner-Dyson level statistics and
eigenstate thermalization (which are formally independent notions of late time
quantum chaos) and their interdependence, across a transition or crossover
between a chaotic and non-chaotic phase. This could be done by tuning the
parameters of a system (say, in a quantum simulator) between such phases, and
measuring PSFFs across different choices of subsystems of different sizes -
analyzing the extent of delocalization of eigenstates in the absence of a ramp
via the plateau value, and additionally the extent of thermalization through
the value of the shift if a ramp or other time-dependent feature is present at
late times. Among the interesting possibilities that have been considered for
such an intermediate regime, which could conceivably be probed with the PSFF,
is the existence of so-called non-ergodic extended states [98, 99, 100, 101,
102, 103] where the eigenstates are incompletely delocalized but do not
thermalize, or alternatives in which the eigenstates thermalize without being
fully delocalized [104].

##  III Partial Spectral Form Factor: Numerical Results

Having discussed features of the PSFF and its connection to the SFF utilizing
Wigner-Dyson random matrix ensembles and the ETH, we now present our numerical
results of PSFFs in locally interacting many-body models, as realized in
quantum simulators. For this purpose, we focus on two examples: the Floquet
model Eq. (2) and the Hamiltonian model Eq. (19). Our results are in agreement
with the analysis of the previous Sec. II, in particular regarding the orders
predicted for the averaged purity $\mathcal{P}_{B}$ and the overlap
$\mathcal{Q}_{B}$ via Eq. (16). We consider the Floquet model in the chaotic
phase and the Hamiltonian model in both the chaotic and MBL phases.

####  III.0.1 Example 1: Floquet system

Figure 3: Results for the Floquet $V_{3}$ model. (a) The SFF and PSFF are
presented for $N=6$, $N_{A}=3,~{}4,~{}5$ in red colors. In gray, we plot the
same quantities in a CUE model. (b) The plateau value $K(\infty)$ multiplied
with the subsystem dimension $D_{A}$ is plotted in black circles and matches
with the averaged purity $\mathcal{P}_{B}$ plotted with red crosses. The
average overlap $\mathcal{Q}_{B}$ and the difference
$\mathcal{P}_{B}-\mathcal{Q}_{B}$ are presented in brown and green
respectively. We observe a perfect match with the respective quantities in CUE
plotted in gray, indicating the same averaged eigenvalue and eigenstate
statistics in CUE and $V_{3}$. In the numerical computation, we have taken
8000 disorder realizations to perform ensemble averaging and the subsystems
$A$ are chosen from the middle of the spin chain.

The Floquet time evolution operator $V_{3}$ has the same quasi-energy
eigenvalue statistics as the CUE random matrix ensemble [81, 58]. As mentioned
in Sec. I the Floquet models are known to thermalize to infinite temperatures
as per RMT and thus we expect the eigenstate statistics to also be the same as
in the corresponding RMT class. To show this, we present in Fig. 3(a)
numerically obtained SFF and PSFF for a total system size of $N=6$ and
subsystem sizes $N_{A}=3,4$ and $5$ for the model $V_{3}$. We plot with gray
lines the corresponding $K_{A}(t)$ in a CUE model where the analytic forms can
be exactly calculated (see Sec. II.1 and App. B). For the PSFF $K_{A}(t)$ at
$N_{A}=3$ and very early times, we notice that the onset of the ramp takes a
few initial periods to set, but eventually the PSFF follows the CUE
prediction.

The closeness between the statistics of CUE and $V_{3}$ can further be seen
from the average overlaps of reduced densities of eigenstates
$\mathcal{P}_{B}$ and $\mathcal{Q}_{B}$. In Fig. 3(b) we present the average
purity and overlaps as functions of subsystem size $N_{A}$. At plateau time,
$t>t_{H}(=D\tau)$ the PSFF becomes
$K_{A}(t\rightarrow\infty)=\mathcal{P}_{B}/D_{A}$, see Eq. (10). We plot
numerically obtained $K_{A}(\infty)D_{A}$ in black circles, and the average
purity $\mathcal{P}_{B}$ with red crosses, they confirm the analytic
expectation. The average overlap $\mathcal{Q}_{B}$ and the difference
$\mathcal{P}_{B}-\mathcal{Q}_{B}$ are plotted in brown and green circles
respectively and match with the CUE data.

To conclude, the SFF, PSFF, averaged purity and overlaps match in the CUE and
$V_{3}$ model and thus we expect the form of the PSFF in Eq. (5) to hold for
the model $V_{3}$, after a small initial time period. We know from Eq. (12),
for large Hilbert space dimensions, that $\mathcal{Q}_{B}\approx 1/{{D}_{B}}$
and $\mathcal{P}_{B}-\mathcal{Q}_{B}\approx 1/D_{A}$. Therefore utilizing, Eq.
(16), we find that $\Delta\mathcal{P}_{B}=0$ and
$\delta\mathcal{P}_{B}=O(1/D_{A})$ for $V_{3}$ and the RMT models. The purity
of the smooth part (of the form of $\mathrm{Tr}[\Delta\rho_{B}^{2}(E)]$)
appears in the ramp part of the PSFF in Eq. (15) and thus we note that the
ramp coefficient is $\sim 1/D^{2}$ for $D_{A}\gg 1$. On the other hand, the
purity of the fluctuating part (of the form of
$\mathrm{Tr}[\delta\rho_{B}^{2}(E)]$) comes in the time-independent term added
to the SFF in Eq. (15), which is to the leading orders $~{}1/D_{A}^{2}$, as
also in the CUE model [Eq. (5)]. To further have another numerical example of
the Floquet model thermalizing according to RMT, we present the example of a
chaotic Floquet model with time-reversal symmetry in App. E.

####  III.0.2 Example 2: Hamiltonian system

As our second example, we consider a transverse field Ising model in presence
of longitudinal local disorders,

$H=J\left(\sum_{\begin{subarray}{c}i,j=1\\\
i<j\end{subarray}}^{N}\frac{1}{(i-j)^{\alpha}}\sigma^{z}_{i}\sigma^{z}_{j}+\sum_{i=1}^{N}\sigma^{x}_{i}\right)+W\sum_{i=1}^{N}h_{i}\sigma^{z}_{i},$
(19)

where $h_{i}$ are drawn uniformly at random from $(-1,1)$. The coefficient $J$
and the exponent $\alpha$ denote the strength and range of the interactions
respectively. The disorder strength $W$ is known to specify the nature of the
dynamics; $W\sim J$ depicts chaotic regime and $W\gg J$ corresponds to the
localized regime (for a similar model see, [64]). In the App. F.1, we present
the adjacent level gap ratio as a function of $W/J$ and $\alpha$ and find that
the chaotic and localized phases exist for short ($\alpha>1$) as well as for
long ($\alpha<1$) range interactions. In this work, we choose $\alpha=1.2$,
and as examples of the chaotic and localized phases, we take $W=J$ and $W=10J$
respectively. In contrast to the presence of the ramp and plateau in the SFF
for chaotic models, the SFF for localized models stays flat for all times
$t\gg 0$. In the numerics, we will find that the PSFF preserves this flat
feature of the SFF, and has a subsystem dependent shift added over the SFF, as
predicted in Sec. II.3. In Fig. 4 and 5 we present numerical results for the
Hamiltonian model (19) in these two phases. For clarity, we have used red
color for the chaotic phase ($W=J$) and blue for the MBL phase ($W=10J$). We
note that the Hamiltonian of Eq. (19) has the time-reversal symmetry of
complex conjugation in the computational ($\sigma^{z}_{i}$) basis [1, 105,
106]. A chaotic Hamiltonian with this symmetry is known to follow the
eigenvalue statistics (or the SFF) of GOE after the Thouless time $t>t_{\rm
Th}$ [24, 27, 1, 105, 106], thus we have also put the results for GOE class in
gray in Fig. 4.

Figure 4: Results for the Hamiltonian model. In a log-log plot we present the
chaotic phase ($W=J$) in red, MBL phase $(W=10J$) in blue, and the GOE in
gray. In both phases the SFF and PSFF are plotted for (sub-)system sizes
$N_{A}=6,~{}7$ and $N_{A}=N=10$. The SFF for the chaotic phase has the
characteristic ramp and plateau and follows the GOE SFF at late times. The
PSFF in this phase also has the shift, ramp and plateau, we plot these in a
focused linear scale plot in Fig. 5(a). The MBL phase shows a flat SFF and
PSFF for all times $t\gg 0$. The mean level spacing (i.e. the Heisenberg time)
in the MBL phase and GOE are numerically rescaled to match to the one in the
chaotic phase. Figure 5: Results for the Hamiltonian model. (a) In linear
scale we present the SFF and PSFF for the chaotic phase ($W=J$). (Sub-)system
sizes $N_{A}=6$, $7$ and $N_{A}=N=10$ are plotted with magenta, green and red
respectively for both the Hamiltonian model (with solid curves) and the GOE
(with dashes). We observe differences in the PSFF for chaotic Hamiltonian and
GOE. These differences are investigated in (b), (c) and (d) through
$\mathcal{P}_{B}$ and $\mathcal{Q}_{B}$. We use red color for the chaotic
phase ($W=J$) and gray for the GOE. For comparison we have also plotted these
quantities in the localized phase ($W=10J$) using blue color. (b) We plot
$K_{A}(\infty)D_{A}$ using black circles which matches with the corresponding
average purity $\mathcal{P}_{B}$ of the MBL and chaotic phase. (c) The average
overlap $\mathcal{Q}_{B}$ for MBL, chaotic and GOE follow closely the behavior
$1/{{D}_{B}}$. (d) The difference
$\mathcal{P}_{B}-\mathcal{Q}_{B}\approx\delta{\mathcal{P}_{B}}$, which encodes
the shift of the PSFF, is larger for large disorders (MBL) compared to small
disorders (chaotic). In the numerical computation, we have taken 200
Hamiltonians to perform ensemble averaging and the subsystems $A$ are chosen
from the middle of the spin chain.

As a side remark, we emphasize at this point that the spectrum of the local
Hamiltonian model, Eq. (19), does not have the same density of states as the
GOE spectrum and thus the Hamiltonian SFF should be compared with an average
of GOE SFFs, each with $t_{H}$ determined by different parts of the
Hamiltonian spectrum. Often, this is circumvented by removing the non-
universal effects arising from the edges of the local Hamiltonian spectrum by
using a filter function such that only the middle part of the spectrum
contributes [69] or considering very large system sizes where the edge effects
are effectively smaller. In our work, we focus on the measurement of chaotic
features through the observation of the ramp, plateau and the shift which can
already be observed without filtering for moderate system sizes, which we
focus on.

In Fig. 4, the SFF and PSFF are presented for the system size $N=10$ and
subsystem sizes $N_{A}=6$ and $7$. In order to have the same Heisenberg time
$t_{H}$, the eigenvalues are numerically rescaled such that the average mean
level spacing for $W=10J$ match with the one for $W=J$. As a guide, we have
plotted in gray the GOE SFF where the $t_{H}$ is determined from the full
width of the chaotic Hamiltonian spectrum and observe that the SFF for the
chaotic phase follows the GOE SFF closely. The PSFF for the chaotic phase,
shifted up compared to the SFF, also shows the ramp and plateau behavior which
are seen better in a linear plot in Fig. 5(a). Here, focused to display
chaotic features, we have used solid lines for the chaotic Hamiltonian and
dashes for the GOE. The different subsystem sizes are shown in different
colors. We note that the PSFF for the chaotic local model and GOE are
different (see the magenta and green curves). These differences arise due to
the differences in eigenstate properties of the local Hamiltonian and GOE.

Further, to concretely discuss second-moments of eigenstates, in Fig. 5(b) we
present the averaged purity $\mathcal{P}_{B}$ using crossed markers. We have
also plotted here the plateau values $K_{A}(\infty)D_{A}$ (in black circles)
for both chaotic and MBL phases which agree with their respective purities
following $K_{A}(t\rightarrow\infty)=\mathcal{P}_{B}/D_{A}$ [see Eq. (10)].
Note that these average purities are consistent with a volume law of
entanglement in the chaotic phase, and an area law in the localized phase
[31]. For the remainder of this section, it is useful to discuss the two
phases $W=J$ and $W=10J$ separately.

For the chaotic phase $W=J$, the average overlaps $\mathcal{Q}_{B}$ and
$\mathcal{P}_{B}-\mathcal{Q}_{B}$ are presented in red in the bottom panel of
Fig. 5 as functions of $N_{A}$. Assuming ETH for the chaotic systems, we have
discussed orders of magnitude of these overlaps in Sec. II.2. Utilizing Eq.
(16) we can comment on the orders of $\Delta\mathcal{P}_{B}$ and
$\delta\mathcal{P}_{B}$ (see App. F.2 for more details on the numerical
extraction of these orders). From $\mathcal{Q}_{B}$ [Fig. 5(c)], we find
$\Delta\mathcal{P}_{B}=O(1/{{D}_{B}})$ and from
$\mathcal{P}_{B}-\mathcal{Q}_{B}=\delta\mathcal{P}_{B}$ [Fig. 5(d)], we find
$\delta\mathcal{P}_{B}\sim O(1/D_{A})$, confirming the ETH predictions for
chaotic systems. We verify that the value of the shift of the PSFF in the
linear ramp region is given in terms of the purity of the fluctuating part
i.e., by $\delta\mathcal{P}_{B}/D_{A}$ in App. F.3. For comparison, we have
plotted the same quantities in a GOE model in gray. We note a difference
between the overlaps (properties of the eigenstates) in the local chaotic
Hamiltonian and GOE, which is not surprising because the statistics of
eigenstates need not be the same in the two models.

Next, we look at the orders of magnitude of the overlaps in the phase $W=10J$,
plotted in blue in the bottom panel of Fig. 5. Following Eq. (16) from the
$\mathcal{Q}_{B}$ [Fig. 5(c)], we find $\Delta\mathcal{P}_{B}=O(1/{{D}_{B}})$
and from $\mathcal{P}_{B}-\mathcal{Q}_{B}=\delta\mathcal{P}_{B}$ [Fig. 5(d)],
we find $\delta\mathcal{P}_{B}\sim O(1)\gg O(1/D_{A})$. The localized phase is
not expected to satisfy ETH, and as discussed in the Sec. II.3, we expect such
large shift in the PSFF in MBL systems. Due to larger $\delta\mathcal{P}_{B}$
in the MBL phase, we notice a larger overall shift of the PSFF in the MBL
phase, shown in blue in Fig. 5(b)-(d).

##  IV Proof of the protocol

In Sec. I.3, we presented our measurement protocol and defined estimators for
the SFF and PSFF [Eqs. (7) and (8)] in terms of the measured bitstrings. In
this section, we prove analytically that these are unbiased estimators of the
SFF and PSFF utilizing the theory of unitary $2$-designs.

###  IV.1 Useful results from unitary $2$-designs

Unitary $n-$designs are ensembles of random unitary matrices, whose averages
of polynomial moments of order up to $n$ coincide with ones of the Haar
measure (or equivalently the CUE) [85]. With the help of Weingarten calculus,
these moments can be expressed analytically [107], allowing us to relate the
statistics of randomized measurements to the quantity that we would like to
measure. Since the measured bitstrings from the protocol are sampled from the
Born probabilities $|\bra{\mathbf{s}}U^{\dagger}T(t)U\ket{\mathbf{0}}|^{2}$
which are polynomial functions of order two in $U$, we restrict ourselves to
Weingarten calculus of order two. Using independent local unitaries
$U=\bigotimes_{i}u_{i}$, one finds for any operator $C$ defined on the ‘two-
copy’ Hilbert space $\mathcal{H}^{\otimes 2}$ [108]

$\mathbb{E}_{U}\left[(U\otimes U)\,C\,(U^{\dagger}\otimes
U^{\dagger})\right]=\sum_{\sigma,\tau}w_{\sigma,\tau}\textrm{Tr}\left[\sigma
C\right]\tau.$ (20)

Here, $\mathbb{E}_{U}$ denotes the average over local unitaries of the form
$U=\bigotimes_{i}u_{i}$ with $u_{i}$ sampled for each $i$ independently from a
unitary $2$-design on the local Hilbert space $\mathbb{C}^{\otimes 2}$.
Further, the sum extends to all two-copy permutation operators
$\sigma=\bigotimes_{i}\sigma_{i}$ and $\tau=\bigotimes_{i}\tau_{i}$ with
$\sigma_{i},\tau_{i}=\mathbb{1}_{i},\mathbb{S}_{i}$. Here, the identity
$\mathbb{1}_{i}$ and the swap operator $\mathbb{S}_{i}$ act as
$\mathbb{1}_{i}\ket{s_{i}}\otimes\ket{s^{\prime}_{i}}=\ket{s_{i}}\otimes\ket{s^{\prime}_{i}}$
and
$\mathbb{S}_{i}\ket{s_{i}}\otimes\ket{s^{\prime}_{i}}=\ket{s^{\prime}_{i}}\otimes\ket{s_{i}}$
on local basis states $\ket{s_{i}}$ and $\ket{s^{\prime}_{i}}$. Finally, the
coefficient
$w_{\sigma,\tau}=\prod_{i}\textup{Wg}^{U(2)}(\sigma_{i}\tau_{i}^{-1})$ is
determined by the Weingarten function $\textup{Wg}^{U(2)}$, with
$\textup{Wg}^{U(2)}(\mathbb{1}_{i})=1/3$ and
$\textup{Wg}^{U(2)}(\mathbb{S}_{i})=-1/6$. The expression above, which is
valid for any operator $C$, is the mathematical backbone of randomized
measurements. In randomized measurement protocols, the goal is then to
identify an operator $C$, whose expectation value can be inferred from the
experimental data, such that the right hand side of the above equation reveals
the quantity of interest.

In order to reconstruct the SFF, it will turn out to be particularly useful to
choose $C=O\otimes\rho_{0}$ with $\rho_{0}=\ket{\mathbf{0}}\bra{\mathbf{0}}$
and

$\displaystyle O$
$\displaystyle=(\ket{0}\bra{0}-\frac{1}{2}\ket{1}\bra{1})^{\otimes
N}=\sum_{\mathbf{s}}(-2)^{|\mathbf{s}|}\ket{\mathbf{s}}\bra{\mathbf{s}}$ (21)

where the sum extends to all bitstrings $\mathbf{s}=(s_{1},\dots,s_{N})$ with
$s_{i}\in\\{0,1\\}$, and $|\mathbf{s}|\equiv\sum_{i}s_{i}$. For this choice,
we obtain

$\mathbb{E}_{U}\left[(U\otimes U)\,(O\otimes\rho_{0})\,(U^{\dagger}\otimes
U^{\dagger})\right]=4^{-N}\mathbb{S}$ (22)

with
$\mathbb{S}=\bigotimes_{i}\mathbb{S}_{i}=\sum_{\mathbf{s},\mathbf{s^{\prime}}}\ket{\mathbf{s^{\prime}}}\bra{\mathbf{s}}\otimes\ket{\mathbf{s}}\bra{\mathbf{s^{\prime}}}$.
The Swap operation $\mathbb{S}$ is the key operation to extract non-trivial
quantities, such as the purity, in randomized measurements [108]. Here, to
access the SFF, it is convenient to take the partial transpose operation
$A\otimes B\to A^{T}\otimes B$ in the above equation, leading to

$\mathbb{E}_{U}\left[(U^{*}\otimes U)(O^{T}\otimes\rho_{0})(U^{T}\otimes
U^{\dagger})\right]=2^{-N}\ket{\Phi_{N}^{+}}\bra{\Phi_{N}^{+}},$ (23)

where
$\ket{\Phi_{N}^{+}}=\bigotimes_{i}\ket{\Phi_{i}^{+}}=2^{-N/2}\sum_{s}\ket{\mathbf{s}}\otimes\ket{\mathbf{s}}$
is a product of Bell pairs
$\ket{\Phi_{i}^{+}}=2^{-1/2}(\ket{0}\otimes\ket{0}+\ket{1}\otimes\ket{1})$.

###  IV.2 Rewriting the SFF in a form suitable for randomized measurements

For clarity, we focus on the measurement of the full SFF $K(t)$, and present
the case of the PSFF in App. G. We first define for a fixed time-evolution
operator $T(t)$

$\displaystyle K_{T(t)}\equiv
4^{-N}\textrm{Tr}\left[T(t)\right]\textrm{Tr}\left[T^{\dagger}(t)\right]$ (24)

such that the ensemble (disorder) average $K(t)=\overline{K_{T(t)}}$ yields
the SFF, according to the definition Eq. (1). Secondly, we show that
$K_{T(t)}$ equals the survival probability of the Bell State
$\ket{\Phi_{N}^{+}}$ under the dynamics generated by $\mathbb{1}\otimes T(t)$,
i.e.

$\displaystyle K_{T(t)}$
$\displaystyle=\braket{\Phi_{N}^{+}}{\mathbb{1}\otimes{T}(t)}{\Phi_{N}^{+}}\braket{\Phi_{N}^{+}}{\mathbb{1}\otimes{T^{\dagger}}(t)}{\Phi_{N}^{+}}.$
(25)

To this end, we use the following identity for any two operators $A,B$ on
$\mathcal{H}$

$\displaystyle\textrm{Tr}\left[AB\right]=2^{N}\braket{\Phi_{N}^{+}}{A^{T}\otimes
B}{\Phi_{N}^{+}},$ (26)

which can be proven by inserting the definition of the Bell state
$\ket{\Phi_{N}^{+}}=2^{-N/2}\sum_{s}\ket{\mathbf{s}}\otimes\ket{\mathbf{s}}$
in terms of computational basis states. Eq. (25) follows directly by choosing
$A=\mathbb{1}$ and $B=T(t)$. We note that the identity Eq. (25) has been
discussed in the context of holographic duality [109]. In this case
generalized finite temperature form factors can be written in terms of
thermofield double-states, which take the form of Bell states in the limit of
infinite temperature. With the help of Eq. (23), we can now replace one Bell
state projector in Eq. (25) with $O\otimes\rho_{0}$ averaged over random
unitaries $U$. We find $K_{T(t)}=\mathbb{E}_{U}\left[K_{T(t),U}\right]$ with
$K_{T(t),U}$ defined as

$\displaystyle K_{T(t),U}$ $\displaystyle\equiv
2^{N}\bra{\Phi_{N}^{+}}U^{*}O^{T}U^{T}\otimes{T}(t)U\rho_{0}U^{\dagger}{T^{\dagger}}(t)\ket{\Phi_{N}^{+}}.$
(27)

Using once more the identity (26), it follows that $K_{T(t),U}$ equals the
expectation values of the operator $O$ in the final state $\rho_{f}(t)$

$\displaystyle K_{T(t),U}$
$\displaystyle=\textrm{Tr}\left[O\,\right.\underbrace{U^{\dagger}T(t)U\rho_{0}U^{\dagger}T^{\dagger}(t)U}_{\rho_{f}(t)}\left.\right]$
$\displaystyle=\sum_{\mathbf{s}}(-2)^{|\mathbf{s}|}|\bra{\mathbf{s}}U^{\dagger}T(t)U\ket{\mathbf{0}}|^{2}.$
(28)

Here, $|\bra{\mathbf{s}}U^{\dagger}T(t)U\ket{\mathbf{0}}|^{2}$ is precisely
the Born probability of finding a bitstring $\mathbf{s}$, in the computational
basis measurement performed at the end of our measurement sequence when the
state $\rho_{f}(t)$ has been prepared [c.f. Sec. I.3]. It follows thus that

$\displaystyle K_{T(t),U}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{QM}\left[(-2)^{|\mathbf{s}|}\right],$ (29)

where $\mathbb{E}_{QM}$ is the quantum mechanical average and ${\mathbf{s}}$
denotes the outcome of the computational basis measurement at the end of the
measurement sequence.

In summary, it follows that for each measured bitstring ${\mathbf{s}}$,
$(-2)^{|{\mathbf{s}}|}$ provides an estimation of the SFF, which in
expectation over ensemble (disorder) average, over random unitaries and
quantum mechanical averaging, yields the SFF

$\displaystyle
K(t)=\overline{\mathbb{E}_{U}\mathbb{E}_{QM}\left[(-2)^{|{\mathbf{s}}|}\right]}.$
(30)

In practice, we repeat our measurement protocol by performing $M$ independent
experimental runs (with independently sampled time evolution operators and
random unitaries), and calculate the empirical average $\widehat{K(t)}$ [Eq.
(7)]. Using Eq. (30), it follows that $\widehat{K(t)}$ converges to $K(t)$ in
the limit $M\to\infty$. For finite $M$, statistical errors are governed by the
variance of $\widehat{K(t)}$, and are discussed in the next section. In the
App. G, we extend our derivation to the case of the PSFF, and illustrate the
mapping between randomized measurements and the (P)SFF graphically.

##  V Statistical errors and imperfections

We have discussed characteristic features of the SFF and PSFF, such as shift,
ramp and plateau. The crucial question arises whether these can be measured in
today’s quantum simulators, utilizing our protocol (Sec. I.3) with a finite
measurement budget (number of experimental runs $M$) and in the presence of
unavoidable experimental imperfections. In the following, we first analyze in
detail statistical errors which arise from a finite number of experimental
runs $M$. These determine the signal-to-noise ratio for a measurement of the
shift of the PSFF (extracted from measurements at a single point in time) and
the slope of the SFF and PSFF (extracted from differences of measurements at
various points in time). Subsequently, we discuss the influence of
experimental imperfections, such as imperfect implementation of our
measurement protocol or decoherence during the time evolution.

###  V.1 Statistical errors

We discuss statistical errors arising from a finite number of experimental
runs $M$. We first consider the estimation of the SFF and PSFF at single point
in time, and secondly the estimation of (the slope of) the ramp from
measurements of the SFF and PSFF at different times.

####  V.1.1 Observing PSFF and SFF

We can bound the statistical errors of the estimator $\widehat{K_{A}(t)}$ [Eq.
(8)] by its variance. As shown in App. H, we find that,

$\displaystyle\text{Var}[\widehat{K_{A}}]$
$\displaystyle=\frac{1}{M}\left(2^{-N_{A}}\sum_{B\subseteq
A}K_{B}-K_{A}^{2}\right)\equiv\frac{\sigma_{A}^{2}}{M}~{},$ (31)

where we have dropped the time argument for brevity. Here, $K_{B}$ denotes the
PSFF defined in the subsystem $B$ and the sum extends over all subsystems
$B\subseteq A$. The variance of $\widehat{K(t)}$ [Eq. (7)] follows by taking
$A$ to be the full system. We obtain an expected relative error
$\mathcal{E}_{A}=\sigma_{A}/(K_{A}\sqrt{M})$ of an estimation
$\widehat{K_{A}(t)}$ with $M$ experimental runs. As it can be rigorously shown
via Chebyshev’s inequality, the required number of measurements to obtain with
high probability an estimate of $K_{A}(t)$ with fixed relative error scales as
$M\sim\sigma_{A}^{2}/K_{A}^{2}$.

Figure 6: Statistical errors in the Floquet model $V_{3}$. (a) The relative
error, $\mathcal{E}_{A}=\sigma_{A}/(K_{A}\sqrt{M})$ is plotted for total
system size $N=6$ and subsystem sizes $N_{A}=3,4,$ and $6$ as a function of
number of measurements $M$ at time $t/\tau=5$. For a fixed $M$ we perform 100
numerical experiments each with $M$ single shots and present the average
$\mathcal{E}_{A}$ using colored lines. The gray lines represent corresponding
errors in a model with CUE dynamics (calculated analytically in App. H). (b)
The relative error $\mathcal{E}_{A}$ as a function of subsystem size $N_{A}$
at two times $t/\tau=5$ and $t/\tau=30$ is shown. (c) Single PSFF with
subsystem size $N_{A}=4$ for two choices $A=[2,3,4,5]$ (red), $A=[1,2,5,6]$
(green dashed), and the average PSFF $K_{|A|}$ (black) follow each other; the
numbers in the $[\cdots]$ denote qubit index. (d) For the observation of the
ramp we plot the SNR of the slope, SNR$[c_{A}(t_{2},t_{1})]$ (in red) and
SNR[$c_{|A|}(t_{2},t_{1})$] (in black). Both SNRs are constructed from a
single data set of $M=10^{6}$. As a guide to the eye, we also present in gray
the SNR[$c_{|A|}(t_{2},t_{1})$] when all the measurements in the averaged PSFF
are done independently, i.e. when
$\sigma_{|A|}=\sigma_{A}\binom{N}{N_{A}}^{-1/2}$. This would require
$M\binom{N}{N_{A}}$ number of independent measurements.

The expected statistical error $\mathcal{E}_{A}$, and hence also the number of
required experimental runs, depends thus on the value of $K_{A}$ itself, as
well as on the PSFF $K_{B}$ of all subsystems $B\subseteq A$. For Hamiltonians
(Floquet-) operators from Wigner-Dyson RMT, we can explicitly evaluate
$\sigma_{A}$ (see App. H). As the worst-case estimate, we find that at the
point of weakest signal, after a single time step $t=\tau$ in Floquet dynamics
$T(t=n\tau)=V^{n}$ with $V$ sampled from CUE where
$K_{A}(t=1\tau)=2^{-2N_{A}}$, the expected relative statistical error is given
by $\mathcal{E}_{A}=\sqrt{10^{N_{A}}/M}$. A total number of measurements
$M\sim 10^{N_{A}}/\epsilon^{2}\approx 2^{3.32N_{A}}/\epsilon^{2}$ is thus
required to obtain a fixed relative error $\epsilon$. This is to be contrasted
with the number of measurements required for quantum process tomography, which
requires, without strong assumptions on the process of interest [88], at least
$r2^{5N_{A}}/\epsilon^{2}$ measurements, with $r=r(N_{A})\geq 1$ being the
Kraus rank of the process [110]. In addition, we can reduce the exponents
associated with the scaling of statistical errors in randomized measurement
protocols further using importance sampling [48, 111, 112, 113].

In Fig. 6(a) we plot the relative error $\mathcal{E}_{A}$ as a function of the
number of experimental runs $M$ in the $V_{3}$ model (2) at time $t/\tau=5$
with total qubits $N=6$. The relative error decays as $\sim 1/\sqrt{M}$ with
increasing $M$, as expected from the central limit theorem. Furthermore, it
decreases with decreasing subsystem size. This is also shown in Fig. 6(b)
where we display, for a fixed $M$, the relative errors as a function of
subsystem size $N_{A}$ at two different times $t/\tau=5$ and $t/\tau=30$. As
expected, we observe that the relative error is largest at early times where
the PSFF is smallest. At early times, the relative error increases with the
subsystem size, thereby requiring more measurements as $N_{A}\rightarrow N$.

####  V.1.2 Observing the ramp in chaotic models

The relative error $\mathcal{E}_{A}=\sigma_{A}/(K_{A}\sqrt{M})$ determines the
required number of measurements to estimate the PSFF at a single point in
time. While this reveals important information on the overall magnitude and in
particular the ‘shift’ of the PSFF, signatures of energy level repulsion are
encoded in the ramp of the SFF and PSFF (see Sec. II). To detect the ramp, we
aim thus to measure the difference $K_{A}(t_{2})-K_{A}(t_{1})$ at two points
in time $t_{2}>t_{1}$, in particular, the slope of $K_{A}$,

$\displaystyle
c_{A}(t_{2},t_{1})=\frac{K_{A}(t_{2})-K_{A}(t_{1})}{t_{2}-t_{1}}.$ (32)

To quantify the experimental effort to resolve $c_{A}(t_{2},t_{1})$, we
introduce its signal-to-noise ratio $\text{SNR}[c_{A}(t_{2},t_{1})]$, which,
for independent measurements of the PSFF at times $t_{2}$ and $t_{1}$, is
given by

$\displaystyle\text{SNR}[c_{A}(t_{2},t_{1})]=\sqrt{M}\frac{K_{A}(t_{2})-K_{A}(t_{1})}{{\sigma_{A}(t_{2})+\sigma_{A}(t_{1})}}.$
(33)

As shown in Secs. II and III, the slope $c_{A}(t_{2},t_{1})$ of the PSFF (i.e.
the signal), is approximately constant as a function of the subsystem size
$N_{A}\gtrsim N/2$. At the same time, the absolute value of the noise, here
$(\sigma_{A}(t_{2})+\sigma_{A}(t_{1}))/\sqrt{M}$, decreases with increasing
$N_{A}$ (as the absolute value of the PSFF decreases). Thus, as shown in Fig.
6(d) (red curve) for the $V_{3}$ model, $\text{SNR}[c_{A}(t_{2},t_{1})]$
typically increases with increasing subsystem size $N_{A}$, reaching a maximum
when the subsystem is the system itself i.e., $N_{A}=N\,(=6$ in the example
here).

In chaotic quantum systems, our protocol enables detection of the ramp with
further improved SNR: First, we note that the order of magnitude of different
features of the PSFF does not depend on the actual choice of the subsystem
$A$, but only on its size $|A|=N_{A}$. Hence, as numerically shown in Fig.
6(c), we can replace the PSFF $K_{A}$ of a specific subsystem $A$ with its
average

$\displaystyle K_{|A|}(t)=\binom{N}{N_{A}}^{-1}\sum_{|A|=N_{A}}K_{A}(t)~{},$
(34)

where we sum over all subsystems $A$ of fixed size $N_{A}$ (including
disconnected subsystems).

Second, we note that from a single experimental data set, taken on the full
system $\mathcal{S}$, we can estimate $K_{A}(t)$ for all subsystems
$A\subseteq\mathcal{S}$, via spatial restriction in the post-processing. Thus,
we can also obtain the average PSFF $K_{|A|}(t)$ and its slope
$c_{|A|}(t_{2},t_{1})$. Since for $N_{A}<N$, there are multiple subsystems $A$
of size $N_{A}$, we can expect an increased $\text{SNR}[c_{|A|}(t_{2},t_{1})]$
for these average quantities.

In Fig. 6(d), we display the numerically determined signal-to-noise-ratio
$\text{SNR}[c_{|A|}(t_{2},t_{1})]$, for the averaged PSFF in black. Indeed,
compared to the SNR for a single subsystem $A$, SNR[$c_{A}(t_{2},t_{1})$] in
red, we observe an enhanced $\text{SNR}[c_{|A|}(t_{2},t_{1})]$ for subsystem
sizes $1<N_{A}<N$. We remark that we do not reach an enhancement
$\binom{N}{N_{A}}^{1/2}$ of the SNR which would result trivially from
$\binom{N}{N_{A}}$ separate experiments (i.e. $\binom{N}{N_{A}}\cdot M$
experimental runs in total, gray line) since the estimations
$\widehat{K_{A}(t)}$ for various subsystems $A$ from a single data set are not
independent. Nevertheless, Fig. 6(d) shows that the average PSFF $K_{|A|}$,
extracted at a subsystem size $N_{A}\approx N/2$ has the largest SNR for
determining the slope of the ramp from a given measurement dataset. Thus, as
compared to the PSFFs $K_{A}(t)$ for fixed subsystems $A$ or the full SFF
$K(t)$, the average PSFF $K_{|A|(t)}$ at half system size provides a favorable
tool to observe the ramp of the (P)SFF, i.e. signatures of level repulsion in
chaotic quantum many systems.

###  V.2 Experimental imperfections

Figure 7: Experimental imperfections and decoherence. We study effects of
measurement errors and decoherence on the estimated SFF and PSFF $K_{A}(t)$
using the example of the kicked spin $V_{3}$ with total system size $N=4$. In
(a), we display the relative error
$\widehat{\epsilon_{\eta}}=(\widehat{K_{A}}-\widehat{(K_{A})}_{\eta})/\widehat{K_{A}}$
of the estimated form factors induced by a decorrelation of local random
unitaries applied before and after the time evolution up to the Heisenberg
time $t_{H}$, with strength $\eta$ (see text). In (b), we display the
estimated SFF $(K)_{\text{dec}}$ (blue dots) and PSFF $(K_{A})_{\text{dec}}$
(blue squares) as function of time in a system subject to global polarization
with strength $p=0.03$ (see text). For this type of decoherence, rescaling
according to Eq. (38), allows to recover the SFF (green dots) and PSFF (green
squares) for unitary dynamics (red line).

First, we consider an imperfect implementation of our measurement protocols,
with errors arising from an erroneous decorrelation of the applied initial and
final local random unitaries. We model such imperfection as the effective
application of a unitary $u_{i}$ before and a unitary
$v_{i}=u^{\dagger}_{i}\exp(-i\eta h_{i})$ after the time evolution, with
$h_{i}$ being a local random Hermitian matrix sampled for each $i$
independently from the GUE [1]. While the case $\eta=0$ corresponds to the
ideal case, we display in Fig. 7(a) the average relative error
$\widehat{\epsilon_{\eta}}=1-\widehat{(K_{A}(t))_{\eta}}/\widehat{K_{A}(t)}$
of the estimated $\widehat{(K_{A}(t))_{\eta}}$ as a function of the error
strength $\eta$, obtained numerically from simulating many experimental runs.
We find $\widehat{\epsilon_{\eta}}$ increases approximately as $\eta^{2}$,
indicating a decrease of the estimated $\widehat{(K_{A}(t))_{\eta}}$.

Secondly, we consider that a measurement of the SFF and PSFF is affected by
decoherence acting during the dynamical evolution of the system. As shown in
the context of other randomized measurement protocols, one can correct the
effect of depolarization errors (or readout errors) based on a randomized
measurement of the purity [35, 37, 36, 11], which allows to extract the value
of the noise strength [37, 114]. Note that if the type of noise is a priori
unknown, one can also mitigate errors with randomized measurements. This is
done via a calibration step that allows to convert randomized measurements
into faithful ‘classical shadows’ estimations of the quantum state [115, 116,
117].

Here for concreteness, we consider a Floquet system with global
depolarization, acting at each time period $\tau$ with strength $p$, i.e. the
final state $\rho_{f}(t)$ at time $t=\tau n$, defined in Sec. I.3, is altered
to $\rho_{\text{dec}}(t)=\alpha_{n}\rho_{f}(t)+(1-\alpha_{n})\,\mathbb{1}/D$
with $\alpha_{n}=(1-p)^{n}$.

Thus, we obtain via our measurement protocol,

$\displaystyle(K_{A})_{\text{dec}}(t)=\alpha_{n}K_{A}(t)+\frac{1-\alpha_{n}}{D_{A}^{2}}\;.$
(35)

With increasing time $t=\tau n$, decoherence leads thus to a smaller measured
value $(K_{A})_{\text{dec}}(t)$ than the actual spectral form factor
$K_{A}(t)$ (see Fig. 7(b), blue dots and squares). However, if we know the
value of $p$, we can rescale our estimator of the SFF. For this purpose, we
can measure the purity of the time evolved state. The purity is,

$\displaystyle P_{n}$ $\displaystyle=$
$\displaystyle{\rm{Tr}}\left[\rho_{\rm{dec}}(t)^{2}\right]=\alpha_{n}^{2}+\frac{1-\alpha_{n}^{2}}{D}~{},$
(36)

which gives [37, 114],

$\alpha_{n}=\sqrt{\frac{DP_{n}-1}{D-1}}~{}.$ (37)

Thus, from a measurement of the purity $P_{n}$ at all times, we can find
$\alpha_{n}$ and rescale the erroneous PSFF (35) to obtain,

$(K_{A})_{\rm{res}}(t)=\frac{(K_{A})_{\rm{dec}}(t)-(1-\alpha_{n})/D_{A}^{2}}{\alpha_{n}}~{}.$
(38)

In Fig. 7(b), using the green color we present this rescaled SFF (using dots)
and PSFF (using squares). We note that using the rescaled (P)SFF (38) we
recover here the (P)SFF of the unitary dynamics (red curve).

In summary, while we have shown in this subsection that we can partially
correct for decoherence effects via independent measurements of decoherence
parameters, we emphasize that imperfections and decoherence discussed in this
section lead to a decay of the estimated $\widehat{K_{A}(t)}$. They, thus can
not cause a false positive detection of the ramp.

##  VI Conclusion and outlook

In this work, we have presented randomized measurement protocols to access the
statistics of energy eigenvalues and energy eigenstates of many-body quantum
systems in present day quantum simulators via (partial) spectral form factors.
The spectral form factor (SFF), $K(t)$ in Eq. (1), is known to be a key
diagnostic of many-body quantum chaos. In chaotic systems, it reveals
universal properties of energy eigenvalue statistics and possesses a
characteristic ramp-plateau structure (see Sec. I.1). In addition, we have
defined partial spectral form factors (PSFFs), $K_{A}(t)$ in Eq. (I.2), which
contain both the statistics of energy eigenvalues and eigenstates (see Sec.
I.2). PSFFs are natural restrictions of the SFF to subsystems
$A\subseteq\mathcal{S}$ of the full system $\mathcal{S}$, such that for
$A=\mathcal{S}$, PSFF and SFF coincide $K_{A=\mathcal{S}}(t)=K(t)$. Utilizing
random matrix theory and the eigenstate thermalization hypothesis (ETH), we
have shown in Sec. II, that PSFFs in generic chaotic quantum many-body systems
possess a characteristic shift-ramp-plateau structure [Eqs. (11) and (15)] and
reveal crucial differences between thermal and non-thermal eigenstates in the
sense of ETH. In Sec. III we investigated the PSFF numerically with examples
of many-body quantum models, discussing, in particular, differences between
chaotic and localized phases.

With our protocol to measure the SFF and PSFF in quantum simulation
experiments, we have extended the toolbox of randomized measurements to access
genuine properties of dynamical quantum evolution, without any reference to
the initial state or measured observable (see Secs. I.3, IV and V). We have
shown that our protocol gives simultaneous access to the SFF and PSFF, thereby
providing a unified testbed of the statistical properties of eigenvalues and
eigenstates. Our protocol can be directly implemented in state-of-the-art
quantum devices, based for instance on trapped ions [4, 6], Rydberg atoms [7]
and superconducting qubits [8, 19], providing crucial experimental tools for
the quantum simulation of many-body quantum chaos and the study of
thermalization in closed quantum systems.

Our work can be generalized in various directions. First, while we have
concentrated here on quantum simulators with local control realizing lattice
spin models, our protocol can be also realized in collective spin systems with
only global operations [118]. Second, while we have considered form factors
which are second-order functionals of the time evolution operators $T(t)$,
partial restrictions of higher-order form factors provide possibilities to
investigate thermalization of quantum many-body systems and emergent
randomness beyond second-order [119, 120]. To access such higher-order
(partial) form factors, our randomized measurement protocols could be readily
combined with the classical shadows framework [42]. Thirdly, we have focused
on determining the properties of unitary quantum dynamics. Beyond that, our
measurement protocol readily extends to the study of noisy quantum channels.
This includes applications in the field of verification and benchmarking of
quantum devices [89, 90, 91, 92, 93, 121, 122], as well as the investigation
of noise-induced quantum many-body phenomena such as entanglement phase
transitions [123, 124, 125, 46]. In addition to the directions listed above,
it will be interesting to explore the PSFF from an analytical perspective
analogous to Ref. [80] to study the physics of thermalization and entanglement
in Hamiltonian many-body systems as well as in quantum gravity, where there
have recently been path integral derivations of the SFF [63].

###### Acknowledgements.

We thank Mikhail Baranov, Amos Chan, Manoj K. Joshi, Barbara Kraus, Rohan
Poojary, Lukas Sieberer and Denis Vasilyev for valuable discussions. Work in
Innsbruck has been supported by the European Union’s Horizon 2020 research and
innovation programme under Grant Agreement No. 817482 (Pasquans) and No.
731473 (QuantERA via QT-FLAG), by the Austrian Science Foundation (FWF, P
32597 N), by the Simons Collaboration on Ultra-Quantum Matter, which is a
grant from the Simons Foundation (651440, P.Z.), and by LASCEM by AFOSR No.
64896-PH-QC. A.E. acknowledges funding by the German National Academy of
Sciences Leopoldina under the grant number LPDS 2021-02. BV acknowledges
funding from the French National Research Agency (ANR-20-CE47-0005, JCJC
project QRand). A.V. and V.G. were supported by US-ARO Contract
No.W911NF1310172, NSF DMR-2037158, and the Simons Foundation.

##  Appendix A Spectral form factor in Wigner-Dyson random matrix ensembles

In this appendix, we review the definition and essential properties of the
Wigner-Dyson random matrix ensembles. Further, we recall the expressions of
the SFF for Hamiltonian and Floquet dynamics modeled with random matrices from
these ensembles.

The Wigner-Dyson ensembles are standard distributions of random matrices used
to model some of the properties of energy or quasi-energy eigenvalues and
eigenstates of chaotic Hamiltonian and Floquet systems [20, 21, 24, 1]. We
work with two classes of the Wigner-Dyson ensembles - the unitary (U) class
for systems that are not time reversal invariant, and the orthogonal (O) class
for some systems with time-reversal invariance (the symplectic (S) class
applies to other systems with time-reversal invariance, but is not relevant
for our examples). We note in particular that nonconventional time-reversal
symmetries should also be considered [1] e.g. invariance under complex
conjugation in some basis (which corresponds to the orthogonal class). Each
class is characterized by a symmetry group comprised of the corresponding set
of similarity transformations (i.e. all unitary or orthogonal
transformations).

For Hamiltonian systems with time evolution operator $T(t)=\exp(-iHt)$, it is
conventional to choose the Gaussian Unitary Ensemble (GUE) of Hermitian
matrices or the Gaussian Orthogonal Ensemble (GOE) of real symmetric matrices
to represent the Hamiltonian $H$ of the appropriate class. In the case of
periodically driven Floquet dynamics with time-evolution operator $T(t=\tau
n)=V^{n},n\in\mathbb{N}$, where $V$ is the unitary Floquet operator
corresponding to a time period $\tau$, the appropriate representative
ensembles for $V$ are the Circular Unitary Ensemble (CUE) of unitary matrices
and the Circular Orthogonal Ensemble (COE) of symmetric unitary matrices.
These ensembles accurately model the local eigenvalue correlations of the
corresponding systems (but not necessarily global eigenvalue features larger
than the inverse Thouless time scale [27, 64] e.g. the smoothened density of
states), and describe an idealization of the eigenstate distribution (which is
generalized by ETH [27, 29]). But for the special case of chaotic Floquet
systems, the eigenstate distribution is seen to be in close agreement with the
Wigner-Dyson ensembles [81, 82, 83, 33, 84, 80].

For these random matrix models, the spectral form factor can be calculated
analytically (see for instance Ref. [97]). For completeness, we recall the
well-known expressions here. For Hamiltonians $H$ from GUE or GOE, one finds

GUE model

$\displaystyle K(t)$
$\displaystyle=r(t)^{2}+\frac{1}{D}\begin{cases}\frac{t}{t_{H}}&\text{ for
}0<t\leq t_{H},\\\ 1&\text{ for }t>t_{H},\end{cases}$ (39)

GOE model

$\displaystyle K(t)$
$\displaystyle=r(t)^{2}+\\!\frac{1}{D}\begin{cases}2\frac{t}{t_{H}}-\frac{t}{t_{H}}\log{\left(1+2\frac{t}{t_{H}}\right)}\\!\\!\\!&\\!\\!\text{
for }0<t\leq t_{H},\\\
2-\frac{t}{t_{H}}\log{\left(\frac{2t+t_{H}}{2t-t_{H}}\right)}\\!\\!\\!&\\!\\!\text{
for }t>t_{H},\end{cases}$ (40)

where $r(t)=t_{H}J_{1}(4Dt/t_{H})/(2Dt)$ with $J_{1}$ denoting the Bessel’s
function of the first kind. The Heisenberg time $t_{H}$, connected to the
inverse spacing of adjacent energy levels, depends on the width of the
Gaussian distribution of the matrix elements and marks the onset time of the
plateau of the SFF. For the results presented in Sec. (III), we fix it
numerically, by matching plateau onset times for the Hamiltonian Eq. (19) and
the GOE model.

For the Floquet operators $V$ from CUE or COE, one finds

CUE model

$\displaystyle K(t)=\frac{1}{D}\begin{cases}\frac{t}{t_{H}},&\text{for
}0<t\leq t_{H},\\\ 1,&\text{for }t>t_{H},\end{cases}$ (41)

COE model

$K(t)=\frac{1}{D}\begin{cases}\frac{2t}{t_{H}}-\frac{t}{t_{H}}\log\left(1+2\frac{t}{t_{H}}\right)&\text{
for }0<t\leq t_{H},\\\
2-\frac{t}{t_{H}}\log\left(\frac{2t/t_{H}+1}{2t/t_{H}-1}\right)&\text{ for
}t>t_{H},\end{cases}$ (42)

Here, $t_{H}=D\tau$ with $\tau$ to be identified with the period of the
Floquet system to be modeled.

##  Appendix B Partial spectral form factor in Wigner-Dyson random matrix
ensembles

In this section, we derive the functional form of the partial spectral form
factors, discussed in Sec. II, for Hamiltonian dynamics (Floquet dynamics)
modeled with the Wigner-Dyson random matrix ensembles GUE, GOE (CUE, COE), as
introduced in App. A.

Let $\mathcal{S}$ be a quantum system with Hilbert space $\mathcal{H}$ of
dimension $D$, and $A\subseteq\mathcal{S}$ a subsystem with dimension $D_{A}$.
Its complement is denoted with $B$ with dimension ${{D}_{B}}$. As discussed in
App. A, we consider

  * •

Hamiltonian dynamics $T(t)=\exp(-iHt)$ with $H$ sampled from the GUE and GOE,
respectively.

  * •

Floquet dynamics with $T(t=\tau n)=V^{n}$ for $n\in\mathbb{N}$ with $V$
sampled from the CUE and COE, respectively.

We can rewrite $T(t)=YD(t)Y^{\dagger}$ with
$D(t)=\text{diag}(e^{-iE_{1}t},\dots,e^{-iE_{D}t})$, the diagonal matrix of
eigenvalues of $T(t)$ and $Y=(y_{1},\dots,y_{D})$ the unitary (GUE, CUE) or
orthogonal (GOE, COE) matrix of eigenvectors of $H$ or $V$. Crucially, we note
that all time-dependence is contained in the diagonal matrix $D(t)$. In the
following, we rely on the fact:

######  Fact 1.

For $H$ from GUE or GOE ($V$ from CUE or COE), the distribution of the
eigenvectors of $H$ ($V$) is independent of the distribution of eigenvalues of
$H$ ($V$). Further, $Y=(y_{1},\dots,y_{D})$ is distributed according to the
Haar measure on the group of unitary matrices $U(D)$ (for GUE, CUE) and the
group of orthogonal matrices $O(D)$ (for GOE, COE).

###### Proof.

This fact relies only on the invariance of the random matrix ensembles under
unitary (GUE, CUE) and orthogonal transformations (GOE, COE). For GUE and GOE,
a proof is given in Ref. [126], Corollary 2.5.4. It generalizes directly to
CUE and COE. ∎

Using this fact, we can carry out the average over eigenvectors in Eq. (I.2)
explicitly (see next subsection). With the identification
$K(t)=D^{-2}\overline{|\textrm{Tr}\left[D(t)\right]|^{2}}$, we find

$K_{A}(t)={c_{A}^{(1)}}+c_{A}^{(2)}K(t)~{},$ (43)

where $\text{for}~{}H\in\text{GUE}\;,\;V\in\text{CUE}$,

$\displaystyle c_{A}^{(1)}$
$\displaystyle=\frac{\left.{{{D}_{B}}}^{2}-1\right.}{{D_{A}}^{2}{{{D}_{B}}}^{2}-1}~{}~{};~{}~{}c_{A}^{(2)}=\frac{{{D}_{B}}^{2}\left({D_{A}}^{2}-1\right)}{{D_{A}}^{2}{{{D}_{B}}}^{2}-1}~{},$
(44)

and $\text{for}~{}H\in\text{GOE}\;,\;V\in\text{COE}$,

$\displaystyle c_{A}^{(1)}$
$\displaystyle=\frac{\left({{D}_{B}}^{2}+{{D}_{B}}-2\right)}{({D_{A}}{{{D}_{B}}}-1)({D_{A}}{{{D}_{B}}}+2)}~{};$
$\displaystyle c_{A}^{(2)}=$
$\displaystyle\frac{{{D}_{B}}\left(D_{A}{{D}_{B}}+{{D}_{B}}+1\right)\left({D_{A}}-1\right)}{({D_{A}}{{{D}_{B}}}-1)({D_{A}}{{{D}_{B}}}+2)}~{}.$
(45)

In particular, $K_{A=\mathcal{S}}(t)=K(t)$ for $D_{A}=D,{{D}_{B}}=1$ and
$K_{A=\emptyset}(t)=1$ for $D_{A}=1,{{D}_{B}}=D$ holds, as expected.

##### Relation to average purity and overlap:

For Hamiltonian $T(t)=\exp(-iHt)$ or Floquet dynamics $T(t=n\tau)=V^{n}$, we
can rewrite the PSFF in terms of the (quasi-) energy eigenvalues and (quasi-)
energy eigenstates [see Eq. (I.2)]. For Hamiltonians $H$ (Floquet operators
$V$) from the Wigner-Dyson random matrix ensembles we can use then fact 1 to
obtain the PSFF in terms of the average purity $\mathcal{P}_{B}$ of reduced
eigenstates and average overlap of distinct reduced eigenstates
$\mathcal{Q}_{B}$ [see Sec. II, in particular Eq. (11)]. Comparing Eq. (11)
with Eq. (43) we find that

$\displaystyle
c_{A}^{(1)}=\frac{\mathcal{P}_{B}-\mathcal{Q}_{B}}{D_{A}}\quad\text{and}\quad
c_{A}^{(2)}=D_{B}\mathcal{Q}_{B}~{}.$ (46)

Using this, Eqs. (44) and Eqs. (45), we obtain Eq. (12) (for GUE, CUE) and the
corresponding expressions for the orthogonal ensembles (GOE, COE),
respectively.

### Proof of Eqs. (43), (44) and (45)

We denote the basis of $\mathcal{H}$ consisting of eigenvectors of $T(t)$ with
$\ket{i}$ ($i=1,\dots,D$). Furthermore, we fix an arbitrary product basis of
$\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ as $\ket{a,b}$ with
$a=1,\dots,D_{A}$ and $b=1,\dots,{{D}_{B}}$. With $T(t)=YD(t)Y^{\dagger}$, we
rewrite Eq. (I.2) in these bases. Using the independence of eigenvalues and
eigenvectors (Fact 1), we find

$\displaystyle DD_{A}K_{A}(t)$
$\displaystyle=\overline{\textrm{Tr}_{B}\left[\textrm{Tr}_{A}\left[T(t)\right]\textrm{Tr}_{A}\left[(T(t)^{\dagger}\right]\right]}$
(47) $\displaystyle=$
$\displaystyle\overline{Y_{(a_{1},b_{1}),i_{1}}(Y^{\dagger})_{i_{1},(a_{1},b_{2})}Y_{(a_{2},b_{2}),i_{2}}(Y^{\dagger})_{i_{2},(a_{2},b_{1})}}\times$
$\displaystyle\overline{D(t)_{i_{1},i_{1}}(D(t)^{\dagger})_{i_{2},i_{2}}}~{},$
(48)

where summation over repeated indices is understood. The ensemble average over
the matrix elements of $Y$ can be carried out using the Weingarten calculus on
the unitary group (GUE and CUE) and orthogonal group (GOE and COE),
respectively.

Figure 8:  Diagrammatic evalution of Eq. (48) for $Y\in U(D)$ (GUE, CUE case).
To perform the average over eigenvectors (green line), we remove the boxes $Y$
and connect white decorations of $Y$ (rhombi) with white decorations of
$Y^{*}$ (rhombi) and black decorations of $Y$ (circles, squares) with black
decorations of $Y^{*}$ (circles, squares) in all possible ways, corresponding
to the pair partitions $\mathbf{m},\mathbf{n}\in\mathcal{M}^{U}(4)$ [127,
108]. Summing over the resulting diagrams, weighted with corresponding value
of the Weingarten function, yields Eq. (44). In all diagrams, each blue loop
contributes a factor $D_{A}$, each red loop a factor ${{D}_{B}}$. Figure 9:
Diagrammatic evalution of Eq. (48) for $Y\in O(D)$ (GOE, COE case). To perform
the average over eigenvectors (green line), we remove the boxes $Y$ and
connect white decorations (rhombi) with white decorations (rhombi) and black
decorations (circles, squares) with black decorations of same type (circles,
squares) in all possible ways, corresponding to all pair partitions
$\mathbf{m},\mathbf{n}\in\mathcal{M}^{O}(4)$ [127, 108]. Summing over the
resulting diagrams, weighted with corresponding value of the Weingarten
function, yields Eq. (45). In all diagrams, each blue loop contributes a
factor $D_{A}$, each red loop a factor ${{D}_{B}}$.

The Weingarten calculus for the unitary group and for orthogonal group can be
formulated in terms of pair partitions, defined as follows.

######  Definition 1 (Pair partitions).

For $n\in\mathbb{N}$, a) we denote with $\mathcal{M}^{O}(2n)$ the set of all
pair partitions of $\\{1,\dots,2n\\}$, partitioning $\\{1,\dots,2n\\}$ into
$n$ distinct pairs. Then, each pair partition
$\mathbf{m}\in\mathcal{M}^{O}(2n)$ can be uniquely expressed as

$\displaystyle\left\\{\\{\mathbf{m}(1),\mathbf{m}(2)\\},\dots\\{\mathbf{m}(2n-1),\mathbf{m}(2n)\\}\right\\}$
(49)

with $\mathbf{m}(1)<\mathbf{m}(3)<\dots<\mathbf{m}(2n-1)$ and
$\mathbf{m}(2i-1)<\mathbf{m}(2i)$ for all $i\in\\{1,\dots n\\}$.
b) we denote with $\mathcal{M}^{U}(2n)\subseteq\mathcal{M}^{O}(2n)$ the set of
all pair partitions of $\\{1,\dots,2n\\}$ which pair elements in
$\left\\{1,\dots,n\right\\}$ with elements $\left\\{n+1,\dots,2n\right\\}$.
Then, each partition $\mathbf{m}\in\mathcal{M}^{U}(2n)$ can be uniquely
expressed as

$\displaystyle\left\\{\\{\mathbf{m}(1),\mathbf{m}(2)\\},\dots\\{\mathbf{m}(2n-1),\mathbf{m}(2n)\\}\right\\}$
(50)

with $\mathbf{m}(1)<\mathbf{m}(3)<\dots<\mathbf{m}(2n-1)$ and
$\mathbf{m}(2i-1)\in\left\\{1,\dots,n\right\\}$ and
$\mathbf{m}(2i)\in\left\\{n+1,\dots,2n\right\\}$ for all $i\in\\{1,\dots
n\\}$.

The following fact is shown in Ref. [128].

######  Fact 2 (Weingarten calculus).

(i) Let $Y$ be distributed according to the Haar measure on the orthogonal
group $O(D)$. With indices $i_{1},\dots,i_{2n}$ and $j_{1},\dots,j_{2n}$ in
$\left\\{1,\dots,D\right\\}$ it holds

$\displaystyle\int_{Y\in O(D)}Y_{i_{1},j_{1}}\cdots
Y_{i_{2n},j_{2n}}\text{d}Y=$
$\displaystyle\smashoperator[]{\sum_{\mathbf{m},\mathbf{n}\in\mathcal{M}^{O}(2n)}^{}}\;\textup{Wg}^{O(D)}(\mathbf{m},\mathbf{n})\prod_{k=1}^{n}\delta_{i_{\mathbf{m}(2k-1)},i_{\mathbf{m}(2k)}}\delta_{j_{\mathbf{m}(2k-1)},j_{\mathbf{m}(2k)}}$
(51)

with $\mathcal{M}^{O}(2n)$ the set of all pair partitions on
$\left\\{1,2,\dots,2n\right\\}$ and $\textup{Wg}^{O(D)}$ the Weingarten
function on the orthogonal group $O(D)$.
(ii) Let $Y$ be distributed according to the Haar measure on the unitary group
$U(D)$. With indices $i_{1},\dots,i_{2n}$ and $j_{1},\dots,j_{2n}$ in
$\left\\{1,\dots,D\right\\}$ it holds

$\displaystyle\int_{Y\in U(D)}Y_{i_{1},j_{1}}\cdots
Y_{i_{n},j_{n}}\,Y^{*}_{i_{n+1},j_{n+1}}Y^{*}_{i_{2n},j_{2n}}\text{d}Y=$
$\displaystyle\smashoperator[]{\sum_{\mathbf{m},\mathbf{n}\in\mathcal{M}^{U}(2n)}^{}}\;\textup{Wg}^{U(D)}(\mathbf{m},\mathbf{n})\prod_{k=1}^{n}\delta_{i_{\mathbf{m}(2k-1)},i_{\mathbf{m}(2k)}}\delta_{j_{\mathbf{m}(2k-1)},j_{\mathbf{m}(2k)}}$
(52)

with $\mathcal{M}^{U}(2n)\subsetneq\mathcal{M}^{O}(2n)$ the set of all pair
partitions on $\left\\{1,2,\dots,2n\right\\}$ which pair elements in
$\left\\{1,\dots,n\right\\}$ with elements $\left\\{n+1,\dots,2n\right\\}$ and
$\textup{Wg}^{U(D)}$ the Weingarten function on the unitary group $U(D)$.

In our case, we are only interested in the case $n=2$. As shown in Ref. [128],
when $\mathbf{m},\mathbf{n}\in\mathcal{M}^{O}(4)$ and $D\geq 2$,

$\displaystyle w^{O}_{\text{eq}}$
$\displaystyle\equiv\text{Wg}^{O(D)}(\mathbf{m},\mathbf{n})=\frac{D+1}{D(D+2)(D-1)}\;\text{for}\;\mathbf{m}=\mathbf{n}$
$\displaystyle w^{O}_{\text{neq}}$
$\displaystyle\equiv\text{Wg}^{O(D)}(\mathbf{m},\mathbf{n})=\frac{-1}{D(D+2)(D-1)}\;\text{for}\;\mathbf{m}\neq\mathbf{n}~{}.$
(53)

Furthermore, it holds for $\mathbf{m},\mathbf{n}\in\mathcal{M}^{U}(4)$ and
$D\geq 2$

$\displaystyle w^{U}_{\text{eq}}$
$\displaystyle\equiv\text{Wg}^{U(D)}(\mathbf{m},\mathbf{n})=\frac{D}{D(D^{2}-1)}\quad\text{for}\quad\mathbf{m}=\mathbf{n}$
$\displaystyle w^{U}_{\text{neq}}$
$\displaystyle\equiv\text{Wg}^{U(D)}(\mathbf{m},\mathbf{n})=\frac{-1}{D(D^{2}-1)}\quad\text{for}\quad\mathbf{m}\neq\mathbf{n}~{}.$
(54)

Using Fact 2 and these expressions, we can perform the average over
eigenvector elements in Eq. (48) explicitly. This is most easily performed
diagrammatically and shown in Figs. 8 and 9.

##  Appendix C Partial spectral form factor in general chaotic systems

Here, we derive the typical behavior of the PSFF for ensembles of chaotic
systems, more general than random matrix ensembles, as considered in Sec. II.2
of the main text. As in Eq. (13), we decompose the reduced density matrix into
a pure trace, a traceless smooth part and a traceless fluctuating part,
$\rho_{B}(E)=D_{B}^{-1}\mathbb{1}+\Delta\rho_{B}(E)+\delta\rho_{B}(E)$. For
the smooth part, we assume that there exists an extrapolation of each matrix
element to a continuous energy variable such that for some (as yet
unspecified) time $t_{\rho}\ll O(D)$,

$(\Delta\tilde{\rho}_{B}(t))_{jk}\equiv\int\text{d}E\
e^{-iEt}(\Delta{\rho_{B}}(E))_{jk}=0,\ \forall\ \lvert t\rvert>t_{\rho}~{}.$
(55)

The remaining energy dependent part of $\rho_{B}(E)$ i.e. the part that
oscillates rapidly and has no low frequency Fourier component (on
extrapolation to continuous energy) will be taken to be the fluctuating part,

$(\delta\tilde{\rho}_{B}(t))_{jk}\equiv\int\text{d}E\
e^{-iEt}(\delta{\rho_{B}}(E))_{jk}=0,\ \forall\ \lvert t\rvert\leq
t_{\rho}~{}.$ (56)

Up to this point, such a decomposition is always possible. We will
additionally take $t_{\rho}$ to be set by the scale of randomization in the
ensemble discussed in Sec. II.2, so that the fluctuating part can be
identified as the part that is completely randomized in the ensemble. We note
that the smooth part may fluctuate between different ensemble realizations,
but can not be randomized in the same sense as the fluctuating part as it is
roughly constant within an energy window of size $t_{\rho}^{-1}$. Similarly,
we will not require randomization of the correlators of $\delta\rho_{B}(E)$
between energies further apart than $\sim t_{\rho}^{-1}$, for which the
correlator may have to be nonvanishing to maintain zero Fourier component of
the fluctuating part at $t\leq t_{\rho}$.

To understand the effect of this decomposition in the PSFF, we will first
perform a prototype calculation with simpler notation. Consider two functions
$f(E)$ and $g(E)$ of a continuous variable $E$, with respective Fourier
transforms $\tilde{f}(t)$ and $\tilde{g}(t)$, both of which potentially vary
over different realizations of the ensemble. We will eventually associate
these functions with (components of) the different parts of the reduced
density matrices of the energy eigenstates. Define the quantity,

$\displaystyle
F(t)=\frac{1}{D^{2}}\overline{\sum_{j,k}e^{i(E_{j}-E_{k})t}f(E_{j})g^{\ast}(E_{k})}$
$\displaystyle=\frac{1}{D^{2}}\int\frac{\text{d}t_{l}}{2\pi}\int\frac{\text{d}t_{r}}{2\pi}\
\overline{\tilde{f}(t_{l})\tilde{g}^{\ast}(t_{r})\sum_{j,k}e^{iE_{j}(t+t_{l})-iE_{k}(t+t_{r})}}~{}.$
(57)

Now, it is convenient to define an ensemble-averaged unequal time SFF
$K(t_{1},t_{2})=D^{-2}\overline{\sum_{j,k}e^{iE_{j}t_{1}-iE_{k}t_{2}}}$, which
reduces to $K(t)$ at equal times $t_{1}=t_{2}=t$. The sum of phases
$D^{-2}\sum_{j,k}e^{iE_{j}(t+t_{l})-iE_{k}(t+t_{r})}$ in Eq. (57) would
fluctuate strongly over different ensemble realizations at large $t_{1},t_{2}$
corresponding to fluctuations of the positions of energy levels, much like the
SFF without ensemble averaging [59]; if we assume the ensemble is such that
these fluctuations are not correlated with those of $f$ and $g$ (i.e. the
reduced energy eigenstates), we can perform the ensemble average over the sum
of phases independently, allowing us to formally replace it with
$K(t+t_{l},t+t_{r})$,

$F(t)=\int\frac{\text{d}t_{l}}{2\pi}\int\frac{\text{d}t_{r}}{2\pi}\
K(t+t_{l},t+t_{r})\overline{\tilde{f}(t_{l})\tilde{g}^{\ast}(t_{r})}~{}.$ (58)

For instance, in a fully chaotic system as we will soon specialize to, this
assumption can be justified by considering the energy eigenstates in an
ensemble realization as sufficiently random superpositions of those of another
ensemble realization (in the spirit of Refs. [2, 94, 95, 96]), which should
then be uncorrelated with the precise positions of the energy levels.

To simplify Eq. (58) further, we need to know the form of $K(t_{1},t_{2})$.
For mathematical simplicity, we assume (fully chaotic) level statistics in the
unitary Wigner-Dyson class. The ensemble-averaged two level correlation
function for nearby energy levels $E_{j}$,$E_{k}$ (closer than $\sim t_{\rm
Th}^{-1}$) in this class takes the universal form [24, 1, 97],

$\displaystyle\overline{\delta\left(E+\frac{\omega}{2}-E_{j}\right)\delta\left(E-\frac{\omega}{2}-E_{k}\right)}$
$\displaystyle=\Omega^{2}(E)\left\\{1+\frac{\delta(\omega)}{\Omega(E)}-\operatorname{sinc}^{2}\left[\omega\pi\Omega(E)\right]\right\\}~{},$
(59)

where $\Omega(E)$ is the smoothened (continuous and ensemble-averaged) density
of states, whose Fourier transform satisfies $\tilde{\Omega}(t\gg t_{\rm
Th})\approx 0$. The ensemble averaged sum over $E_{j},E_{k}$ in the definition
of $K(t)$ can then be replaced by an integral weighted by the two level
correlation in Eq. (59). Using methods analogous to the calculation of $K(t)$
for this correlation function in Ref. [97], we obtain the following late time
behavior for $t_{1},t_{2}\gg t_{\rm Th}$,

$\displaystyle
K(t_{1},t_{2})=\frac{1}{D^{2}}\begin{dcases}\tilde{\Omega}(\tau_{12}),&T_{12}>2\pi\Omega(E)\
\forall E,\\\ \frac{\lvert
T_{12}\rvert}{\beta\pi}\tilde{\Theta}_{\Omega}(\tau_{12}),&T_{12}<2\pi\Omega(E)\
\forall E,\end{dcases}$ (60)

where $\beta=2$ for the unitary Wigner-Dyson class, and we have introduced the
shorthand symbols $T_{12}=(t_{1}+t_{2})/2$, $\tau_{12}=t_{2}-t_{1}$.
$\tilde{\Theta}_{\Omega}(t)$ is the Fourier transform of the unit step
function $\Theta(\Omega(E))$, the latter being $1$ where $\Omega(E)>0$ and
zero elsewhere. Essentially, the unequal time SFF is generally negligible for
(large) unequal times, with a small spread around $t_{1}=t_{2}$ determined by
the variation of the density of states; as noted earlier, it reduces to the
SFF at precisely equal times. We also identify $2\pi\Omega(E)$ with the
Heisenberg time $t_{H}$, assuming that $\Omega(E)$ is at least of the same
order of magnitude throughout the spectrum. In the orthogonal and symplectic
Wigner-Dyson classes, there are significant corrections (relative to the
unitary class) to the form of the equal time SFF $K(t)$ near $t\sim t_{H}$.
But for $t\ll t_{H}$, virtually the same results hold with $\beta=1$ for the
orthogonal class and $\beta=4$ for the symplectic class [97] (of course, the
plateau behavior for $t\gg t_{H}$ is generally independent of such specifics).
Analogously, we expect similar replacements (the appropriate value of $\beta$,
and focusing on the $T_{12}\gg t_{H}$ and $T_{12}\ll t_{H}$ regimes) to work
for the unequal time SFF in Eq. (60) as well. With this expectation, we write

$K(t_{1},t_{2})=\frac{1}{D^{2}}\begin{dcases}\tilde{\Omega}(\tau_{12}),&T_{12}\gg
t_{H},\\\ \frac{\lvert
T_{12}\rvert}{\beta\pi}\tilde{\Theta}_{\Omega}(\tau_{12}),&T_{12}\ll
t_{H},\end{dcases}$ (61)

for $t_{1},t_{2}\gg t_{\rm Th}$ in any Wigner-Dyson symmetry class.

Using the decomposition of $\rho_{B}(E)$ with these definitions then gives
several terms for $K_{A}(t)$ of the form of Eq. (57), where $f$ and $g$
independently go over each of $D_{B}^{-1}$, $\Delta\rho_{B}$ and
$\delta\rho_{B}$, with an additional trace of the product over the $B$
subspace. Now, we will argue that all cross terms with $f\neq g$ may be taken
to vanish. When $f=D_{B}^{-1}$, the overlap becomes
$\textrm{Tr}_{B}\left[fg\right]=D_{B}^{-1}\textrm{Tr}_{B}\left[g\right]$,
which is zero when $g=\Delta\rho_{B},\delta\rho_{B}$, which are both
traceless. When say, $f$ is $\Delta\rho_{B}$ and $g$ is $\delta\rho_{B}$, the
cross term vanishes due to the assumption that ensemble averaging randomizes
$\delta\rho_{B}$.

Dropping the cross terms for the above reasons gives the form of Eq. (14) in
the main text, $K_{A}(t)=K(t)+\Delta K_{A}(t)+\delta K_{A}(t)$, where $K(t)$
is the full SFF, and

$\displaystyle\Delta K_{A}(t)$
$\displaystyle=\sum_{j,k}\frac{\overline{e^{i(E_{j}-E_{k})t}\textrm{Tr}_{B}\left[\Delta\rho_{B}(E_{j})\Delta\rho_{B}(E_{k})\right]}}{DD_{A}},$
(62) $\displaystyle\delta K_{A}(t)$
$\displaystyle=\sum_{j,k}\frac{\overline{e^{i(E_{j}-E_{k})t}\textrm{Tr}_{B}\left[\delta\rho_{B}(E_{j})\delta\rho_{B}(E_{k})\right]}}{DD_{A}}.$
(63)

In the main text, it is argued that $\delta K_{A}(t\gg t_{\rho})$ amounts to a
constant shift after ensemble averaging due to the randomization of
$\delta\rho_{B}(E)$. Here, we will complete the evaluation of $\Delta
K_{A}(t)$ using the prototype Eq. (58) with $f=g=(\Delta\rho_{B})_{ab}$ and
the expression in Eq. (61) with $t_{1}=t+t_{l}$, $t_{2}=t+t_{r}$. As the
definition of $\Delta\rho_{B}$ sets $t_{l},t_{r}<t_{\rho}$, we have $\lvert
T_{12}\rvert=\lvert t\rvert+\operatorname{sgn}(t)(t_{l}+t_{r})/2$ at large
times (i.e. $t\gg t_{\rm Th},t_{\rho}$). For $t\ll t_{H}$ in this regime, this
gives,

$\displaystyle\Delta K(t:t_{\rm Th},$ $\displaystyle t_{\rho}\ll t\ll
t_{H})=\frac{1}{DD_{A}}\int\frac{\text{d}t_{l}}{2\pi}\int\frac{\text{d}t_{r}}{2\pi}$
$\displaystyle\left[\frac{1}{\beta\pi}\left(\lvert
t\rvert+\operatorname{sgn}(t)\frac{t_{l}+t_{r}}{2}\right)\tilde{\Theta}_{\Omega}(t_{l}-t_{r})\right.$
$\displaystyle\left.\overline{\left(\sum_{a,b}\left(\Delta\tilde{\rho}_{B}(t_{l})\right)_{ab}\left(\Delta\tilde{\rho}_{B}^{\ast}(t_{r})\right)_{ab}\right)}\right].$
(64)

The Hermiticity of $\Delta\rho_{B}$ implies that
$\left(\Delta\tilde{\rho}_{B}(-t)\right)_{ab}=\left(\Delta\tilde{\rho}_{B}^{\ast}(t)\right)_{ba}$.
Consequently, making the integration variable transformation $t_{l}\to-t_{r}$,
$t_{r}\to-t_{l}$ in Eq. (64), we see that inside the parentheses in the second
line the $\lvert t\rvert$ term is unaltered but the $\operatorname{sgn}(t)$
term transforms to its negative, while all factors outside the parentheses
remain unaltered. It follows that the contribution from the
$\operatorname{sgn}(t)$ term actually evaluates to zero, leaving only a linear
ramp term from $\lvert t\rvert$. For $t\gg t_{H}$, we directly obtain only a
plateau contribution. Now, it is straightforward to Fourier transform back to
the energy variable $E$,

$\displaystyle\Delta K_{A}(t\gg t_{\rm Th},t_{\rho})$
$\displaystyle=\frac{1}{DD_{A}}\int\text{d}E\
\begin{dcases}\overline{\Omega(E)\textrm{Tr}_{B}\left[\Delta\rho_{B}^{2}(E)\right]}~{},&\
t\gg t_{H},\\\
\frac{t}{\beta\pi}\overline{\Theta(\Omega(E))\textrm{Tr}_{B}\left[\Delta\rho_{B}^{2}(E)\right]}~{},&\
t\ll t_{H}.\end{dcases}$ (65)

For ease of interpretation, we can convert $E$ back to a discrete energy
variable from its present continuous form via the following correspondence
relations for sums over energy levels: $\sum_{i}\leftrightarrow\int\text{d}E\
\Omega(E)$ and $\sum_{i}\Omega^{-1}(E_{i})\leftrightarrow\int\text{d}E\
\Theta(\Omega(E))$, which become equalities on ensemble averaging. Then we get
the expression,

$\displaystyle\Delta K_{A}(t\gg t_{\rm Th},t_{\rho})$
$\displaystyle=\frac{1}{DD_{A}}\
\begin{dcases}\overline{\sum_{i}\textrm{Tr}_{B}\left[\Delta\rho_{B}^{2}(E_{i})\right]}~{},&\
t\gg t_{H},\\\
\frac{t}{\beta\pi}\overline{\sum_{i}\Omega^{-1}(E_{i})\textrm{Tr}_{B}\left[\Delta\rho_{B}^{2}(E_{i})\right]}~{},&\
t\ll t_{H}.\end{dcases}$ (66)

Together with the expression for the full SFF [$t_{1}=t_{2}$ in Eq. (61)] and
the constant contribution from the fluctuating part, this directly leads to
Eq. (15) in the main text.

##  Appendix D Constraints from eigenstate thermalization

In this Appendix, we discuss the constraints on the spectrum and ensemble
averaged PSFF parameters, $\mathcal{P}_{B}$ (purity of reduced density
matrices), $\delta\mathcal{P}_{B}$ (fluctuating part) and
$\Delta\mathcal{P}_{B}$ (smooth part), as measures of the extent of
delocalization and thermalization of energy eigenstates. In App. D.1, we
discuss these constraints based on a qualitative picture of subsystem ETH,
paying particular attention to thermalization as a distinct phenomenon from
delocalization. We justify this qualitative picture in the subsequent section,
first in terms of a version of the original conjecture of subsystem ETH [29]
for fully delocalized states in App. D.2.1, and argue for its extension to
eigenstates of arbitrary delocalization in App. D.2.2.

###  D.1 PSFF as a probe of thermalization and delocalization

We begin with a qualitative discussion of thermalization (in the sense of
subsystem ETH) and delocalization. We work in a ‘physical basis’ - one whose
basis vectors are close to pure states in most physically accessible (e.g.
local [30]) subsystems, such as a product basis of qubits. Thermalization then
corresponds to a significant overlap of the macroscopic features of
eigenstates of nearby energies whose individual components are sufficiently
random (and therefore, macroscopically similar), whereas non-thermal behavior
is seen when nearby eigenstates do not have a large overlap. This is to be
distinguished from the extent of delocalization of an eigenstate, which is the
number of bases states $\ell\leq D$ that it has a significant probability of
being found in.

It is useful to introduce an effective dimension $D_{A}^{\rm{eff}}\leq
D_{A},\ell$ of the Hilbert space of subsystem $A$, corresponding to the
typical number of degrees of freedom of subsystem $A$ over which the
eigenstate is delocalized within its support in the physical basis. In
particular, $D_{A}^{\rm{eff}}=D_{A}$ if the eigenstates appear completely
delocalized over subsystem $A$, and more generally $D_{A}^{\rm{eff}}$ is
typically larger for larger $D_{A}$ (up to $\ell$). For instance,
$D_{A}^{\rm{eff}}$ is a monotonically increasing function of $D_{A}$ when the
latter is varied by successively choosing larger subsystems $A$ containing the
previous one; additionally, it increases from $D_{A}^{\rm{eff}}=1$ for
$D_{A}=1$, to $D_{A}^{\rm{eff}}=\ell$ for $D_{A}=D$. We also use the notation
$O(x)$ to mean a non-negative number whose magnitude is at most of the order
of magnitude of $x$, to leading order when $x\gg 1$. In particular, we will
take $D\gg D_{A},D_{B}\gg 1$.

Assuming that $D_{A}^{\rm{eff}}$ is typical for $A$ throughout the spectrum,
the purity in subsystem $B$ satisfies,

$\mathcal{P}_{B}=(D_{A}^{\rm{eff}}/\ell)+O(D_{A}^{\rm{eff}}/\ell)+O(1/D_{A}^{\rm{eff}}),$
(67)

subject to
$\mathcal{P}_{B}\gtrsim(D_{A}^{\rm{eff}}/\ell),(1/D_{A}^{\rm{eff}})$. The
first two terms are due to the eigenstate being delocalized in subsystem $B$
with effective dimension $(\ell/D_{A}^{\rm{eff}})$, with the second term
containing larger scale variations of its components. We will call this, the
‘macroscopic’ contribution, which grows with $D_{A}^{\rm{eff}}$. The last term
is due to the randomness of the eigenstate components i.e. the ‘microscopic’
contribution, which decays with $D_{A}^{\rm{eff}}$ (and is also typically
bounded from below by $(1/D_{A}^{\rm{eff}})$). Being a linear combination of
the macroscopic and microscopic contribution, the purity shows an initial
decay with $D_{A}^{\rm{eff}}$ for small values of the latter, and eventually a
growth for larger values of $D_{A}^{\rm{eff}}\gtrsim\sqrt{\ell}$. Both
$D_{A}^{\rm{eff}}=1,\ell$ correspond to pure states with $\mathcal{P}_{B}=1$.

The parameters $\delta\mathcal{P}_{B}$, $\Delta\mathcal{P}_{B}$ satisfy the
following order-of-magnitude inequalities,

$\displaystyle\delta\mathcal{P}_{B}$ $\displaystyle\gtrsim
O(1/D_{A}^{\rm{eff}}),$ (68) $\displaystyle D_{B}^{-1}+\Delta\mathcal{P}_{B}$
$\displaystyle\lesssim(D_{A}^{\rm{eff}}/\ell)+O(D_{A}^{\rm{eff}}/\ell).$ (69)

The first inequality is the statement that the fluctuating part must include
at least the randomness of eigenstate components; the second says that the
smooth part or overlap of such eigenstates can at most contain all their
macroscopic features. They are also subject to the constraint
$\mathcal{P}_{B}=D_{B}^{-1}+\Delta\mathcal{P}_{B}+\delta\mathcal{P}_{B}$,
which can be interpreted in the present context as follows: the macroscopic
contribution to the purity must be distributed in some manner between the
smooth and fluctuating parts (with the exception of the maximally mixed part
$D_{B}^{-1}$); the microscopic contribution is however completely contained in
the fluctuating part.

According to ETH, the only difference between thermal eigenstates of nearby
energies is in their microscopic random fluctuations, with all their
macroscopic features completely contained in their overlap. This means that
the inequalities in Eqs. (68) and (69) are satisfied as equalities for thermal
eigenstates. In particular, $\delta\mathcal{P}_{B}$ can only decay with
increasing $D_{A}^{\rm{eff}}$ \- a fact that is responsible for the nearly
identical dynamics of observables in subsystem $B$ (for large $D_{A}$) in such
eigenstates. In contrast, non-thermal eigenstates have at least some of the
macroscopic contribution included in the fluctuating part, and therefore
satisfy Eqs. (68) and (69) much further from equality. In this case, the
macroscopic contribution to the fluctuating part may even show up as a growth
of $\delta\mathcal{P}_{B}$ with $D_{A}^{\rm{eff}}$ if the latter is
sufficiently large (analogous to the behavior of the purity), for choices of
subsystems where the incomplete overlap of neighboring eigenstates remains
‘visible’. At the same time, all eigenstates trivially satisfy
$\delta\mathcal{P}_{B}=\Delta\mathcal{P}_{B}=0$ for $D_{A}=D$.

We conclude that $\mathcal{P}_{B}$ is a measure of delocalization of
eigenstates, while $\delta\mathcal{P}_{B}$ and $\Delta\mathcal{P}_{B}$ are
probes of thermalization. Setting $\ell=D$ gives the results discussed in the
main text for chaotic systems with fully delocalized eigenstates (Sec.
II.2.2). For fully localized systems, $\ell=O(1)$ gives
$D_{A}^{\rm{eff}}=O(1)$, with $\mathcal{P}_{B}=O(1)$ and
$\delta\mathcal{P}_{B}=O(1)\lesssim(1-D_{B}^{-1})$, automatically implying a
lack of thermalization (Sec. II.3). Additionally, the same results hold when
the PSFF is defined only over a portion of the spectrum, where the parameters
merely become averages over that portion of the spectrum. This suggests that
such a filtered [64] PSFF can access equivalent information about the
properties of a smaller set of eigenstates of interest.

###  D.2 Subsystem ETH constraints

####  D.2.1 Fully delocalized eigenstates

Subsystem ETH [29] is a hypothesis concerning the behavior of energy
eigenstates in a chaotic system, applying in its original version to fully
delocalized eigenstates. It states that the eigenstates are of such a form as
to lead to the thermal behavior of all observables on subsystem $B$, when it
is a physically accessible subsystem - in the sense of diagonal and off-
diagonal ETH (e.g. as presented in the reviews [27, 28]). Denoting the
eigenstates by $\lvert E\rangle$, there are two statements of the hypothesis:
the diagonal statement stating that the reduced density matrix
$\rho_{B}(E)=\textrm{Tr}_{A}\left[\lvert E\rangle\langle E\rvert\right]$ is
close to some smooth density matrix $P_{B}(E)$ that does not vary rapidly with
energy, and the off-diagonal statement requiring the reduced transition
operators $q_{B}(E_{1},E_{2})=\textrm{Tr}_{A}\left[\lvert E_{1}\rangle\langle
E_{2}\rvert\right]$ with $E_{1}\neq E_{2}$ to be small. We will adapt these
statements, in their subsystem dependent version (which doesn’t need the
restriction $D_{B}\ll D_{A}$ to few-body subsystems), for our present context
as follows:

$\displaystyle\textrm{Tr}_{B}\left[\left(\rho_{B}(E)-P_{B}(E)\right)^{2}\right]$
$\displaystyle=O(D_{A}^{-1}),$ (70)
$\displaystyle\textrm{Tr}_{B}\left[q_{B}^{2}(E_{1},E_{2})\right]$
$\displaystyle=O(D_{A}^{-1}),$ (71)

where we use the notation $x^{2}=xx^{\dagger}$ for an operator $x$ for
simplicity. Eqs. (70) and (71) should be considered leading order constraints
on the order of magnitude of these quantities when $D_{A},D_{B}\gg 1$, as
noted in the main text. They are also slightly different in some minor
technical details from the main statements of Ref. [29], which we will refer
to as the ‘original conjecture’ in this appendix, and we will now comment on
these differences.

We replace the density of states $\Omega(E)$ with its $O(D)$ scaling behavior
in all subsequent discussions though the original conjecture is stated in
terms of $\Omega(E)$. This is justified by assuming an $O(1)$ spectral width
for the $D$ energy levels and that $\Omega(E)$ is of a comparable order of
magnitude throughout the spectrum (consistent with e.g. $t_{H}=O(D)$ in fully
chaotic systems). As the PSFF involves averages over the entire spectrum, it
is only this scaling behavior that is of interest to us rather than
$\Omega(E)$-dependent variations in smaller regions of the spectrum.

The smallness of $(\rho_{B}-P_{B})$ and $q_{B}$ are enforced above by
requiring the trace of their squares $\textrm{Tr}_{B}\left[x^{2}\right]$
(which we will generally call purity) to be $O(D_{A}^{-1})$. However, the
original conjecture is stated in terms of the trace norm
$(1/2)\textrm{Tr}_{B}\left[(x^{2})^{1/2}\right]$ restricted to be
$O(\sqrt{D_{B}/D_{A}})$. As Ref. [29] notes, on account of the inequality
$\left\\{\textrm{Tr}_{B}\left[(x^{2})^{1/2}\right]\right\\}^{2}\leq
D_{B}\textrm{Tr}_{B}\left[x^{2}\right]$ the constraints in terms of purity
would imply the original conjecture but are also slightly stronger, and it is
in fact these stronger constraints that they verify numerically. We use the
stronger statement because it is more convenient for our purposes, and also
because there appears to be no compelling theoretical reason to rule out such
stronger statements in general. For instance, Ref. [29] motivates the diagonal
statement of the original conjecture in terms of the trace norm based on
analogous canonical typicality [129, 130] constraints for the thermalization
of Haar-random superpositions of energy eigenstates derived in Refs. [131,
132]; but in the process of the derivation in the latter, constraints in terms
of purity similar to Eq. (70) are also seen to hold. We also note that the