Probing Many-Body Quantum Chaos with Quantum Simulators

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 the 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 nonintegrable, but also the ability to measure novel observables beyond standard low-order correlation functions [4][5][6][7][8]. It includes recent measurements of the growth of entanglement entropies in quantum many-body systems [9][10][11][12] as well as of the decay of out-of-time-ordered correlation functions [13][14][15][16][17][18][19]. In this work, our interests lie in developing experimentally feasible probes of universal RMT predictions for the statistics of energy eigenvalues [1,[20][21][22][23][24] and predictions of ETH for the statistics of energy eigenstates [2,3,[25][26][27][28][29] of quantum chaotic manybody 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 ETH is violated.

In this paper, we identify the spectral form factor (SFF) and its generalization to partial SFF (PSFF) as observables of interest to reveal energy level and eigenstate statistics. The SFF is defined in terms of the time-evolution operator of the quantum many-body system of interest and provides us with statistics of energy levels [1]. The PSFF will be defined in terms of the time-evolution operator restricted to a subsystem of the many-body system and contains information on both the statistics of energy eigenvalues and energy eigenstates. We derive analytic expressions for the PSFF in Wigner-Dyson random matrix ensembles. More generally, in chaotic quantum many-body systems, 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][36][37][39][40][41][42][43][44][45][46][48][49][50][51][52][53][54][55][56][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 timeevolution operator. We emphasize that the present protocol is ancilla-free. This is in contrast to Ref. [58], where a measurement scheme for the SFF is 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 one 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 a 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 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 fE j g 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 define 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 quasienergies, 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 V 3 ¼ e -iH ðxÞ τ=3 e -iH ðyÞ τ=3 e -iH ðzÞ τ=3 :

Here, the Hamiltonians H ðx;y;zÞ contain nearest-neighbor interactions with strength J ¼ 3τ -1 and disordered transverse fields with strength h ðx;y;zÞ i ∈ ½-J; J: H ðx;y;zÞ ¼ J X N- ; and σ a [a ∈ ðx; y; zÞ] denote the Pauli matrices. We denote the number of spins with N such that D ¼ 2 N . An ensemble average is naturally performed by averaging over many FIG. 1. Illustration of the characteristic properties of the SFF and PSFF using the chaotic spin-1=2 Floquet model V 3 . (a) We display the SFF KðtÞ for the Floquet model V 3 with N ¼ 6 qubits as a function of time t. We observe characteristic features such as the ramp between t ∼ τ to t ¼ t H ¼ 2 N τ and a plateau for t > t H . (b) For the PSFF K A ðtÞ, we observe a ramp, a plateau and, in particular, a constant, additive shift of the PSFF compared to the SFF, which depends on the subsystem size N A of the subsystem A. We have chosen subsystems A from the middle of the total system. In both, the colored lines show the numerically calculated SFF and PSFFs, averaged over 8000 disorder realizations. In addition, we illustrate our measurement protocol (see Sec. I 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 (a) sketches the form of the SFF generically expected in a many-body localized model.

instances of Tðt ¼ nτÞ ¼ V n 3 , each with local disorder potentials h ðaÞ i sampled independently from the uniform distribution on ½-J; J.

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 the Appendix 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. Figure 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 manybody 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 quasienergy 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 are 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 ½jE i ihE i j denotes the reduced density matrix obtained after partial trace of the eigenstate jE i i of the Hamiltonian H (the Floquet time-evolution operator V) with energy (quasienergy) 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 jE i 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 one to extract these purities and overlaps, averaged over spectrum and ensemble, i.e., allows one 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 the 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 subsystemdependent 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 exhibit 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 fewbody 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 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

This can be noted for the CUE in 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.

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 ¼ j0ih0j with j0i ≡ j0i ⊗N . (ii) On this initial state, we apply local random unitaries

, 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) single-qubit 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 U † resulting in the final state ρ f ðtÞ ¼ U † TðtÞUρ 0 U † T † ðtÞU. (v) Last, we perform a single-shot measurement in the computational basis with outcome bit string s ¼ ðs 1 ; …; s N Þ with s i ∈ f0; 1g 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 bit strings s ðrÞ with r ¼ 1; …; M is collected.

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 Sec. 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 around 2 5N A experiments to reconstruct the full timeevolution 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 FIG. 2. Probing SFF and PSFF using randomized measurements. We present our protocol for the measurement of the SFF and PSFF using statistical correlations of local random unitaries applied at different times in a single experiment. We begin with a product state ρ 0 ¼ j0ih0j ⊗N . Before and after the time evolution TðtÞ, we apply random local rotations U ¼ ⊗ i u i and U † , respectively, where local unitaries u i 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 the Floquet evolution operator for time period τ. In the last step, a single-shot measurement is performed in the z basis to collect a bit string of the form s ¼ ðs 1 ; s 2 ; …; s N Þ with s i ∈ f0; 1g. This procedure is repeated M times, and M bit strings are collected to estimate the SFF and PSFF using Eqs. (7) and (8). The gray shaded region shows one possible choice of subsystem A.

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.

Last, 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. Section 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. Last, we summarize in Sec. VI with some concluding remarks and future directions.

II. PARTIAL SPECTRAL FORM FACTOR: ANALYTIC RESULTS

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

This section is organized as follows. In Sec. II A, we analyze the PSFF in standard Wigner-Dyson random matrix ensembles (see Appendix 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 involves only overlaps of distinct energy levels.

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 jtj > t ρ and ½δρ B ðtÞ jk ¼ 0 for jtj < t ρ (with some additional details in Appendix 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 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 B ½δρ

ρ , fluctuating around an average of zero (we do not require this behavior to persist over larger energy scales jE i -E j 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 or 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 size comparable to 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 much 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 A ðt ≫ t ρ Þ ¼ δP B =D A , where δP B ≡ D -1 P i Tr B ½δρ 2 B ðE i Þ (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 nontrivial 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 Appendix 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 γ ¼ P 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 Gaussian orthogonal ensemble (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πÞ. In Eq. ( 15), we also introduce 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

, giving the expected plateau value of P B =D A in Eq. (15). There are also two competing timescales 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

We note that this definition of Q B is equivalent to that in Sec. II A for random matrix ensembles, where the ensembleaveraged overlaps between distinct states are independent of their energies. Section III directly uses P B and Q B , with the implicit assumption that f Δ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 see that, at late times, the PSFF preserves the characteristic features of the SFF, such as the ramp and the Heisenberg time [as in Eq. ( 15) for fully chaotic systems]. However, there are non-negative subsystem-dependent parameters P B , δP B , and ΔP B (∼ g Δ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 show 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 nonthermal 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 turns out to be the primary nontrivial indicator of the properties of eigenstates in the PSFF.

In Appendix 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 in Sec. III. If these eigenstates are thermal, subsystem ETH requires the smooth and fluctuating parts to satisfy

Nonthermal 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

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) and ( 18), on relating the two using Eq. ( 16). In this context, we note that Ref. [80] observes subleading corrections to the random matrix prediction for eigenstates in 1D Floquet quantum circuits.

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 time-dependent feature due to the existence of local level correlations. The PSFF inherits the ramp, but the ramp couples only to the smooth, slowly varying part of the reduced energy eigenstates. The rapidly fluctuating part is left over as a nearly timeindependent 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 nonthermal 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 or 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 latetime 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 latetime quantum chaos) and their interdependence, across a transition or crossover between chaotic and nonchaotic phases. This could be done by tuning the parameters of a system (say, in a quantum simulator) between such phases and measuring PSFFs across different choices of subsystems of different sizes-analyzing the extent of delocalization of eigenstates in the absence of a ramp via the plateau value and, additionally, the extent of thermalization through the value of the shift if a ramp or other timedependent 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 nonergodic 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].

A. Example 1: Floquet system

The Floquet time-evolution operator V 3 has the same quasienergy 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 Appendix B). For the PSFF K A ðtÞ at N A ¼ 3 and very early times, we notice that the onset of the ramp takes a few initial periods to set, but eventually the PSFF follows the CUE prediction.

The closeness between the statistics of CUE and V 3 can further be seen from the average overlaps of reduced densities of eigenstates 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.

In conclusion, the SFF, PSFF, averaged purity, and overlaps match in the CUE and V 3 model, and, thus, we expect the form of the PSFF in Eq. ( 5) to hold for the model V 3 , after a small initial time period. We know from Eq. ( 12), for large Hilbert space dimensions, that Q B ≈ 1=D B and P B -Q B ≈ 1=D A . Therefore, utilizing Eq. ( 16), we find that ΔP B ¼ 0 and δP B ¼ Oð1=D A Þ for V 3 and the RMT models. The purity of the smooth part (of the form of Tr½Δρ 2

B ðEÞ) appears in the ramp part of the PSFF in Eq. ( 15), and, thus, we note that the ramp coefficient is approximately 1=D 2 for D A ≫ 1. 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 order 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 Appendix E.

B. Example 2: Hamiltonian system

As our second example, we consider a transverse field Ising model in the 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 the chaotic regime, and W ≫ J corresponds to the localized regime (for a similar model, see Ref. [64]). In Appendix 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 find that the PSFF preserves this flat feature of the SFF and has a subsystemdependent shift added over the SFF, as predicted in Sec. II C. In Figs. 4 and5, we present numerical results for the Hamiltonian model (19) in these two phases. For clarity, we use red 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 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 nonuniversal 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 shift, FIG. 3. Results for the Floquet V 3 model. (a) The SFF and PSFF are presented for N ¼ 6, N A ¼ 3, 4, 5 in red colors. In gray, we plot the same quantities in a CUE model. (b) The plateau value Kð∞Þ multiplied with the subsystem dimension D A is plotted in black circles and matches with the averaged purity P B plotted with red crosses. The average overlap Q B and the difference P B -Q B are presented in brown and green, respectively. We observe a perfect match with the respective quantities in CUE plotted in gray, indicating the same averaged eigenvalue and eigenstate statistics in CUE and V 3 . In the numerical computation, we take 8000 disorder realizations to perform ensemble averaging, and the subsystems A are chosen from the middle of the spin chain. 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 matches with the one for W ¼ J. As a guide, we plot 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 use 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.

Furthermore, to concretely discuss second moments of eigenstates, in Fig. 5(b), we present the averaged purity P B using crossed markers. We also plot 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 in Fig. 5 as functions of N A . Assuming ETH for the chaotic systems, we discuss 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 Appendix F 2 for more details on the numerical extraction of these orders). From Q B [Fig. 5(c)], we find ΔP B ¼ Oð1=D B Þ, and from P B -Q B ¼ δP B [Fig. 5(d)], we find δP B ∼ Oð1=D A Þ, confirming 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 Appendix F 3. For comparison, we plot 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 in Fig. 5. Following Eq. ( 16) from the Q B [Fig. 5(c)], we find ΔP B ¼ Oð1=D B Þ, and from P B -Q B ¼ δP B [Fig. 5(d)], we find δP B ∼ Oð1Þ ≫ Oð1=D A Þ. The localized phase is not expected to satisfy ETH, and, as discussed in the Sec. II C, we expect such a large shift in the PSFF in MBL FIG. 5. Results for the Hamiltonian model. (a) In linear scale, we present the SFF and PSFF for the chaotic phase (W ¼ J). (Sub)system sizes N A ¼ 6, 7 and N A ¼ N ¼ 10 are plotted with magenta, green, and red, respectively, for both the Hamiltonian model (with solid curves) and the GOE (with dashes). We observe differences in the PSFF for chaotic Hamiltonian and GOE. These differences are investigated in (b)-(d) through P B and Q B . We use red for the chaotic phase (W ¼ J) and gray for the GOE. For comparison, we also plot these quantities in the localized phase (W ¼ 10J) using blue. (b) We plot K A ð∞ÞD A using black circles, which matches with the corresponding average purity P B of the MBL and chaotic phase. (c) The average overlap Q B for MBL, chaotic, and GOE follow closely the behavior 1=D B . (d) The difference P B -Q B ≈ δP B , which encodes the shift of the PSFF, is larger for large disorders (MBL) compared to small disorders (chaotic). In the numerical computation, we take 200 Hamiltonians to perform ensemble averaging, and the subsystems A are chosen from the middle of the spin chain. FIG. 4. Results for the Hamiltonian model. In a log-log plot, we present the chaotic phase (W ¼ J) in red, MBL phase (W ¼ 10J) in blue, and the GOE in gray. In both phases, the SFF and PSFF are plotted for (sub)system sizes N A ¼ 6, 7 and

The SFF for the chaotic phase has the characteristic ramp and plateau and follows the GOE SFF at late times. The PSFF in this phase also has the shift, ramp, and plateau; we plot these in a focused linear scale plot in Fig. 5(a). The MBL phase shows a flat SFF and PSFF for all times t ≫ 0. The mean level spacings (i.e., the Heisenberg time) in the MBL phase and GOE are numerically rescaled to match to the one in the chaotic phase.

A. Useful results from unitary 2-designs

Unitary n-designs are ensembles of random unitary matrices, whose averages of polynomial moments of the order of 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 bit strings from the protocol are sampled from the Born probabilities jhsjU † TðtÞUj0ij 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

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 . Furthermore, the sum extends to all two-copy permutation operators

Here, the identity 1 i and the swap operator S i act as

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 turns out to be particularly useful to choose

where the sum extends to all bit strings s ¼ ðs 1 ; …; s N Þ with s i ∈ f0; 1g and jsj ≡ P i s i . For this choice, we obtain

with S ¼ ⊗ i S i ¼ P s;s 0 js 0 ihsj ⊗ jsihs 0 j. The swap operation S is the key operation to extract nontrivial 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 Appendix G. We first define for a fixed time-evolution operator TðtÞ K TðtÞ ≡ 4 -N Tr½TðtÞTr½T † ðtÞ ð24Þ

such that the ensemble (disorder) average KðtÞ ¼ K TðtÞ yields the SFF, according to the definition Eq. ( 1). Second, we show that K TðtÞ equals the survival probability of the Bell state jΦ þ N i under the dynamics generated by 1 ⊗ TðtÞ, i.e.,

To this end, we use the following identity for any two operators A and B on H:

which can be proven by inserting the definition of the Bell state jΦ þ N i ¼ 2 -N=2 P s jsi ⊗ jsi in terms of computational basis states. Equation (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 K TðtÞ ¼ E U ½K TðtÞ;U with K TðtÞ;U defined as

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, jhsjU † TðtÞUj0ij 2 is precisely the Born probability of finding a bit string s, in the computational basis measurement performed at the end of our measurement sequence when the state ρ f ðtÞ has been prepared (cf. 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 bit string s, ð-2Þ jsj 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 d KðtÞ [Eq. (7)]. Using Eq. ( 30), it follows that d KðtÞ converges to KðtÞ in the limit M → ∞. For finite M, statistical errors are governed by the variance of d KðtÞ and are discussed in the next section. In Appendix G, we extend our derivation to the case of the PSFF and illustrate the mapping between randomized measurements and the (P)SFF graphically.

Observing PSFF and SFF

We can bound the statistical errors of the estimator d K A ðtÞ [Eq. (8)] by its variance. As shown in Appendix H, we find that

where we drop the time argument for brevity. Here, K B denotes the PSFF defined in subsystem B, and the sum extends over all subsystems B ⊆ A. The variance of d KðtÞ [Eq. (7)] follows by taking A to be the full system. We obtain an expected relative error

with M experimental runs. As 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 Hamiltonian (Floquet) operators from Wigner-Dyson RMT, we can explicitly evaluate σ A (see Appendix 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 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= ffiffiffiffi ffi M p 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

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 ½σ

, 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 subsystem A but only on its size jAj ¼ 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 dataset, 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 jAj ðtÞ and its slope c jAj ð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 jAj ðt 2 ; t 1 Þ for these average quantities.

In Fig. 6(d), we display the numerically determined signalto-noise-ratio SNR½c jAj ð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 jAj ðt 2 ; t 1 Þ for subsystem sizes 1 < N A < N. We remark that we do not reach an enhancement ð

Þ is plotted for total system size N ¼ 6 and subsystem sizes N A ¼ 3, 4, and 6 as a function of the number of measurements M at time t=τ ¼ 5. For a fixed M, we perform 100 numerical experiments, each with M single shots, and present the average E A using colored lines. The gray lines represent corresponding errors in a model with CUE dynamics (calculated analytically in Appendix H). (b) The relative error E A as a function of subsystem size N A at two times t=τ ¼ 5 and t=τ ¼ 30 is shown. (c) Single PSFF with subsystem size N A ¼ 4 for two choices A ¼ ½2; 3; 4; 5 (red line), A ¼ ½1; 2; 5; 6 (green dashed line), and the average PSFF K jAj (black line) follow each other; the numbers in the ½Á Á Á denote qubit index. (d) For the observation of the ramp, we plot the SNR of the slope, SNR½c A ðt 2 ; t 1 Þ (in red) and SNR[c jAj ðt 2 ; t 1 Þ] (in black). Both SNRs are constructed from a single dataset of M ¼ 10 6 . As a guide to the eye, we also present in gray the SNR[c jAj ðt 2 ; t 1 Þ] when all the measurements in the averaged PSFF are done independently, i.e., when

from a single dataset are not independent. Nevertheless, Fig. 6(d) shows that the average PSFF K jAj 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 jAjð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.

B. Experimental imperfections

First, we consider an imperfect implementation of our measurement protocols, with errors arising from an erroneous decorrelation of the applied initial and final local random unitaries. We model such imperfection as the effective application of a unitary u i before and a unitary v i ¼ u † 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 Gaussian unitary ensemble (GUE) [1]. While the case η ¼ 0 corresponds to the ideal case, we display in Fig. 7 Second, 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 one 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 one 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 green, 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 show 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 d K A ðtÞ. They, thus, cannot 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 SFF [KðtÞ in Eq. ( 1)] is known to be a key diagnostic of many-body quantum chaos. In chaotic FIG. 7. Experimental imperfections and decoherence. We study effects of measurement errors and decoherence on the estimated SFF and PSFF K A ðtÞ using the example of the kicked spin V 3 with total system size N ¼ 4. In (a), we display the relative error

the estimated form factors induced by a decorrelation of local random unitaries applied before and after the time evolution up to the Heisenberg time t H , with strength η (see the text). In (b), we display the estimated SFF ðKÞ dec (blue dots) and PSFF ðK A Þ dec (blue squares) as a function of time in a system subject to global polarization with strength p ¼ 0.03 (see the text). For this type of decoherence, rescaling according to Eq. ( 38) allows one to recover the SFF (green dots) and PSFF (green squares) for unitary dynamics (red line). 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 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 ETH, we have shown in Sec. II that PSFFs in generic chaotic quantum many-body systems possess a characteristic shift-ramp-plateau structure [Eqs. (11) and (15)] and reveal crucial differences between thermal and nonthermal 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 test bed of the statistical properties of eigenvalues and eigenstates. Our protocol can be directly implemented in state-of-the-art quantum devices, based, for instance, on trapped ions [4,6], Rydberg atoms [7], and superconducting qubits [8,19], providing crucial experimental tools for the quantum simulation of many-body quantum chaos and the study of thermalization in closed quantum systems.

Our work can be generalized in various directions. First, while we have concentrated here on quantum simulators with local control realizing lattice spin models, our protocol can be also realized in collective spin systems with only global operations [118]. Second, while we have considered form factors which are second-order functionals of the time-evolution operators TðtÞ, partial restrictions of higherorder 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]. Third, we have focused on determining the properties of unitary quantum dynamics. Beyond that, our measurement protocol readily extends to the study of noisy quantum channels. This includes applications in the field of verification and benchmarking of quantum devices [89][90][91][92][93]121,122], as well as the investigation of noise-induced quantum many-body phenomena such as entanglement phase transitions [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].

(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 timescale [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 the following.

GUE model:

for 0 < t ≤ t H ;

GOE model:

where rðtÞ ¼ t H J 1 ð4Dt=t H Þ=ð2DtÞ with J 1 denoting the Bessel function of the first kind. The Heisenberg time t H , connected to the inverse spacing of adjacent energy levels, depends on the width of the Gaussian distribution of the matrix elements and marks the onset time of the plateau of the SFF. For the results presented in Sec. III, we fix it numerically, by matching plateau onset times for the Hamiltonian Eq. ( 19) and the GOE model. For the Floquet operators V from CUE or COE, one finds the following.

CUE model:

COE model:

Here, t H ¼ Dτ with τ to be identified with the period of the Floquet system to be modeled.

Proof of Eqs. (B1)-(B3)

We denote the basis of H consisting of eigenvectors of TðtÞ with jii (i ¼ 1; …; D). Furthermore, we fix an arbitrary product basis of H ¼ H A ⊗ H B as ja; bi with a ¼ 1; …; D A and b ¼ 1; …; D B . With TðtÞ ¼ YDðtÞY † , we rewrite Eq. ( 4) in these bases. Using the independence of eigenvalues and eigenvectors (Fact 1), we find

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 the orthogonal group can be formulated in terms of pair partitions, defined as follows.

Definition 1 (pair partitions).-For n ∈ N, (a) we denote with M O ð2nÞ the set of all pair partitions of f1; …; 2ng, partitioning f1; …; 2ng into n distinct pairs. Then, each pair partition m ∈ M O ð2nÞ can be uniquely expressed as ffmð1Þ; mð2Þg; …fmð2n -1Þ; mð2nÞgg ðB7Þ with mð1Þ < mð3Þ < Á Á Á < mð2n -1Þ and mð2i -1Þ < mð2iÞ for all i ∈ f1; …ng.

(b) We denote with M U ð2nÞ ⊆ M O ð2nÞ the set of all pair partitions of f1; …; 2ng which pair elements in f1; …; ng with elements fn þ 1; …; 2ng. Then, each partition m ∈ M U ð2nÞ can be uniquely expressed as ffmð1Þ; mð2Þg; …fmð2n -1Þ; mð2nÞgg ðB8Þ with mð1Þ < mð3Þ < Á Á Á < mð2n -1Þ and mð2i -1Þ ∈ f1; …; ng and mð2iÞ ∈ fn þ 1; …; 2ng for all i ∈ f1; …ng.

The following fact is shown in Ref. [127]. B6) for Y ∈ UðDÞ (GUE and CUE case). To perform the average over eigenvectors (green line), we remove the boxes Y and connect white decorations of Y (rhombi) with white decorations of Y Ã (rhombi) and black decorations of Y (circles and squares) with black decorations of Y Ã (circles and squares) in all possible ways, corresponding to the pair partitions m; n ∈ M U ð4Þ [108,128]. Summing over the resulting diagrams, weighted with the corresponding value of the Weingarten function, yields Eq. (B2). In all diagrams, each blue loop contributes a factor D A , and each red loop a factor D B .

Wg UðDÞ ðm; nÞ Y n k¼1 δ i mð2k-1Þ ;i mð2kÞ δ j mð2k-1Þ ;j mð2kÞ ðB10Þ with M U ð2nÞ ⊊ M O ð2nÞ the set of all pair partitions on f1; 2; …; 2ng which pair elements in f1; …; ng with elements fn þ 1; …; 2ng and Wg UðDÞ the Weingarten function on the unitary group UðDÞ.

In our case, we are interested only in the case n ¼ 2. As shown in Ref. [127], when m; n ∈ M O ð4Þ and D ≥ 2,

Furthermore, it holds for m; n ∈ M U ð4Þ and D ≥ 2 that w U eq ≡ Wg UðDÞ ðm; nÞ

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

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 in 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:

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), is taken to be the fluctuating part:

Up to this point, such a decomposition is always possible. We 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 cannot 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 do not require randomization of the correlators of δρ B ðEÞ between energies much 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 first perform a prototype calculation with simpler notation. Consider two functions fðEÞ and gðEÞ of a continuous variable E, with respective Fourier transforms fðtÞ and gðtÞ, both of which potentially vary over different realizations of the ensemble. We 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 P j;k e iE j t 1 -iE k t 2 , which reduces to KðtÞ at equal times t 1 ¼ t 2 ¼ t. The sum of phases D -2 P j;k e iE j ðtþt l Þ-iE k ðtþt r Þ in Eq. (C3) fluctuates strongly over different ensemble realizations at large t 1 and 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 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 and E k (closer than the scale of 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 satisfies Ωðt ≫ t Th Þ ≈ 0. The ensemble-averaged sum over E j and 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 : 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 argue that all cross terms with f ≠ g may be taken to vanish. When f ¼ D -1 B , the overlap becomes Tr B ½fg ¼ D -1 B Tr B ½g, 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 complete the evaluation of ΔK A ðtÞ using the prototype Eq. (C4) with f ¼ g ¼ ðΔρ B Þ ab and the expression in Eq. (C7) with

As the definition of Δρ B sets t l ; t r < t ρ , we have jT 12 j ¼ jtj þ 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 jtj 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 jtj. 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:

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.

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 l ≤ D that it has a significant probability of being found in.

It is useful to introduce an effective dimension D eff A ≤ D A ; l 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 l). 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

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 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 ðl=D eff A Þ, with the second term containing larger scale variations of its components. We 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

Being a linear combination of the macroscopic and microscopic contributions, the purity shows an initial decay with D eff A for small values of the latter and eventually a growth for larger values of

l correspond to pure states with P B ¼ 1. The parameters δP B and ΔP B satisfy the following order-of-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

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 decay only 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, nonthermal 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

We conclude that P B is a measure of delocalization of eigenstates, while δP B and ΔP B are probes of thermalization. Setting l ¼ D gives the results discussed in the main text for chaotic systems with fully delocalized eigenstates (Sec. II B 2). For fully localized systems, l ¼ Oð1Þ gives D eff A ¼ Oð1Þ, with P B ¼ Oð1Þ and δP B ¼ Oð1Þ ≲ ð1 -D -1 B Þ, 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 Refs. [27,28]). Denoting the eigenstates by jEi, there are two statements of the hypothesis: the diagonal statement stating that the reduced density matrix ρ B ðEÞ ¼ Tr A ½jEihEj 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 ½jE 1 ihE 2 j with E 1 ≠ E 2 to be small. We adapt these statements, in their subsystem-dependent version (which does not 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. Equations (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 in Ref. [29], which we refer to as the "original conjecture" in this Appendix, and we 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 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ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D B =D A p Þ. As Ref. [29] notes, on account of the inequality fTr B ½ðx 2 Þ 1=2 g 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 Haar-random superpositions of energy eigenstates derived in Refs. [131,132]; but, in the process of the derivation in the latter, constraints in terms of purity similar to Eq. (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 has to restrict the subsystem-dependent 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 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 comparable to t -1 s . We 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

of Eq. ( 13). Equation (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

ρ (e.g., neighboring levels) so that the second term is approximately Tr A ½Δρ 2

A ðE 1 Þ. 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 g Δ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ð ffiffiffiffiffiffi ffi D A p =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 jE 1 i and jE 2 i 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 Þ ¼ jE 1 ihE 1 j and ρðE 2 Þ ¼ jE 2 ihE 2 j 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 two-point correlation, which crucially requires even largescale nonuniformities 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ð ffiffiffiffiffiffi ffi D A p =DÞ, justifying the above constraint. Similarly, the operator qðE 1 ; E 2 Þ ¼ jE 1 ihE 2 j for such eigenstates has OðD -1 Þ matrix elements with random phases, giving Oð ffiffiffiffiffiffi ffi D A p =DÞ matrix elements after the partial trace and, therefore, the off-diagonal statement Eq. (D5). Such a picture of random energy projector matrix elements of comparable magnitudes is reminiscent of Berry's conjecture for chaotic wave functions [134] (as well as other related statements, e.g., Refs. [2,[94][95][96]), which can be 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 jE 1 i that is randomly (but not necessarily uniformly) distributed only over a set of approximately l ≤ D physical basis vectors, with negligible support outside this set. Its density matrix ρðE 1 Þ then has an l × l block (after suitably permuting rows and columns) of nonvanishing elements each of typical magnitude Oðl -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 and columns.

The behavior of ρðE 1 Þ under a partial trace depends on the choice of the subsystem A. We 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 l 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 (l ¼ 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 l 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 l 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 call such a subsystem A an unbiased subsystem (from the point of view of the eigenstate of interest). Then, ρ B ðE 1 Þ has an approximately ðl=D A Þ × ðl=D A Þ nonvanishing block with non-negative diagonal elements of magnitude OðD A =lÞ and off-diagonal elements of typical magnitude Oð ffiffiffiffiffiffi ffi D A p =lÞ 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 =lÞ 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 approximately l=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 ½ðD A =lÞ þ OðD A =lÞ, giving

Thus, we can extract information about the extent of delocalization l 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 jE 2 i that is also distributed only across l basis vectors (but not necessarily the same ones or in the same way as jE 1 i) 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 l 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 sufficient only for fully, uniformly delocalized random eigenstates as in Appendix B).

We now consider two illustrative extreme cases of fully overlapping (thermal) and fully nonoverlapping (nonthermal) eigenstates. In both cases, we are interested in

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.

(i) Thermal eigenstates.-If jE 1 i and jE 2 i are distributed in a similar manner across the same basis vectors, then ½ρðE 1 Þ -ρðE 2 Þ again has an l × l block structure, with random Oðl -1 Þ off-diagonal elements within the block. However, the diagonal elements, being differences of random Oðl -1 Þ non-negative numbers, also have at most Oðl -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 ½ρ

consistent with diagonal ETH [Eq. (D4)] in subsystem B. For the overlap Tr B ½ρ B ðE 1 Þρ B ðE 2 Þ, the positivity and normalization of the diagonal matrix elements of each reduced density matrix ensure that their contribution is of the form ½ðD A =lÞ þ OðD A =lÞ. The products of the offdiagonal matrix elements add up with random phases, leading to a negligible Oðl -1 Þ contribution. We therefore have

which is the analog of off-diagonal ETH [Eq. (D5)] for subsystem A. (ii) Nonthermal eigenstates.-In the nonthermal case, jE 1 i and jE 2 i 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

A Þ with some OðD A =lÞ 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 l × l 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 < l, due to the assumption on subsystem A. The reduced energy eigenstates in subsystem B are pure basis states when D A ¼ l 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 ≤ l. Both bounds are saturated by thermal eigenstates.

For greater completeness of the present discussion, we should account for a more typical choice of subsystem Aone that mixes the zero and nonzero elements of these eigenstate reduced density matrices on performing the partial trace over A. We consider such a typical subsystem to have an effective dimension D eff A ≤ D A , corresponding to the typical number of nonzero 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 ¼ l=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) and (D18) in terms of D eff A gives the constraints discussed in Appendix D 1.

APPENDIX E: TIME-REVERSAL SYMMETRIC FLOQUET THERMALIZATION

In this Appendix, we consider another Floquet model of a periodically kicked spin-1=2 system. We consider one period of duration τ to be , 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].

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 are 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 secondorder 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 the 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 approximately 1=D 2 A , as in the RMT models.

Comparison of the PSFF shift and δP B

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 now argue that it is reasonable to take

for our purposes. In Eq. (F3), by subsystem ETH, the former is OðD

. 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 approximately ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Dt H =t 0 p 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 > ffiffiffiffi D p as well. A minor additional effect that improves this approximation is that, for large D A , the coefficient D B g ΔP B of the second term would be small, though still Oð1Þ, from Fig. 12 (as D B g ΔP B ∼ ½D B Q B -1). 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.

APPENDIX G: DERIVATION OF THE MEASUREMENT PROTOCOL

In this Appendix, we show that our measurement protocol indeed allows us to measure the PSFF and SFF. We generalize the proof for the SFF in the main text (Sec. IV) to the PSFF K A ðtÞ and provide additional mathematical details. Our aim is to prove that d K A ðtÞ, as defined in Eq. (8), is an unbiased estimator of K A ðtÞ, i.e., E½ d K A ðtÞ ¼ K A ðtÞ, where E comprises the expectation value taken over the ensemble of time-evolution operators (the disorder average) E T , the local random unitaries E U , and projective measurements E QM . Note that we use in this Appendix E T in place of Á Á Á to denote the expectation over an ensemble of time-evolution operators.

We consider a quantum system S consisting of N qubits with Hilbert space H ¼ ðC 2 Þ ⊗N of dimension D ¼ 2 N and A ⊆ S of N A qubits with dimension D A ¼ 2 N A . From the r ¼ 1; …; M (single-shot) repetitions of our protocol with outcome bit strings fs ðrÞ g r¼1;…;M , we obtain the estimator d K A ðtÞ as defined in Eq. (8). For simplicity of notation, we drop the time argument in the following in this Appendix.

As a first step, it is most convenient to reformulate Eq. ( 8) as an average over r ¼ 1; …; M single-shot estimates ôðrÞ of an observable O ¼ ⊗ i O i with O i ¼ j0ih0j -1=2j1ih1j for i ∈ A and O i ¼ 1 i for i ∉ A:

Second, we note that the outcome bit strings fs ðrÞ g r¼1;…;M of the M repetitions of our measurement protocol are identically and independently distributed by construction: For each experimental run, a set of local unitaries fu r i g i¼1;…;N and time-evolution operator T is independently sampled and applied according to the experimental sequence shown in Fig. 2. Last, a single-shot computational basis measurement is taken. We, thus, have

for an arbitrary r ∈ f1; …; Mg. We drop the superscript (r) in the following.

To 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 ½ ôjU; 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 † TUρ 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 hΦ þ N jA T ⊗ BjΦ þ N i for any two operators A and B. Here, we define jΦ þ N i ¼ ⊗ i jΦ þ i i as the tensor product of Bell states jΦ þ i i ¼ 2 -1=2 ðj00i þ j11iÞ 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 ¼ ⊗

where R 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

To arrive at the last line, we use ρ i ¼ j0ih0j and that O i ¼ j0ih0j -1=2j1ih1j for i ∈ A and O i ¼ 1 i for i ∉ A. Inserting this into Eq. (G6), we find

Taking finally the ensemble (disorder) average over timeevolution operators, we find

Thus, we see that c K A is an unbiased estimator of K A ðtÞ.

APPENDIX H: STATISTICAL ERROR ANALYSIS AND REQUIRED MEASUREMENT BUDGET

As described in the main text [Eq. ( 8)], we obtain an estimate of the PSFF K A (for notational simplicity, we drop the time argument in this Appendix) from r ¼ 1; …; M (single-shot) repetitions of our protocol with outcome bit strings s ðrÞ via

Here, ôðrÞ is a single-shot estimate of an observable O ¼ ⊗ i O i with O i ¼ j0ih0j -1=2j1ih1j for i ∈ A and O i ¼ 1 i for i ∉ A, as defined in Appendix G. We show in Appendix G that c K A is an unbiased estimator of the PSFF K A , i.e., that E½ c K A ¼ K A with the expectation value taken over the ensemble of time-evolution operators (the disorder average) E T , the local random unitaries E U , and projective measurements E QM . The statistical error of c K A and its convergence to K A is controlled by its variance

for any r ¼ 1; …; M. Here, we use that the individual single-shot estimates ôðrÞ are statistically independent and FIG. 14. Diagrammatic proof of the measurement protocol. We use the diagrammatic notation and calculus developed in Ref. [128] (see also Ref. [108]). With the definitions of the text, we have

, and accordingly for subsystem B. From the second to the third line, we use the 2-design identities of the local random unitaries u i [see also Eq. (G8) and Refs. [108,128]].

identically distributed by construction. We drop the superscript (r) in the following. We can evaluate Eq. (H2) using the law of total variance [136]:

To arrive at the second expression, we employ the definition of the (conditional) variance Var½XjY ¼ E½X 2 jY -E½XjY 2 for any two random variables X and Y and use then that various terms cancel out. As shown in Appendix G, the last term in Eq. (H3) simply yields

We, thus, concentrate on the first term in Eq. (H3). The quantum mechanical expectation value E QM ½ ô2 jT; U of the squared single-shot estimate ô evaluates, for fixed T and U, to

Next, we evaluate the average over local random unitaries.

With O replaced by O 2 , we follow the calculation presented in Appendix G-we first rewrite Eq. (H5) as an expectation value on two copies:

Factorizing the average over local random unitaries, we find

with R du i the Haar integral over the unitary group Uð2Þ. Using Eq. (G8), for O i → O 2 i , we find

Inserting this into Eq. (H6), we obtain

Taking the ensemble (disorder) average over time-evolution operators T, we get

This finally yields

Given the variance Var½ c K A ¼ Var½ ô=M of our estimator, Chebyshev's inequality asserts that

for any ϵ > 0. This allows one to rigorously obtain an estimate for the required number of measurements M to achieve a certain relative error (for a similar treatment, see, e.g., Refs. [42,43]). Proposition 1.-Consider a subsystem A ⊆ S with N A ≤ N qubits. Our aim is to estimate the PSFF K A using the estimator c K A defined in Eq. ( 4). Then, for any ϵ, δ > 0, a total of

experimental runs (single-shot estimates) suffice to ensure that the relative error of the estimator c K A obeys j c K A =K A -1j ≤ ϵ with probability 1 -δ. Here, we define the rescaled variance

where the sum extends over all subsystems B ⊆ A containing N B ≤ N A qubits.

For the random matrix ensembles considered in Appendix 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 ¼ τ, ṼA ¼ 10 N A -1.