Probing many-body quantum chaos with quantum simulators

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.

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-timeordered 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 [1,[20][21][22][23][24] 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 * These authors contributed equally.
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 [11, 35-37, 39-46, 48-57] and (out-of-time-ordered) correlation functions of Heisenberg operators [18,38], 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.

A. Spectral form factor
The SFF in a many-body quantum system with timeindependent 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 Here, we normalize K(t) such that K(0) = D −2 Tr [1] 2 = 1, with D the Hilbert space dimension and have defined the unitary time-evolution operator T (t) ≡ 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τ ) = V n (n ∈ N) and V the Floquet time evolution operator for a single period τ [60].
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 ∈ N of driving period τ is given by T (t = nτ ) = V n 3 with, Here, the Hamiltonians H (x,y,z) contain nearestneighbor interactions with strength J = 3τ −1 and dis- Illustration of the characteristic properties of the SFF and PSFF using the chaotic spin-1/2 Floquet model V3.
(a) We display the SFF K(t) for the Floquet model V3 with N = 6 qubits as a function of time t. We observe characteristic features such as the ramp between t ∼ τ to t = tH = 2 N τ and a plateau for t > tH . (b) For the PSFF KA(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 NA 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 C) by simulating M = 2 × 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.
and σ a [a ∈ (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τ ) = V n 3 , each with local disorder potentials h 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/τ ≈ D = 2 N . This ramp-plateau structure of the SFF is a characteristic feature of quantum chaotic systems [1,77,78], 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 = τ n) = V n , with V from the circular unitary ensemble (CUE), yields Here, the slope of the ramp and the onset of the plateau is determined by the Heisenberg (or plateau) time t H which is connected to the mean inverse spacing of adjacent (quasi-) energies. It typically scales with the Hilbert space dimension t H /τ ∼ D -for V from CUE, t H /τ = D [1,64]. Thus, the SFF is expected to drop with increasing Hilbert space dimension D = 2 N , as D −2 at times 1 t/τ D and as D −1 at times t/τ 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 Th [65]. For the model V 3 we note that t Th ≈ 5τ (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 [30,[32][33][34]. This corresponds to a flat SFF without a ramp which is, after an initial transient regime, constant in time [1], K(t 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 [58,69,70].

B. 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 ⊆ S of the total system S with complement B (A∪B = S) and Hilbert space dimensions D A and D B respectively (D = D A D B ), we define the PSFF as where ρ B (E i ) = Tr A [|E i E i |] denotes the reduced density matrix obtained after partial trace of the eigenstate |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) = Tr B Tr A [1] 2 /(DD A ) = 1. Hence, the SFF and PSFF coincide when A = S, i.e. K A=S (t) = K(t). We emphasize that for A ⊂ S, the PSFF K A (t) contains nontrivial contributions from the eigenstates |E i : We obtain terms of the form Tr(ρ B (E i ) 2 ) and Tr(ρ B (E i )ρ B (E j )) (i = 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 symmetryprotected 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 N/2 = 3, which holds more generally for D A 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 = τ n) = V n , with V from the CUE, and sufficiently large subsystems A, B, (D A , D B 1), 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 [33,[80][81][82][83][84]. 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 Tr B (ρ B (E) 2 ) and of the overlaps of reduced eigenstates Tr B (ρ B (E i )ρ 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 B 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 B. For chaotic models, the shift can be identified as the time independent constant during the linear ramp phase, and for D A D it is approximated by 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 2 A ). 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.

C. Randomized measurements of spectral form factors
Initially, randomized measurements have been proposed and experimentally implemented to characterize many-body quantum states [11, 35-37, 39-46, 48-57] and (out-of-time-ordered) correlation functions of Heisenberg operators [18,38]. 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.
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 ρ0 = |0 0| ⊗N . Before and after the time evolution T (t), we apply random local rotations U = i ui and U † , respectively, where local unitaries ui are sampled from a unitary 2−design. Here, T (t) can be generated as Hamiltonian evolution, T (t) = exp(−iHt), or Floquet dynamics, T (t = nτ ) = V n , n ∈ N, where V is Floquet evolution operator for time period τ . In the last step, a single-shot measurement is performed in the z−basis to collect a bitstring of the form s = (s1, s2, ..., sN ) with si ∈ {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.

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 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 Tr [T (t)] of the time evolution operator and its adjoint uniformly. For instance, in the trivial case T (t = 0) = 1, we obtain that the 'time-evolved' state always matches to the initial random state corresponding to D −1 Tr [T (0)] = 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.
In detail, the experimental recipe reads as follows: (i) We begin with a product state ρ 0 = |0 0| with |0 ≡ |0 ⊗N . (ii) On this initial state, we apply local random unitaries U = N i=1 u i where u i are the local unitaries independently sampled from a unitary 2-design [85,86] on the local Hilbert space 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 C 2 include the (discrete) singlequbit Clifford group as well as uniformly distributed unitary 2 × 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 Lastly, we perform a single-shot measurement in the computational basis with outcome bitstring s = (s 1 , . . . , s N ) with s i ∈ {0, 1} for i = 1, . . . , 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 s (r) with r = 1, . . . M is collected.

Estimation formulas and illustrations
The statistics of the measured bitstrings s (r) , r = 1, . . . 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 where |s| ≡ i s i . As we show in Sec. IV, K(t) yields an (unbiased) estimate of the SFF for a finite number M of experimental runs and converges to K(t) when M → ∞.
Remarkably, from the same measurement data s (r) , we have also access to the PSFF K A (t) for arbitrary subsystems A ⊆ S via post-processing. To this end, we simply project the measured bitrings on the subsystem A of interest, i.e., define s A = (s i ) i∈A , and use which gives an (unbiased) estimate for K A (t) for finite M and converges to K A (t) when M → ∞ (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 × 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.

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 ∼ 10 N A ≈ 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 ∼ 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 = 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 summa-rize 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 A, 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 B, 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 C, and summarize our main conclusions for all cases in Sec. II D.
Common to all these cases is the fact that the timeindependent part of the PSFF in Eq. (4) is given by the plateau value, which depends only on the eigenstate purities (assuming no degeneracies) i.e. K A (t → ∞) = P B /D A , where is the (spectrum-and ensemble-)averaged purity of the reduced energy eigenstates. For later reference, we separate out this time-independent plateau value, and note that the time-dependent second term only involves overlaps of distinct energy levels.

A. Random matrix ensembles
To understand the essential features of the PSFF we first analyze it in RMT, allowing for an exact determina-tion of the PSFF. We choose Hamiltonians H (Floquet operators V ) from the canonical Wigner-Dyson random matrix ensembles [1,20,21,24], yielding time evolution operators T (t) = exp(−iHt) [T (t = τ 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 where and 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 Q B and a constant subsystem dependent shift (P B − 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 P B and Q B explicitly using Wigner-Dyson RMT for the eigenstates of The analogous expressions for orthogonal Wigner-Dyson ensembles can be found in App. B. In both symmetry classes at D A , D B 1, we find that, P B − Q B ≈ 1/D A and Q B ≈ 1/D B . Thus, in this limit, the PSFF has a constant shift of 1/D 2 A 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) ≈ K(t) + 1/D 2 A [see also Eq. (5)].

B. 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, Here, the first term is a constant corresponding to a maximally mixed reduced density matrix; ∆ρ B (E) is traceless and a smooth function of E, while δρ 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 ρ O(D), we take their respective Fourier transforms to satisfy (∆ρ B (t)) jk = 0 for |t| > t ρ , and (δρ B (t)) jk = 0 for |t| < t ρ (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 −1 ρ containing several levels, while the fluctuating part varies only over scales much smaller than t −1 ρ . We will further assume that δρ 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 Tr ρ , fluctuating around an average of zero (we do not require this behavior to persist over larger energy scales |E i − E j | t −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 ∼ t −1 ρ , effectively randomizing the fluctuations δρ B (E) faster than this scale, while ∆ρ 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.
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. (4) 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 δρ B (E) for terms involving overlaps of the smooth and fluctuating part for t t ρ . We can write this as, where ∆K A (t) involves only overlaps of the form Tr B [∆ρ B (E i )∆ρ B (E j )] and similarly, δK A (t) involves only those of the form Tr B [δρ B (E i )δρ B (E j )]. On decomposing δK A (t) in a manner analogous to Eq. (10), it follows that its time dependent part for t t ρ (which sees contributions only from variations of the overlaps of fluctuating parts within energy windows smaller than ∼ t −1 ρ ) vanishes on ensemble averaging, an important consequence of the randomization of δρ B (E). This leaves only a constant contribution from the purity of the fluctuating part, δK (here we use 'purity' to generally mean Tr x 2 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 ∆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 t ρ . We note that ∆K A (t) is modulated only by a smooth function of two energy variables varying over scales larger than t −1 ρ . For t t ρ , 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 t Th , t ρ , Here, β = 1, 2 respectively for the orthogonal and unitary classes, while γ = i Ω −1 (E i ) is the range of energies in the spectrum with Ω(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 ∼ t H [see, for instance, the exact form of the GOE SFF in Eq. (A2)]; we focus instead on the t t H regime where the ramp appears linear for all values of β and profiles of Ω(E), and the t t H regime with a constant plateau. However, both expressions are exact throughout the range of times when β = 2 with constant density of states Ω(E) = t H /(2π). We have also defined two ensemble-averaged quantities corresponding to slightly different spectrum averages of the purity of the smooth part, , 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 P B = D −1 B + ∆P B + δP B , giving the expected plateau value of P B /D A in Eq. (15). There are also two competing time scales for the onset of the ramp, t Th and t ρ -the former entirely determines the behavior of K(t) but the latter appears in ∆K A (t) and δK A (t).
For direct comparison with numerics, it is useful to define the ensemble averaged overlap of adjacent states, Using Eq. (13), we note that, which follow from the assumption of uncorrelated δρ B (E) in the ensemble, and taking ∆ρ B (E i ) ∆ρ B (E i+1 ). We note that this definition of Q B is equivalent to that in Sec. II A for random matrix ensembles, where the ensemble averaged overlaps between distinct states are independent of their energies. Sec. III will directly use P B and Q B , with the implicit assumption that ∆P B is of similar order of magnitude to ∆P B (due to Ω(E) being of a similar order of magnitude throughout the spectrum) and is therefore similarly well represented by Q B .

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 subsystemdependent parameters P B , δP B and ∆P B (∼ ∆P B ) that respectively influence the plateau value, the magnitude of the shift and the magnitude i.e., slope of the ramp. The purity 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 δP B and ∆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], which cannot be less than D −1 B as well as D −1 A . 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, Non-thermal eigenstates are characterized by much larger fluctuations, δP B O(D −1 A ), with ∆P B being correspondingly smaller so as to satisfy the constraint P B = D −1 B + ∆P B + δ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 A); 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.

C. 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 P B through the PSFF. Fully localized states are essentially nearly pure states with P B ∼ O(1) (more precisely, following an area law of entanglement [31]), and additionally have 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 nontrivial correlations between nearby energy eigenvalues emerge in the spectrum, leading to a time dependence of K(t) for t > t ρ , δ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.

D. 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 ∼ 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 timedependent 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 sizesanalyzing 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 P B and the overlap 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. DA is plotted in black circles and matches with the averaged purity PB plotted with red crosses. The average overlap QB and the difference PB − QB 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 V3. 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.

Example 1: Floquet system
The Floquet time evolution operator V 3 has the same quasi-energy eigenvalue statistics as the CUE random matrix ensemble [58,81]. 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 A 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 den-sities of eigenstates P B and 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τ ) the PSFF becomes K A (t → ∞) = P B /D A , see Eq. (10). We plot numerically obtained K A (∞)D A in black circles, and the average purity P B with red crosses, they confirm the analytic expectation. The average overlap Q B and the difference P B − 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) 15) and thus we note that the ramp coefficient is On the other hand, the purity of the fluctuating part (of the form of Tr[δρ 2 B (E)]) comes in the time-independent term added to the SFF in Eq. (15), which is to the leading orders 1/D 2 A , 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.

Example 2: Hamiltonian system
As our second example, we consider a transverse field Ising model in presence of longitudinal local disorders, where h i are drawn uniformly at random from (−1, 1). The coefficient J and the exponent α denote the strength and range of the interactions respectively. The disorder strength W is known to specify the nature of the dynamics; W ∼ J depicts chaotic regime and W 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 α and find that the chaotic and localized phases exist for short (α > 1) as well as for long (α < 1) range interactions. In this work, we choose α = 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 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 C. 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 (σ 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 Th [1,24,27,105,106], thus we have also put the results for GOE class in gray in Fig. 4.  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 P B using crossed markers. We have also plotted here the plateau values K A (∞)D A (in black circles) for both chaotic and MBL phases which agree with their respective purities following K A (t → ∞) = 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 Q B and P B − 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 B. Utilizing Eq. (16) we can comment on the orders of ∆P B and δP B (see App. F 2 for more details on the numerical extraction of these orders).
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 δ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 The localized phase is not expected to satisfy ETH, and as discussed in the Sec. II C, we expect such large shift in the PSFF in MBL systems. Due to larger δ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

IV. PROOF OF THE PROTOCOL
In Sec. I C, 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.

A. 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 | s| U † T (t)U |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 = i u i , one finds for any operator C defined on the 'two-copy' Hilbert space (20) Here, E U denotes the average over local unitaries of the form U = i u i with u i sampled for each i independently from a unitary 2-design on the local Hilbert space C ⊗2 . Further, the sum extends to all two-copy permutation operators σ = i σ i and τ = i τ i with σ i , τ i = 1 i , S i . Here, the identity 1 i and the swap operator S i act as 1 i |s i ⊗ |s i = |s i ⊗ |s i and S i |s i ⊗ |s i = |s i ⊗ |s i on local basis states |s i and |s i . Finally, the coefficient is determined by the Weingarten function Wg U (2) , with Wg U (2) (1 i ) = 1/3 and Wg U (2) (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 ⊗ρ 0 with ρ 0 = |0 0| and where the sum extends to all bitstrings s = (s 1 , . . . , s N ) with s i ∈ {0, 1}, and |s| ≡ i s i . For this choice, we obtain with S = i S i = s,s |s s| ⊗ |s s |. The Swap operation 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 ⊗ B → A T ⊗ B in the above equation, leading to

B. 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) such that the ensemble (disorder) average K(t) = 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 |Φ + N under the dynamics generated by 1 ⊗ T (t), i.e.
To this end, we use the following identity for any two operators A, B on H which can be proven by inserting the definition of the Bell state |Φ + N = 2 −N/2 s |s ⊗ |s in terms of computational basis states. Eq. (25) follows directly by choosing A = 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 ⊗ ρ 0 averaged over random unitaries U . We find 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 ρ f (t) Here, | s| U † T (t)U |0 | 2 is precisely the Born probability of finding a bitstring s, in the computational basis measurement performed at the end of our measurement sequence when the state ρ f (t) has been prepared [c.f. Sec. I C]. It follows thus that where E QM is the quantum mechanical average and 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 s, (−2) |s| provides an estimation of the SFF, which in expectation over ensemble (disorder) average, over random unitaries and quantum mechanical averaging, yields the SFF 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 K(t) [Eq. (7)]. Using Eq. (30), it follows that K(t) converges to K(t) in the limit M → ∞. For finite M , statistical errors are governed by the variance of 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 C) 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-tonoise 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.

A. 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.

Observing PSFF and SFF
We can bound the statistical errors of the estimator K A (t) [Eq. (8)] by its variance. As shown in App. H, we find that, 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 ⊆ A. The variance of K(t) [Eq. (7)] follows by taking A to be the full system. We obtain an expected relative error E A = σ A /(K A √ M ) of an estimation 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 ∼ σ 2 A /K 2 A .  The expected statistical error 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 ⊆ A. For Hamiltonians (Floquet-) operators from Wigner-Dyson RMT, we can explicitly evaluate σ A (see App. H). As the worst-case estimate, we find that at the point of weakest signal, after a single time step t = τ in Floquet dynamics T (t = nτ ) = V n with V sampled from CUE where K A (t = 1τ ) = 2 −2N A , the expected relative statistical error is given by E A = 10 N A /M . A total number of measurements M ∼ 10 N A / 2 ≈ 2 3.32N A / 2 is thus required to obtain a fixed relative error . 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 / 2 measurements, with r = r(N A ) ≥ 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 E A as a function of the number of experimental runs M in the V 3 model (2) at time t/τ = 5 with total qubits N = 6. The relative error decays as ∼ 1/ √ 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/τ = 5 and t/τ = 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 → N .

Observing the ramp in chaotic models
The relative error 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 , To quantify the experimental effort to resolve c A (t 2 , t 1 ), we introduce its signal-to-noise ratio SNR[c A (t 2 , t 1 )], which, for independent measurements of the PSFF at times t 2 and t 1 , is given by 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 N/2. At the same time, the absolute value of the noise, here (σ A (t 2 ) + σ A (t 1 ))/ √ 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, 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 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 S, we can estimate K A (t) for all subsystems A ⊆ S, via spatial restriction in the postprocessing. 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 SNR[c |A| (t 2 , t 1 )] for these average quantities.
In Fig. 6(d), we display the numerically determined signal-to-noise-ratio 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 SNR[c |A| (t 2 , t 1 )] for subsystem sizes 1 < N A < N . We remark that we do not reach an enhancement N N A 1/2 of the SNR which would result trivially from N N A separate experiments (i.e. N N A · M experimental runs in total, gray line) since the estimations 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 ≈ 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. First, we consider an imperfect implementation of our measurement protocols, with errors arising from an erro-neous 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 † i exp(−iη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 η = 0 corresponds to the ideal case, we display in Fig. 7(a) the average relative error η = 1 − (K A (t)) η / K A (t) of the estimated (K A (t)) η as a function of the error strength η, obtained numerically from simulating many experimental runs. We find η increases approximately as η 2 , indicating a decrease of the estimated (K A (t)) η .

B. Experimental imperfections
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 [11,[35][36][37], 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 τ with strength p, i.e. the final state ρ f (t) at time t = τ n, defined in Sec. I C, is altered to ρ dec (t) = α n ρ f (t) + (1 − α n ) 1/D with α n = (1 − p) n .
Thus, we obtain via our measurement protocol, With increasing time t = τ n, decoherence leads thus to a smaller measured value (K A ) 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, which gives [37,114], Thus, from a measurement of the purity P n at all times, we can find α n and rescale the erroneous PSFF (35) to obtain, 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 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 A). In addition, we have defined partial spectral form factors (PSFFs), K A (t) in Eq. (4), which contain both the statistics of energy eigenvalues and eigenstates (see Sec. I B). PSFFs are natural restrictions of the SFF to subsystems A ⊆ S of the full system S, such that for A = S, PSFF and SFF coincide K A=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-rampplateau 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 C, 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 manybody 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-93, 121, 122], as well as the investigation of noise-induced quantum many-body phenomena such as entanglement phase transitions [46,[123][124][125]. 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]. 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 [1,20,21,24]. 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 timereversal 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 = τ n) = V n , n ∈ N, where V is the unitary Floquet operator corresponding to a time period τ , 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 [33,[80][81][82][83][84].
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: GOE model: 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 COE model: Here, t H = Dτ with τ 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 S be a quantum system with Hilbert space H of dimension D, and A ⊆ S a subsystem with dimension D A . Its complement is denoted with B with dimension D B . As discussed in App. A, we consider • Floquet dynamics with T (t = τ n) = V n for n ∈ N with V sampled from the CUE and COE, respectively.
We can rewrite 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. (4) explicitly (see next subsection). With the identification K(t) = D −2 |Tr [D(t)] | 2 , we find where for H ∈ GUE , V ∈ CUE, and for H ∈ GOE , V ∈ COE, 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. The Weingarten calculus for the unitary group and for orthogonal group can be formulated in terms of pair partitions, defined as follows. The following fact is shown in Ref. [128].  In our case, we are only interested in the case n = 2. As shown in Ref. [128], when m, n ∈ M O (4) and D ≥ 2, Furthermore, it holds for m, n ∈ M U (4) and D ≥ 2 for m = n . (B12) Using Fact 2 and these expressions, we can perform the average over eigenvector elements in Eq. (B6) 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 B 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, ρ B (E) = D −1 B 1 + ∆ρ B (E) + δρ 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 ρ O(D), The remaining energy dependent part of ρ 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, (C2) Up to this point, such a decomposition is always possible. We will additionally take t ρ to be set by the scale of randomization in the ensemble discussed in Sec. II B, 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 −1 ρ . Similarly, we will not require randomization of the correlators of δρ B (E) between energies further apart than ∼ t −1 ρ , for which the correlator may have to be nonvanishing to maintain zero Fourier component of the fluctuating part at t ≤ t ρ .
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 transformsf (t) andg(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, Now, it is convenient to define an ensemble-averaged unequal time SFF K(t 1 , t 2 ) = D −2 j,k e iEj t1−iE k t2 , which reduces to K(t) at equal times t 1 = t 2 = t. The sum of phases D −2 j,k e iEj (t+t l )−iE k (t+tr) in Eq. (C3) 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 ), 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. (C4) 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 ∼ t −1 Th ) in this class takes the universal form [1,24,97], where Ω(E) is the smoothened (continuous and ensembleaveraged) density of states, whose Fourier transform sat-isfiesΩ(t t Th ) ≈ 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. (C5). 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 t Th , where β = 2 for the unitary Wigner-Dyson class, and we have introduced the shorthand symbols T 12 = (t 1 +t 2 )/2, τ 12 = t 2 − t 1 .Θ Ω (t) is the Fourier transform of the unit step function Θ(Ω(E)), the latter being 1 where Ω(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πΩ(E) with the Heisenberg time t H , assuming that Ω(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 ∼ t H . But for t t H , virtually the same results hold with β = 1 for the orthogonal class and β = 4 for the symplectic class [97] (of course, the plateau behavior for t t H is generally independent of such specifics). Analogously, we expect similar replacements (the appropriate value of β, and focusing on the T 12 t H and T 12 t H regimes) to work for the unequal time SFF in Eq. (C6) as well. With this expectation, we write for t 1 , t 2 t Th in any Wigner-Dyson symmetry class. Using the decomposition of ρ B (E) with these definitions then gives several terms for K A (t) of the form of Eq. (C3), where f and g independently go over each of D −1 B , ∆ρ B and δρ B , with an additional trace of the product over the B subspace. Now, we will argue that all cross terms with f = g may be taken to vanish. When , which is zero when g = ∆ρ B , δρ B , which are both traceless. When say, f is ∆ρ B and g is δρ B , the cross term vanishes due to the assumption that ensemble averaging randomizes δρ B .
Dropping the cross terms for the above reasons gives the form of Eq. (14) in the main text, K A (t) = K(t) + ∆K A (t) + δK A (t), where K(t) is the full SFF, and In the main text, it is argued that δK A (t t ρ ) amounts to a constant shift after ensemble averaging due to the randomization of δρ B (E). Here, we will complete the evaluation of ∆K A (t) using the prototype Eq. (C4) with f = g = (∆ρ B ) ab and the expression in Eq. (C7) with t 1 = t + t l , t 2 = t + t r . As the definition of ∆ρ B sets t l , t r < t ρ , we have |T 12 | = |t| + sgn(t)(t l + t r )/2 at large times (i.e. t t Th , t ρ ). For t t H in this regime, this gives, The Hermiticity of ∆ρ B implies that (∆ρ B (−t)) ab = (∆ρ * B (t)) ba . Consequently, making the integration variable transformation t l → −t r , t r → −t l in Eq. (C10), we see that inside the parentheses in the second line the |t| term is unaltered but the sgn(t) term transforms to its negative, while all factors outside the parentheses remain unaltered. It follows that the contribution from the sgn(t) term actually evaluates to zero, leaving only a linear ramp term from |t|. For t t H , we directly obtain only a plateau contribution. Now, it is straightforward to Fourier transform back to the energy variable E, 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: i ↔ dE Ω(E) and i Ω −1 (E i ) ↔ dE Θ(Ω(E)), which become equalities on ensemble averaging. Then we get the expression, Together with the expression for the full SFF [t 1 = t 2 in Eq. (C7)] 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, P B (purity of reduced density matrices), δP B (fluctuating part) and ∆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 a, and argue for its extension to eigenstates of arbitrary delocalization in App. D 2 b.

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 nonthermal 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 ≤ D that it has a significant probability of being found in.
It is useful to introduce an effective dimension D eff A ≤ D A , 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 eff A = D A if the eigenstates appear completely delocalized over subsystem A, and more generally D eff A is typically larger for larger D A (up to ). For instance, D eff A 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 eff A = 1 for D A = 1, to D eff A = 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 1. In particular, we will take D D A , D B 1.

Assuming that D eff
A is typical for A throughout the spectrum, the purity in subsystem B satisfies, . The first two terms are due to the eigenstate being delocalized in subsystem B with effective dimension ( /D eff A ), with the second term containing larger scale variations of its components. We will call this, the 'macroscopic' contribution, which grows with D eff A . The last term is due to the randomness of the eigenstate components i.e. the 'microscopic' contribution, which decays with D eff A (and is also typically bounded from below by (1/D eff A )). Being a linear combination of the macroscopic and microscopic contribution, the purity shows an initial decay with D eff A for small values of the latter, and eventually a growth for larger values of D eff . Both D eff A = 1, correspond to pure states with P B = 1.
The parameters δP B , ∆P B satisfy the following orderof-magnitude inequalities, 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 P B = D −1 B + ∆P B + δ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 −1 B ); 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. (D2) and (D3) are satisfied as equalities for thermal eigenstates. In particular, δP B can only decay with increasing D eff A -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. (D2) and (D3) much further from equality. In this case, the macroscopic contribution to the fluctuating part may even show up as a growth of δP B with D eff A 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 δP B = ∆P B = 0 for D A = D.
We conclude that P B is a measure of delocalization of eigenstates, while δP B and ∆P B are probes of thermalization. Setting = D gives the results discussed in the main text for chaotic systems with fully delocalized eigenstates (Sec. II B 2). For fully localized sys- automatically implying a lack of thermalization (Sec. II C). 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.

Subsystem ETH constraints a. 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 |E , there are two statements of the hypothesis: the diagonal statement stating that the reduced density matrix ρ B (E) = Tr A [|E E|] 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 ) = Tr A [|E 1 E 2 |] with E 1 = E 2 to be small. We will adapt these statements, in their subsystem dependent version (which doesn't need the restriction D B D A to few-body subsystems), for our present context as follows: where we use the notation x 2 = xx † for an operator x for simplicity. Eqs. (D4) and (D5) should be considered leading order constraints on the order of magnitude of these quantities when D A , D B 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 Ω(E) with its O(D) scaling behavior in all subsequent discussions though the original conjecture is stated in terms of Ω(E). This is justified by assuming an O(1) spectral width for the D energy levels and that Ω(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 Ω(E)-dependent variations in smaller regions of the spectrum.
The smallness of (ρ B − P B ) and q B are enforced above by requiring the trace of their squares Tr B x 2 (which we will generally call purity) to be O(D −1 A ). However, the original conjecture is stated in terms of the trace norm (1/2)Tr B (x 2 ) 1/2 restricted to be O( D B /D A ). As Ref. [29] notes, on account of the inequality Tr B (x 2 ) 1/2 2 ≤ D B Tr B x 2 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 Haarrandom 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. (D4) are also seen to hold. We also note that the purity constraints remain < O(1) for D B > D A , whereas the corresponding constraints on the trace norm (which cannot be greater than 1 for differences of density matrices [133]) are > O(1) and therefore meaningless in this regime. The original conjecture had to restrict the subsystemdependent form to D B < D A (in our notation) for this reason. However, in Sec. III of the main text, we find numerical support for the validity of Eqs. (D4) and (D5) even for D B > D A .
Finally, we note that the smooth reduced density matrix P B (E) is not precisely characterized in Ref. [29] but it is also unnecessary to be too precise in specifying it as Eq. (D4) is only an order-of-magnitude constraint.
Here, in analogy with Eq. (C1), we will define P B (E) to be that part of ρ B (E) that varies slower than some rate t s , effectively amounting to a weighted average of ρ B (E) over energy windows of size ∼ t −1 s . We will assume Eq. (D4) is satisfied for any choice of t s larger than some minimum magnitude ∼ t ETH O(D) (intuitively, because the more the smooth part is allowed to fluctuate, the more closely it can approximate ρ B (E)). Then, if our ensemble is such that t ρ t ETH , we can choose t s = t ρ . This allows the identification P B (E) = D −1 B + ∆ρ B (E) in the decomposition ρ B (E) = D −1 B + ∆ρ B (E) + δρ B (E) of Eq. (13). Eq. (D4) then gives, The constraint δP B = O(D −1 A ) then follows directly from here.
To similarly obtain a condition from Eq. (D5) that applies directly to the PSFF, we note that this equation can be rewritten in terms of reduced density matrices of the complementary subsystem A as On taking the ensemble average, and using the expansion of ρ A (E) in terms of its smooth and fluctuating parts, the contribution from the fluctuating part δρ A (E) to the left hand side vanishes due to the randomization assumption in Sec. II B. We are then left with

neighboring levels) so that the second term is approximately Tr
. From this, we get the smooth purity constraint ∆P A = O(D −1 A ) on taking the appropriate spectrum averages. In the context of ∆P B (and ∆P B ) in the main text, these purities are evaluated in subsystem B rather than A, and the corresponding constraints are therefore consequences of off-diagonal subsystem ETH, Eq. (D5), applied to subsystem A instead of B.

b. Extension to partially delocalized eigenstates
We begin with a complementary approach to that of the previous subsection, to argue that the purity based expressions of subsystem ETH should generally hold for chaotic systems with fully delocalized eigenstates. Consider requiring each matrix element of ρ B (E) to differ from the corresponding matrix element of P B (E) only by a small amount O( √ D A /D), as a stronger diagonal statement that implies Eq. (D4) (a weaker, D A -independent version of such a statement is also considered in Ref. [29]). To justify this constraint, we consider the following situation. Let |E 1 and |E 2 be two 'typical' nearby eigenstates that are completely delocalized over the D basis vectors (in some 'physical' product basis of subsystems A and B) with random (real or complex) phases. Their density operators ρ(E 1 ) = |E 1 E 1 | and ρ(E 2 ) = |E 2 E 2 | have matrix elements of the schematic form The difference ρ(E 1 ) − ρ(E 2 ), after a partial trace over A, can be taken to represent the fluctuations of ρ B (E) around P B (E). Given our above assumptions on the eigenstates, the matrix elements of ρ(E 1 ) − ρ(E 2 ) are typically ∼ O(D −1 ) in magnitude with random signs or phases (i.e. with zero 2-point correlation, which crucially requires even large-scale non-uniformities in the magnitudes to agree up to random fluctuations). The sum of D A such matrix elements in the partial trace over subsystem A then has magnitude O( √ D A /D), justifying the above constraint. Similarly, the operator q(E 1 , E 2 ) = |E 1 E 2 | for such eigenstates has O(D −1 ) matrix elements with random phases, giving O( √ D A /D) matrix elements after the partial trace and therefore the offdiagonal statement Eq. (D5). Such a picture of random energy projector matrix elements of comparable magnitudes is reminiscent of Berry's conjecture for chaotic wavefunctions [134] (as well as other related statements e.g. Refs. [2,[94][95][96]), which has been interpreted as the origin of eigenstate thermalization in chaotic systems [2,3].
Using an analogous argument for eigenstates that are not necessarily delocalized over all D basis vectors, we can clearly highlight the difference between delocalization and thermalization, and the distinct information contained in the overall purities as opposed to the smooth and fluctuating parts of the reduced density matrices. For this purpose, consider an eigenstate |E 1 that is randomly (but not necessarily uniformly) distributed only over a set of ∼ ≤ D 'physical' basis vectors, with negligible support outside this set. Its density matrix ρ(E 1 ) then has an × block (after suitably permuting rows and columns) of non-vanishing elements each of typical magnitude O( −1 ), and all elements outside this block may be taken to vanish. As always, all the diagonal elements are strictly non-negative and add to 1, while the independent off-diagonal elements could have arbitrary signs or phases (which are typically random). Thus, we have the schematic form, where Θ(x) = 1 if x is true and 0 otherwise, and p k denotes the index corresponding to k after a permutation p of rows/columns. The behavior of ρ(E 1 ) under a partial trace depends on the choice of the subsystem A. We will choose subsystems which can be traced out by factorizing the chosen basis (which means the basis states are pure states within the subsystem). This identifies a class of subsystems which are sensitive to the specific extent of delocalization of eigenstates; in a more general basis in the Hilbert space, the eigenstates may appear delocalized by an arbitrary extent, including fully localized in the energy eigenbasis and generically fully delocalized ( = D) in a Haar random basis according to canonical typicality [130,131]. An equivalent, more physically motivated viewpoint is that the extent of delocalization of eigenstates should be determined by their minimum such delocalization in bases comprised of nearly pure states (e.g. a product basis) in most physically accessible subsystems -so that a small subset of eigenstates may be treated as if they each have independent random components (neglecting the global constraint of orthonormality) under a (sufficiently small) partial trace.
For convenience, we first consider the case where the eigenstate looks fully delocalized in subsystem A within its support on the physical basis -in other words, the partial trace over A does not mix the zero and nonzero elements of ρ(E 1 ). In this appendix, we will call such a subsystem A an unbiased subsystem (from the point of view of the eigenstate of interest). Then ρ B (E 1 ) has an ∼ ( /D A ) × ( /D A ) non-vanishing block with nonnegative diagonal elements of magnitude O(D A / ), and off-diagonal elements of typical magnitude O( √ D A / ) in the case of an eigenstate with random phases (as long as the partial trace combines several basis vectors where the eigenstate has comparable magnitudes). Now, we can evaluate the purity Tr B ρ 2 B (E 1 ) , which sees a net contribution of O(D A / ) from the diagonal elements and O(D −1 A ) from the off-diagonal elements. Additionally, normalization requires that the diagonal elements must add up to 1, therefore the sum of their squares is greater than or equal to ∼ /D A -the inverse of the number of diagonal elements. Their contribution to the purity can then be written in a more descriptive form as Thus, we can extract information about the extent of delocalization, , by looking at the subsystem size dependence of the purity. We note that the purity can also be written as Tr A ρ 2 A (E 1 ) from the viewpoint of subsystem A giving an additional lower bound of D −1 A , which is mostly contained in the O(D −1 A ) term for D A 1 (as the diagonal contribution to purity from ρ A (E 1 ) is primarily due to contributions from the off-diagonal elements of ρ B (E 1 )).
A nearby eigenstate |E 2 that is also distributed only across basis vectors (but not necessarily the same ones or in the same way as |E 1 ) again shows a subsystem purity of the form of Eq. (D11). The two eigenstates thermalize if their reduced density matrices do not differ significantly, in small enough subsystems that trace out a lot of the independent eigenstate components. This would be the case if these two eigenstates are distributed across roughly the same basis vectors in a largely similar manner (up to random fluctuations). From this point of view, subsystem ETH is a qualitative identification of the thermalization of a set of otherwise random-looking eigenstates with the extent of their 'overlap' within subsystems, rather than merely with entanglement as represented by their individual purities (the latter being the canonical typicality approach that is only sufficient for fully, uniformly delocalized random eigenstates as in App. B).
We now consider two illustrative extreme cases of fully overlapping (thermal) and fully non-overlapping (non-thermal) eigenstates. In both cases, we will be interested in Tr B (ρ B (E 1 ) − ρ B (E 2 )) 2 as a representative of the size of the fluctuating part [ρ B (E) − P B (E)] of reduced energy eigenstates in subsystem B, as well as the (real-valued) overlap Tr B [ρ B (E 1 )ρ B (E 2 )] which is equal to the norm of off-diagonal operators Tr A [q B (E 1 , E 2 )q B (E 2 , E 1 )] in subsystem A. These are complementary quantities, being related to the subsystem purities of the individual eigenstates via This relation quantifies the identification of thermalization with overlap.
• Thermal eigenstates: If |E 1 and |E 2 are distributed in a similar manner across the same basis vectors, then again has (ρ(E 1 ) − ρ(E 2 )) an × block structure, with random O( −1 ) off-diagonal elements within the block. However, the diagonal elements, being differences of random O( −1 ) non-negative numbers, also have at most O( −1 ) magnitudes with random signs (if large scale nonuniformities match), and largely cancel each other out in a partial trace. After the partial trace, all matrix elements of (ρ in magnitude, and we have which is the analog of off-diagonal ETH, Eq. (D5), for subsystem A.
• Non-thermal eigenstates: In the non-thermal case, |E 1 and |E 2 are distributed in completely different ways, and the diagonal elements of (ρ(E 1 ) − ρ(E 2 )) do not have completely random signs among elements with comparable magnitudes. Consequently, there is no longer a significant cancellation of the diagonal elements in a partial trace for a general choice of A. The fluctuating part Tr B (ρ B (E 1 ) − ρ B (E 2 )) 2 is then typically much larger than O(D −1 A ) with some O(D A / ) contribution, and the overlap is correspondingly smaller. In the extreme case of the two eigenstates being distributed across completely different basis vectors, (ρ(E 1 ) − ρ(E 2 )) has two different × blocks, and the reduced difference in subsystem B also has the structure of two independent blocks. We then obtain behavior analogous to the subsystem purities, while the overlap for this case vanishes entirely, We note that these trends hold only for D A < , due to the assumption on subsystem A. The reduced energy eigenstates in subsystem B are pure basis states when D A = , and behave accordingly on a further partial trace.
The fluctuations in reduced energy eigenstates and their overlaps therefore contain information about eigenstate thermalization that is not visible to the purity alone, which is merely an indicator of eigenstate delocalization. We also see that, at least for 'typical' eigenstates, diagonal subsystem ETH should be understood (in a coarse, order of magnitude sense) as a lower bound relation, while off-diagonal subsystem ETH is a complementary upper bound relation, related through Eq. (D12) to each other and the subsystem purities. In place of Eqs. (D4), and (D5), we can therefore write the more general relations for partially delocalized eigenstates, when A is an unbiased subsystem, with D A ≤ . Both bounds are saturated by thermal eigenstates. For greater completeness of the present discussion, we should account for a more typical choice of subsystem A -one that would mix the zero and non-zero elements of these eigenstate reduced density matrices on performing the partial trace over A. We will consider such a typical subsystem to have an effective dimension D eff A ≤ D A , corresponding to the typical number of non-zero density matrix elements added together in the partial trace. This can be thought of as a generalization of the notion of effective dimension, discussed for the case of infinite dimensional Hilbert spaces in Ref. [29]. We ignore the more complicated case where the number of matrix elements added together is not approximately uniform for all nonzero matrix elements (and therefore, no effective subsystem dimension exists), with the belief that it would not significantly alter our qualitative conclusions. When the effective dimension does exist, all the above conclusions hold for any system but with D A replaced by the smaller quantity D eff A . As an aside, Eqs. (D17) and (D18) continue to hold even without this replacement, but are then not necessarily saturated by thermal eigenstates unless A is an unbiased subsystem.
As a simple example, if subsystem A is unbiased with respect to a set of eigenstates of interest, then its complementary subsystem B has effective dimension D eff B = /D A (note that B is not unbiased). Using this, we can finally write off-diagonal ETH for subsystem B and diagonal ETH for subsystem A as follows, More generally, expressing Eqs. (D17), (D18) in terms of D eff A gives the constraints discussed in App. D 1.

Appendix E: Time-reversal symmetric Floquet thermalization
In this section, we consider another Floquet model of periodically kicked spin-1/2 system. We consider one period of duration τ to be, At multiples t = nτ (n ∈ N) the time evolution of this model is governed by the Floquet time evolution operator where the local disorder potentials h (y,z) i are uniformly and independently sampled from [−J, J]. We fix the driving frequency to τ −1 = J/2. With these parameters, the time evolution operator V n 2 is known to have COE eigenvalue statistics after a few initial kicks [81]. DA is plotted in black circles and matches with the averaged purity PB plotted with red crosses. The average overlap QB and the difference PB − QB are presented in brown and green respectively. We observe a match with the respective quantities in COE plotted in gray, signaling the same averaged eigenvalue and eigenvector statistics in COE and V2.
We present numerically obtained SFF and PSFF for a total system size of N = 6 and subsystem sizes N A = 3, 4 and 5 in Fig. 10(a) and observe a shift in the PSFF in addition to the characteristic chaotic features; the ramp and the plateau. We plot with gray lines the corresponding K A (t) in a COE model where the analytic forms have been exactly calculated [see Eq. (B3)], and observe a good match between V 2 and COE. The match between the statistics of COE and V 2 can further be explored using the second-order moments of the reduced density of eigenstates.
In Fig. 10(b) we present the overlaps P B and Q B as functions of subsystem size N A . We plot numerically obtained K A (∞)D A in black circles, and the average purity P B with red crosses and note a good match between the two. Note that unlike the SFF in the unitary class where the transition to plateau at the Heisenberg time is sharp, the transition to a constant plateau takes a long time in an orthogonal model. This is why we observe slight differences in the numerically calculated plateau value and the purity in Fig. 10(b). The average overlap Q B and the difference P B − Q B are plotted in brown and green circles respectively and match with those in COE. Therefore, similar to RMT models and the model V 3 (in Sec. III), the Floquet dynamics V 2 also has ∆P B = 0. Thus, numerically we confirm that the reduced densities in the Floquet system V 2 also thermalize to infinite temperature and the ramp is governed entirely by maximally mixed part of ρ B (E). From the plots of P B − Q B (in green) we conclude that the constant term added to the SFF is ∼ 1/D 2 A , as in the RMT models. In the main text, we considered the Ising Hamiltonian in Eq. (19) as an example of local many-body models. In this Appendix, we provide some supporting data which were used in the main section. We first begin with analyzing the interesting set of parameters for which we observe chaotic and localized phases in the Hamiltonian model. In Sec. II, we derived the orders for the purity and overlaps of the reduced density matrices on the basis of ETH, and presented them numerically in Sec. III. Here, we provide some additional information on the numerics used to extract the orders for the Hamiltonian model. In the last subsection, we numerically cross-check the shift in the PSFF data with the shift δP B calculated using Eq. (16), where in the latter we directly use the reduced densities of eigenstates.

Chaotic and MBL regimes in Ising Hamiltonian
We explain our choice of parameters in the Ising Hamiltonian Eq. (19). The Hamiltonian contains ZZ interactions with strength J and the range of interactions is given by α. It has a transverse field with strength J and a longitudinal local random disordered field with strength W . Our interests lie in the parameters such that the Hamiltonian dynamics is either in the chaotic phase or in the localized phase. For this purpose, we analyze the energy level statistics, using the adjacent energy gap ratio. From the sorted energy eigenvalues E 1 < E 2 < · · · < E D we compute the energy gaps ∆E m = E m+1 − E m . Then W/J < l a t e x i t s h a 1 _ b a s e 6 4 = " h 2 U U M V m R Y D H h E h 8 7 K X u t r D e 6 O 4 c = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e 5 G Q Y 9 B L + I p o n l A s o T Z y S Q Z M j u 7 z P Q K Y c k n e P G g i F e / y J t / 4 y T Z g y Y W N B R V 3 X R 3 B b E U B l 3 3 2 8 m t r K 6 t b + Q 3 C 1 v b O 7 t 7 x f 2 D h o k S z X i d R T L S r Y A a L o X i d R Q o e S v W n I a B 5 M 1 g d D P 1 m 0 9 c G x G p R x z H 3 A / p Q I m + Y B S t 9 N A 8 u + s W S 2 7 Z n Y E s E y 8 j J c h Q 6 x a / O r 2 I J S F X y C Q 1 p u 2 5 M f o p 1 S i Y 5 J N C J z E 8 p m x E B 7 x t q a I h N 3 4 6 O 3 V C T q z S I / 1 I 2 1 J I Z u r v i Z S G x o z D w H a G F I d m 0 Z u K / 3 n t B P t X f i p U n C B X b L 6 o n 0 i C E Z n + T X p C c 4 Z y b A l l W t h b C R t S T R n a d A o 2 B G / x 5 W X S q J S 9 8 3 L l / q J U v c 7 i y M M R H M M p e H A J V b i F G t S B w Q C e 4 R X e H O m 8 O O / O x 7 w 1 5 2 Q z h / A H z u c P t 5 2 N b A = = < / l a t e x i t > ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " u A K V O m s h A h r 5 j m W q u r Z Y S R U q Z l g = " > A A A B 7 n i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k V 9 F j 0 4 r G C / Y A 2 l M l 2 0 y 7 d b M L u R i i h P 8 K L B 0 W 8 + n u 8 + W / c t D l o 6 4 O B x 3 s z z M w L E s G 1 c d 1 v Z 2 1 9 Y 3 N r u 7 R T 3 t 3 b P z i s H B 2 3 d Z w q y l o 0 F r H q B q i Z 4 J K 1 D D e C d R P F M A o E 6 w S T u 9 z v P D G l e S w f z T R h f o Q j y U N O 0 V i p 0 0 e R j L E 8 q F T d m j s H W S V e Q a p Q o D m o f P W H M U 0 j J g 0 V q H X P c x P j Z 6 g M p 4 L N y v 1 U s w T p B E e s Z 6 n E i G k / m 5 8 7 I + d W G Z I w V r a k I X P 1 9 0 S G k d b T K L C d E Z q x X v Z y 8 T + v l 5 r w x s + 4 T F L D J F 0 s C l N B T E z y 3 8 m Q K 0 a N m F q C V H F 7 K 6 F j V E i N T S g P w V t + e Z W 0 6 z X v s l Z / u K o 2 b o s 4 S n A K Z 3 A B H l x D A + 6 h C S 2 g M I F n e I U 3 J 3 F e n H f n Y 9 G 6 5 h Q z J / A H z u c P w 0 e P M A = = < / l a t e x i t > FIG. 11. Density plot for the mean adjacent gap ratio rm . The behavior of the mean adjacent gap ratio rm as a function of disorder strength W/J and range of interactions α is presented. We notice the presence of chaotic ( rm ∼ 0.53) and localized ( rm ∼ 0.39) phases for a wide range of α. We have worked with α = 1.2, W = J (chaotic), and W = 10J (MBL).
we find, the adjacent energy gap ratio Integrable systems are characterized by a mean ratio of r m ≈ 0.39 whereas the chaotic systems with timereversal symmetry, obeying GOE Wigner-Dyson energy level statistics have a mean r m ≈ 0.53. We use this mean value of r m to choose the parameters for the chaotic and localized phase in our Hamiltonian model. In a density plot of the mean r m as a function of W/J and α in Fig. 11, we notice that the chaotic and localized phases exist for both the short (α > 1) and the long (α < 1) range of interactions. In this work, to discuss the two phases, we have chosen the parameters to be α = 1.2, W/J = 1 (chaotic) and W/J = 10 (localized).

Orders of magnitude of ∆PB and δPB
The subsystem ETH specifies the orders of magnitude for ∆P B and δP B to be O(1/D B ) and O(1/D A ) respectively for the chaotic models. For the localized models, which are known to not satisfy ETH, we concluded that the shift coefficient, δP B O(1/D A ). We note that these orders for the chaotic phase, expressed in terms of P B and Q B as deviations from the RMT prediction, amount to, and In Fig. 12, we plot these quantities for the Hamiltonian model, Eq. (19), for a total of N = 10 qubits. We find that the chaotic phase W = J (in red) satisfies the ETH results whereas for the localized phase W = 10J (in blue), the shift coefficient P B −Q B 1/D A , as predicted in the Sec. II.  . 12. Validating the orders. We test the validity of Eq. (F2) for the chaotic and MBL phases of the Hamiltonian. As expected, the chaotic phase, in red, satisfies the predicted order and the MBL phase, in blue, violates it (in the right plot for PB − QB).

Comparison of the PSFF shift and δPB
Here, we numerically verify the prediction of Eq. (15) for the constant late time shift of the PSFF in the chaotic phase -namely, that the shift is given by δP B /D A in the ramp region. For this purpose, we subtract the full SFF from the PSFF at some time t 0 in the linear ramp region, satisfying t Th , t ρ t 0 t H , which gives This difference has two contributions -the first term is the additive shift which we are presently interested in, but the second term is due to the differing slopes of the linear ramp, from the excess purity of the smooth part of the reduced density matrix. We will now argue that it is reasonable to take for our purposes. In Eq. (F3), by subsystem ETH, the former is O(D −2 A ) while the latter is ∼ O(D −1 (t 0 /t H )) (taking γ ∼ O(1), consistent with t H ∼ O(D)). The second term is therefore negligible if t 0 /t H O(D B /D A ). This is immediately satisfied for any t 0 in the linear ramp region if D A < D B ; conversely, for a given choice of t 0 , D A can be as large as ∼ Dt H /t 0 while maintaining the validity of Eq. (F4). As t 0 t H in general, we expect Eq. (F4) to be a reasonable approximation for a range of values of D A > √ D as well. A minor additional effect that improves this approximation is that for large D A , the coefficient D B ∆P B of the second term would be small, though still O(1), from On the basis of Eq. (F4) and the relation δP B = P B −Q B from Eq. (16), we compare D A (K A (t 0 )−K(t 0 )) for some suitably chosen t 0 to P B − Q B in Fig. 13 and observe good agreement, especially for smaller D A as expected. We note that this agreement is much closer than, for instance, the difference between P B − Q B for the Hamiltonian system and the corresponding RMT prediction in Fig. 5(d), which is considerable evidence that the origin of the shift is indeed the randomization of the fluctuating part of the reduced energy eigenstates, as discussed in Sec. II B. FIG. 14. Diagrammatic proof of the measurement protocol. We use the diagramatic notation and calculus developed in Ref. [127] (see also Ref. [108]). With the definitions of the text, we have UA = i∈A ui, ρA = i∈A = ρi, OA = i∈A Oi, and accordingly for subsystem B. From second to third line, we use the 2-design identities of the local random unitaries ui (see also Eq. (G8) and Refs. [108,127]).
construction: For each experimental run, a set of local unitaries {u r i } i=1,...,N and time-evolution operator T is independently sampled and applied according to the experimental sequence shown in Fig. 2. Lastly, a single-shot computational basis measurement is taken. We thus have for an arbitrary r ∈ {1, . . . , M }. We drop the superscript (r) in the following.
To make progress, we evaluate the expectation E [ô] over the ensemble of time evolution operators (the disorder average) E T , the local random unitaries E U and projective measurements (the quantum mechanical expectation value) E QM step-by-step using the law of total expectation Here, E QM [ô|U, T ] denotes the quantum mechanical expectation value of the single-shot estimatorô for a fixed unitary U and a fixed time-evolution operator T . By definition, this is just the quantum expectation value of the observable O in the output state ρ f = U † T U ρ 0 U † T † U of our protocol, The key part of the proof is the evaluation of the expectation value over the local random unitaries U = i u i for a fixed time evolution operator T , As also visualized in Fig. 14, this requires several steps: We first rewrite as an expectation value of two 'virtual copies' of qubit i, using the identity Tr [AB] = 2 N Φ + N |A T ⊗ B|Φ + N for any two operators A and B. Here, we have defined |Φ + N = i |Φ + i as the tensor product of Bell states |Φ + i = 2 −1/2 (|00 + |11 ) on the doubled Hilbert space C 2 ⊗ C 2 . We now use the independence of the local random unitaries u i to completely factorize the expectation value E U over the local random unitaries U = i u i where du i denotes the Haar integral over the unitary group U (2). As shown in Refs. [41,107,135] [and also follows directly from Eq. (B10)], we can use the 2-design identities of the applied local random unitaries u i to evaluate the Haar integral. We find, Proposition 1. Consider a subsystem A ⊆ S with N A ≤ N qubits. Our aim is to estimate the PSFF K A using the estimator K A defined in Eq. (4). Then, for any , δ > 0, a total of M ≥Ṽ A δ 2 (H13) experimental runs (single shot estimates) suffice to ensure that the relative error of the estimator K A obeys | K A /K A − 1| ≤ with probability 1 − δ. Here, we defined the rescaled variancẽ where the sum extends over all subsystems B ⊆ A containing N B ≤ N A qubits. For the random matrix ensembles considered in App. B, we can determineṼ A explicitly. In particular, for CUE dynamics T (t = nτ ) with V from CUE, we find at the point of weakest signal, i.e. the dip time t = τ , V A = 10 N A − 1.