Probing scrambling using statistical correlations between randomized measurements

We propose and analyze a protocol to study quantum information scrambling using statistical correlations between measurements, which are performed after evolving a quantum system from randomized initial states. We prove that the resulting correlations precisely capture the so-called out-of-time-ordered correlators and can be used to probe chaos in strongly-interacting, many-body systems. Our protocol requires neither reversing time evolution nor auxiliary degrees of freedom, and can be realized in state-of-the-art quantum simulation experiments.


I. INTRODUCTION
Recent developments in quantum simulation have enabled the remarkable ability to interrogate and control atomic, molecular, and ionic degrees of freedom (d.o.f.) in lattice experiments with single-site resolution [1][2][3][4]. In atomic Hubbard models with bosonic and fermionic atoms in optical lattices, a quantum gas microscope provides us with single-shot spatial-and spin-resolved images of atomic densities. By averaging over many images, this allows one to extract spatial and spin equal-time correlation functions, which reveal unique properties of (non)equilibrium quantum phases [5][6][7][8]. In spin models, as realized with trapped ions [9,10], Rydberg atoms [11][12][13][14], and superconducting qubits [15][16][17][18], the state of the spins (qubits) can be measured in a given standard basis in single-shot measurements with high fidelity and with high repetition rates. Building on these tools, we develop below quantum protocols to measure many-body observables from analyzing statistical cross-correlations between such quantum images representing different runs of an experiment. We see that this provides us with simple, generic, and robust techniques to extract many-body observables, which are challenging to access otherwise within existing experimental setups. In particular, we develop novel protocols for out-of-time-ordered correlators (OTOCs), which are time-dependent quantities that cannot be measured directly as a standard time-ordered correlation function. OTOCs represent a key quantity to diagnose quantum chaos and enable one to understand how quantum information propagates, and "scrambles" [19], in close connection to the notion of entanglement spreading [20][21][22][23].
OTOCs have been introduced to characterize quantum dynamics, described by a unitary time-evolution operator UðtÞ in terms of the complexity of Heisenberg operators WðtÞ ¼ U † ðtÞWUðtÞ. For chaotic dynamics, even an initially "simple" and local Hermitian operator W rapidly becomes complex and nonlocal. As a consequence, after a short time, WðtÞ is delocalized and no longer commutes with an initially nonoverlapping local operator V. The degree of noncommutativity, or equivalently the scrambling of WðtÞ, is quantified by the OTOC, which takes the form with ρ the initial quantum state. Note that this definition ensures that OðtÞ ¼ 1 when WðtÞ and V commute. In the following, we focus on the "infinite-temperature" OTOC [22][23][24][25][26] for which ρ ¼ I=N H , with I the identity operator and N H the Hilbert space dimension. In Sec. V, we discuss extensions of our approach to thermal states. The time dependence of OðtÞ can distinguish between different classes of scrambling, ranging from "fast scrambling" in models with holographic duals [27][28][29][30][31][32] and chaotic many-body spin systems [19,22,23,33,34], to "slow scrambling" characteristic of many-body localization (MBL) [24][25][26]35]. These theoretical insights raise the question of how to experimentally measure OðtÞ, despite the peculiar time order inherent in its definition. A first option to measure OðtÞ consists in implementing time-reversal operations [21,[36][37][38][39] or using auxiliary quantum systems [40]. The first measurements of OTOCs were realized using this approach in systems with few d.o.f. [41], but also in a trapped ion setup with infinite-range interactions [42]. However, protocols based on realizing time-reversal operations remain an experimental challenge for many experimental platforms-like Hubbard systems or systems with local interactions. For such protocols, recent studies have also shown that decoherence can "mimick" the effect of scrambling, and they have developed [39,43] and realized [44] implementations involving auxiliary d.o.f. to distinguish the two effects.
In contrast, a unique feature of our protocols to measure OTOCs via statistical correlations is that they do not rely on time-reversal operations nor the presence of an ancillary system. In addition, we show that OTOCs extracted from statistical averages provide the key advantage to be naturally robust against various forms of decoherence and experimental noise, including depolarization and readout errors. As a consequence, our protocols can be realized in any state-of-the-art AMO [2,3,45] or superconducting qubit platforms [4], and they can be used as experimental probes of scrambling in many-body systems.
Our work presents two key results related to two protocols. First, we introduce the global protocol and demonstrate the exact equivalence between the OTOC OðtÞ, as defined above, and the statistical correlations obtained from initial states, which are randomized with a global unitary operator u for the total many-body system. We then present the local protocol, which consists of an experimentally simpler approach for spin systems, where the initial state is randomized with local unitary operations and where the statistical correlations also give access to OðtÞ.
This paper is organized as follows. Section II presents the main results of this paper by describing the two protocols to measure OðtÞ with local and random unitaries. Section III gives different physical examples, accessible to current AMO and solid-state platforms. Finally, we discuss in Sec. IV the role of statistical errors and imperfections, and we identify, in particular, the different situations (depolarization, readout errors) where the protocol is not affected.

II. PROTOCOLS MAPPING STATISTICAL CORRELATIONS TO OTOCs
In this section, we present and illustrate both the global and local protocols to measure OTOCs via random measurements. We consider a system S associated with a Hilbert space of dimension N H . This system can be, for example, a set of atoms described by a Hubbard model or an ensemble of spin-1=2 as shown in Fig. 1(a). In the following, we also assume the operator V to be unitary and the operator W to be Hermitian and traceless [TrðWÞ ¼ 0]. Note that these conditions do not restrict the ability of (a) (c) OTOCs to describe scrambling: For spin systems, we consider, as examples, W and V to be local Pauli operators, which are particularly relevant in this context [22,23]. We also give examples below that are relevant to probe scrambling in Hubbard systems. Our two protocols are illustrated in Fig. 1. The first protocol described in Sec. II A relies on global random unitaries u from the circular unitary ensemble CUEðN H Þ [46] or from a unitary 2-design [47]. As illustrated below for a Bose-Hubbard (BH) system, such random unitaries can be realized in generic interacting models using timedependent disorder [48][49][50]. The second protocol presented in Sec. II B considers spin systems with individual spin control, which allows us to drastically simplify the experimental task by replacing the global random unitaries by local random unitaries u ¼ u 1 ⊗ u 2 ⊗ … ⊗ u N , u i ∈ CUEð2Þ, acting on single spins i ¼ 1; …; N. Note that such local random unitaries have recently been realized with high fidelity with trapped ions [10,51].
We find it convenient to present both protocols as experimental recipes to measure the statistical correlations. In each case, we then mathematically relate the correlations to OðtÞ.

A. Global protocol
Experimental protocol.-The protocol consists of the following steps, as illustrated in Fig. 1(a).
(i) We prepare an arbitrary state jk 0 i and apply a global random unitary u to obtain jψ u;k 0 i ¼ ujk 0 i. The randomized state jψ u;k 0 i is the starting point for two independent experiments: (ii.a) In the first experiment, we evolve the system in time with UðtÞ and perform a measurement of the expectation value of W. The time-evolution operator UðtÞ can be generated, for instance, from a static Hamiltonian UðtÞ ¼ e −iHt , from Floquet evolution, or from a quantum circuit operating on qubits, depending on the type of system under study. We repeat steps (i) and (ii.a) with the same random unitary u to measure hWðtÞi u;k 0 ¼ hψ u;k 0 jWðtÞjψ u;k 0 i, as illustrated in Fig. 1(a). (ii.b) In the second experiment, we prepare again jψ u;k 0 i and apply the unitary V. This operation is followed by the time evolution with UðtÞ and a measurement of W. We repeat this sequence to obtain hV † WðtÞVi u;k 0 ¼ hψ u;k 0 jV † WðtÞVjψ u;k 0 i, as shown in Fig. 1(a). (iii) Finally, we repeat steps (i) and (ii) for different random unitaries. The OTOC OðtÞ, as defined in Eq. (1), is then obtained from the statistical correlations between the measurement outcomes hWðtÞi u;k 0 and hV † WðtÞVi u;k 0 of (ii.a) and (ii.b), respectively. Here, Á Á Á denotes the ensemble average over random unitaries u, and D ðGÞ ¼ hWðtÞi 2 u;k 0 is a normalization term. Proof and illustration.-Equation (2) can be proven using the 2-design identities, which provide analytical expressions for the statistical correlations between the matrix elements of u [52], with δ the Kronecker delta. In order to simplify the proofs, in this work we use a diagrammatic representation of Eq. (3), where the contractions of the different indices are represented by lines [53], as shown in Fig. 2(a). This representation allows us to prove in Fig. 2  . Finally, to conclude our proof, we use the same identity with V → I to prove that the denominator in Eq. (2) reduces to cTr½W 2 ðtÞ.
As an illustration, we present in Fig. 1(b) the intuitive physical picture behind this result [based on realizing UðtÞ for the kicked Ising model; see caption and text below]. At is not affected by the operator V ¼ σ z 1 , which distinguishes the two initial states jψ u;k 0 i and Vjψ u;k 0 i. Indeed, ½V;Wð0Þ ¼ 0 and hence hWð0Þi u;k 0 ¼ hV † Wð0ÞVi u;k 0 , implying perfect correlations, i.e., Oð0Þ ¼ 1 [up to shot noise errors; see below and Fig. 1(b), left panel]. At later times (middle panel), when ½V; WðtÞ ≠ 0 due to the spreading of WðtÞ, the value of hV † WðtÞVi u;k 0 becomes decorrelated from hWðtÞi u;k 0 . Note that in analogy to our approach, the distance between a quantum state and a (physical) copy, which is perturbed at t ¼ 0 by the operator V, has been proposed to numerically detect scrambling [54,55].
Measurement budget and statistical errors.-In an experiment, a finite number of N u random unitaries is realized to measure OðtÞ based on Eq. (2). Furthermore, the operator W is measured via a finite number of projective measurements N M per unitary. The total number of measurements leading to an estimation of the OTOC is therefore N meas ¼ N u N M . The finite value of N meas will lead to statistical errors. To illustrate this aspect, we compare in Fig. 1(b) the time evolution of OðtÞ with the estimation obtained from a finite realistic number of measurements [10] (circles with statistical error bars). We analyze the role of statistical errors in more detail in Sec. IV.
Implementation of random unitaries.-One of the experimental challenges in the protocol presented above consists in generating, with high fidelity, global random unitaries u satisfying the 2-design properties. In quantum simulators implementing Hubbard or spin models, this can be done using random quenches based on time-dependent disorder potentials (see Refs. [49,50] and the example in Sec. III A). We now proceed to describe an experimentally significantly simpler protocol for spin systems, which only requires us to generate local random unitaries acting on individual spins. These unitaries can be realized by combining local rotations along a fixed axis of the Bloch sphere, say, the z axis, with global rotations along an orthogonal direction, for instance, the x axis [10,51]; they can therefore be implemented in present qubit experiments with single-site control, e.g., with trapped ions [2], Rydberg atoms [3], or superconducting qubits [4].

B. Local protocol
The protocol.-We now describe our protocol based on local unitaries. The main difference compared to the protocol presented in Sec. II A is that we need to consider an ensemble E n ¼ fjk 0 i; …g of initial states, instead of a single one jk 0 i, in order to obtain a mapping between statistical correlations and OTOCs. The states jk s i ¼ jk ð1Þ s ; k ð2Þ s ; …i that we consider are written as product states in a standard fixed basis, i.e., k ðiÞ s ¼ ↑; ↓, and can be easily prepared in an experiment with single-site control.
The protocol consists of the following steps: (i) We prepare the first initial state jk 0 i, (s ¼ 0), and apply the local unitary u ¼ u 1 ⊗ …u N , cf. Fig. 1(c), to prepare jψ u;k 0 i. Here, each local unitary u i is drawn independently from a unitary 2-design, such as the circular unitary ensemble CUE(2). (ii.a) In the first experiment, we measure W after evolving the system in time with UðtÞ.
Repeating steps (i) and (ii.a) with the same random unitary u gives us access to In the second experiment, we proceed as in (ii.a), but we apply V before the time evolution with UðtÞ. This method gives us access to hV † WðtÞVi u;k 0 . (iii) We repeat steps (i) and (ii.a) for the other initial states jk s i ∈ E n , using the same random unitary u, in order to obtain hWðtÞi u;k s , for each k s ∈ E n . (iv) Finally, we repeat steps (i)-(iii) for different random unitaries. From these measurements, we can construct the statistical correlations with D ðLÞ n ¼ P k s ∈E n c k s hWðtÞi u;k s hWðtÞi u;k 0 , and weights c k s . Measurement budget and convergence aspects.-For this protocol, the number of measurements is N meas ¼ N u N M ðjE n j þ 1Þ, with jE n j the cardinal of E n , i.e., the number of initial states k s . In the following, we show how to choose the ensembles E n and weights c k s so that the correlations O n ðtÞ represent a converging series, indexed by n, of "modified OTOCs" approximating OðtÞ. The low-order OTOCs O 0;1;… ðtÞ, which correspond to small numbers of initial states (jE n j ¼ 1; 2; …) to sample and are thus the easiest quantities to access experimentally, provide generically good approximations of OðtÞ. We also prove, in the other limit n ¼ N, the exact relation O N ðtÞ ¼ OðtÞ.
Introducing the modified OTOCs.-Here, for concreteness, we consider V to be a Pauli operator on the first site i ¼ 1. For each value of n ¼ 0; …; N, we then define the ensemble E n as the set of all 2 n configurations jk s i, such that only the states of the first n spins can differ from the ones of the reference state jk 0 i, i.e., k ði>nÞ s ¼ k ðiÞ 0 for i > n. For instance, for n ¼ 0 (n ¼ 1, respectively), which we study in detail below, the ensemble is represented by a single state (just two states). As proven in the Appendix A, by choosing the weights c k s ¼ ð−1=2Þ d½k 0 ;k s , with d½k 0 ; k s the Hamming distance (the number of spin flips between jk 0 i and jk s i), we can relate the statistical correlations O n ðtÞ to modified OTOCs, which converge to OðtÞ for n ¼ N [56]. Here, the sums in the numerator are performed over all partitions A, which include the set B n ¼ f1; …; ng of the first n spins (for is "resolved" so that we can expect a better approximation to OðtÞ, and so on for n ¼ 2; 3; …. Note that our construction of the sets E n can be generalized easily to other positions of V and also to multisite operators. An illustration of this protocol is shown in Fig. 1(d) [compare to Fig. 1(b)], where we show O n ðtÞ (solid lines) and the corresponding statistical correlations obtained by simulating the protocol numerically. For n ¼ 0, the modified OTOC captures the scrambling time as OðtÞ but saturates to a nonzero value at long times. For n ¼ 1; 2; 3; …, the values of O n ðtÞ are in good quantitative agreement with OðtÞ. We also note that for short times, the local protocol has an advantage compared to the global protocol in terms of statistical errors, as we explain below.

III. REALIZATIONS OF THE PROTOCOLS IN DIFFERENT PHYSICAL SCENARIOS
This section is devoted to physical examples that can be accessed with our protocol. In Sec. III A, we show how to apply the global protocol in an atomic Bose-Hubbard system. We then focus on the local protocol and analyze, for chaotic (Sec. III B), many-body localized (Sec. III C), and long-range spin models (Sec. III D), the behavior of modified OTOCs O n ðtÞ. We analyze, in particular, both via analytical models and numerical simulations, the convergence properties of O n ðtÞ to OðtÞ.
A. Implementation of the global protocol in a Bose-Hubbard chain We now present an example to illustrate the different steps of the protocol with global unitaries and consider the situation of scrambling dynamics of the BH chain [21], with UðtÞ ¼ expð−iH BH tÞ, with a i bosonic operators and n i ¼ a † i a i . We consider here the OTOC dynamics for the unitary V ¼ expð−iπa † 1 a 1 Þ and the traceless observable W ¼ n jþ1 − n j .
We first illustrate the mapping Eq. (2)  Assuming here no projection noise N M → ∞, we obtain a very good agreement between the two quantities. Note, as in the case of the spin models described below, the is sampled from a uniform distribution of width 2J, and the quench time is T ¼ 1=J. For all panels, we consider N ¼ 8 lattice sites with jk 0 i ¼ j10101010i in the number basis, and U int ¼ 2J. In panels (b)-(d), error bars at 2 standard deviations are calculated using the jackknife resampling method. Here, we consider N M → ∞. characteristic "scrambling time" for the OTOCs varies essentially linearly with the position of the operators, showing the existence of a "butterfly" velocity.
We now discuss the physical realization of the global random unitaries u. Unitaries satisfying the required 2-design properties can be generated using the same Hamiltonian H BH subject to a sequence of η random quenches [49,50], with Δ ðmÞ j a random disorder potential, which is reinitialized for each quench m, and the quench time T. As shown in Refs. [49,50], such random quenches efficiently generate, in each particle number sector, unitaries satisfying the required 2-design properties of the CUE, after a time ηT ≈ N. This case is illustrated in Figs. 3(c) and 3(d) for two different operator W positions: For η ≥ N, the generated u converges to the CUE (with respect to the required 2-design properties), and therefore, the measured statistical correlations coincide within the statistical error bars with OðtÞ.

B. Chaotic dynamics with modified OTOCs
In the rest of this section, we focus on spin models with OTOCs measured by local unitaries. To illustrate the ability of the modified OTOCs O n ðtÞ, Eq. (6), to approximate OðtÞ, we first consider scrambling in chaotic spin models, which is characterized by two key features: The support of an operator WðtÞ that is initially localized grows ballistically with a "butterfly velocity" v B , and the operator front traveling at v B broadens diffusively [19,22,23]. The ballistic growth can be captured by a simple phenomenological model, which assumes that UðtÞ takes the form UðtÞ ¼ U½LðtÞ 1 ⊗ … ⊗ U½LðtÞ N=LðtÞ , where the Haar random unitaries U½LðtÞ ∈ CUEð2 LðtÞ Þ describe scrambling on a linearly growing scale LðtÞ ¼ 1 þ floorðv B tÞ. For W ¼ σ z j and V ¼ σ z 1 , we obtain, in leading order in LðtÞ ≫ 1 [33], and, as shown in Appendix B, Here, OðtÞ represents the average OTOC over the unitaries U, whereas the expression for O n ðtÞ corresponds to including the sampling over U in the ensemble average Á Á Á [57]. All OTOCs OðtÞ ¼ O 0;1 ðtÞ ¼ 1 thus coincide at short times when WðtÞ and V commute exactly. At the scrambling time t B ¼ r=v B (r ¼ j − 1), OðtÞ and O 0;1 ðtÞ exhibit a sharp drop. For t > t B , OðtÞ and O 1 ðtÞ are exponentially suppressed, while O 0 ðtÞ converges to 1=3.
We now confront these analytical predictions with numerical simulations of the kicked Ising model, which is an example of a chaotic spin model [22], with m a positive integer and T the period of the Floquet system. Throughout this work, we use open boundary conditions (OBC). The results are shown in Fig. 4, where we compare the modified OTOCs for n ¼ 0, 1, 2, 3, 4 with OðtÞ. In panel (a), corresponding to an operator position j ¼ 4, the scrambling dynamics described by OðtÞ is fast in the sense that it occurs at a time t B ∼ j=J. This behavior is qualitatively captured by the first modified OTOC, with fast decay at the scrambling time and saturation at the predicted value 1=3 of our phenomenological model. Interestingly, the convergence of O n ðtÞ to OðtÞ is achieved for small values, n ≳ 1. In panel (b), we represent the same quantities for a distant operator j ¼ 7. In this case, the dynamics of the OTOC includes a long-time slow behavior, which we attribute to the diffusive character of the operator WðtÞ [22,23]. This additional complexity of the operators [compared to ballistic spreading as in panel (a)] is captured for spatial resolutions n ≳ 4. In panel (c), we show another example of deviation from ballistic scrambling, with strong oscillations of the OTOCs, which are quantitatively captured for n ≥ 2; 3. Finally, we show in panel (d) the different regimes of scrambling as a function of the transverse field h x for a fixed large time Jt ¼ 35.2. All modified OTOCs identify a region of fast scrambling around h x ∼ J. Interestingly, the modified OTOCs with n ≥ 1 also provide an excellent approximation to OðtÞ in a regime of slow scrambling, h x ≥ 1.25J. Finally, for h x ≤ 0.75J, an increased resolution is necessary to access OðtÞ. This example shows that the required resolution n to access OðtÞ up to a given error depends on the type of evolution realized by UðtÞ and is generically small. Note that these results also suggest that the convergence properties of the series O n ðtÞ can be useful to identify different regimes of scrambling.

C. MBL dynamics with modified OTOCs
As a second example, with the opposite type of scrambling, we consider MBL, which is the paradigmatic example of a closed quantum system that does not thermalize [58,59]. As a key signature, the decay of OðtÞ is slow: It occurs at a characteristic time t B that scales exponentially with the distance r between the support of W and V at t ¼ 0 [24,26], which has to be contrasted to the linear increasing t B ∝ r of chaotic systems.
To begin our analysis, we first calculate the modified OTOC O 0;1 ðtÞ for the phenomenological l-bit model [60,61]

described by the Hamiltonian
with h z i random fields, J R;ij random interaction strengths that are taken uniformly in ½−J z ; J z , and ξ the localization length. Here, we only consider 2-body interaction terms, which is sufficient to show that MBL exhibits slow scrambling [24,26]

D. Information scrambling by long-range interactions with modified OTOCs
So far, we have considered examples where interactions are local, with analytical models supporting the statement that low-order n modified OTOCs provide good approximations of OðtÞ. To conclude, we numerically study scrambling in a long-range interacting model. This situation is particularly interesting as the decay of the OTOCs is not necessarily controlled by a butterfly velocity associated with the presence of a light cone [63,64].
Here, we consider a long-range XY model, realized, for instance, in trapped ion experiments [2], with timeevolution operator UðtÞ ¼ expð−iH LR tÞ, and with α a positive number controlling the range of the interactions. We show in Fig. 6 the space-time expansion of OðtÞ and O 3 ðtÞ for two different values of α. For α ¼ 1.5, the system satisfies Lieb-Robinson bounds [63,64], meaning that the system behaves effectively as if the interactions are local. This result is manifested by a light-cone spreading of the OTOCs. Conversely, for α ¼ 0.5, the characteristic decay time of OðtÞ is superlinear with respect to the operator position j.
In such a model with long-range interactions, the convergence of the modified OTOCs is slightly slower than for the examples shown before with nearest-neighbor interactions. However, the two types of light cones are well resolved by the modified OTOC O 3 ðtÞ. Within the light cone, O 3 ðtÞ oscillates weakly with time. These oscillations fade out with increasing n.

IV. ERRORS AND IMPERFECTIONS
We conclude our manuscript by presenting a study of errors and imperfections. From this analysis, we can draw the conclusion that OTOCs can be measured with good precision in AMO and superconducting qubit experiments with current technology, and with the total number of measurements compatible with state-of-the-art repetitions rates [10].

A. Statistical errors
Statistical errors arise in an experiment from a finite number N M of measurements per random unitary u to access the expectation values hWðtÞi u;k 0 and hV † WðtÞVi u;k 0 , and from a finite number of random unitaries N u used to estimate the correlation coefficients. This results in deviations E ¼ j½OðtÞ e − OðtÞj between estimated and exact correlation coefficients.
The scaling of statistical errors can be explained best in terms of the phenomenological model for scrambling Eqs. (9) and (10). Accordingly, using global unitaries, the typical value of hWðtÞi u;k 0 ∼ 1= ffiffiffiffiffi ffi 2 N p is suppressed exponentially. The number of projective measurements N M required to access OðtÞ up to a given error thus scales as 2 N . The protocol based on local unitaries accessing O n ðtÞ has a crucial advantage, provided that low resolution n ≪ N is sufficient to approximate OðtÞ. The typical value of the expectation values hWðtÞi u;k s ∼ 1= ffiffiffiffiffiffiffiffi ffi 2 LðtÞ p scales instead with the effective complexity 2 LðtÞ of the operator. Accordingly, the early-time dynamics of large chaotic systems subject to Lieb-Robinson bounds [22,23], but also the long time evolution of a MBL system, both of which correspond to small scrambling lengths LðtÞ ≪ N, are accessible with a moderate number of measurements, about 2 LðtÞ . These findings are confirmed by the numerical simulations of Fig. 7

B. Imperfections and decoherence
Our two protocols have a certain robustness against various types of experimental imperfections and decoherence. First, let us consider the effect of imperfect implementations of the local random unitaries u, which we can parametrize with a characteristic mismatch angle dθ in the realization of spin rotations. Our error analysis presented in Appendix F 1 shows that our estimations of the OTOCs are only affected in second order in dθ. The absence of a linear correction in dθ indeed illustrates a robustness mechanism inherent to protocols based on random measurements: OTOCs are estimated via an ensemble average over many random unitaries, which tends to "average out" the effect of errors from individual measurements. In particular, we expect the global protocol to have a similar robustness mechanism with respect to small errors in the preparation of global unitaries.
The two protocols are also robust against readout errors. In the case of decoherence, depolarization noise only rescales the values of the measurements of W while leaving statistical correlations unaltered [see Appendix F 2 and Fig. 7(c)]. For other sources of decoherence, the values of the correlations can be affected but in a way that can be clearly distinguished from unitary scrambling. While scrambling generically leads to a decay of statistical correlations, decoherence increases correlations. This behavior, which is the opposite of the case of protocols based on the reversal of time evolution [39,43,44], can be understood by noting that, for a Markovian dissipative evolution, the distance between two different states always decreases with time [65]. This case is illustrated in Fig. 7(d) for the estimation of the modified OTOC ½O 0 ðtÞ e .

V. CONCLUSION AND OUTLOOK
In the present work, we provide novel protocols to measure OTOCs for spin models, based on statistical correlations between measurement outcomes obtained from random initial states from both global and local random unitaries. These protocols can be implemented in state-ofthe-art quantum simulation experiments on various physical platforms, in particular, in Rydberg atoms, trapped ions, or superconducting qubits, which provide high repetition rates. The paradigm of extracting nonstandard correlation functions of quantum many-body systems from statistical correlations points to several interesting future developments, including extensions to measure modified OTOCs for Hubbard models [45].
While we have focused on second-order statistical correlations based on cross-correlating two random measurements, higher-order moments can also be used-at the price of higher requirements w.r.t. statistical errors-to access higher-order OTOCs [33] or OTOCs of arbitrary states ρ. In particular, we show in Appendix G 1 that the protocol can be adapted to access "symmetrized" OTOCs of arbitrary states, in particular, of thermal and ground states. This method could be used to extract the crucial temperature dependence of the Lyapunov exponents in models of high-energy physics such as the Sachdev-Ye-Kitaev models [27][28][29][30][31][32].

APPENDIX A: CORRESPONDENCE BETWEEN STATISTICAL CORRELATIONS AND OTOCs WITH LOCAL UNITARIES
In this section, we prove the relation between statistical correlations and OTOCs for the second protocol with local random unitaries u ¼ u 1 ⊗ Á Á Á ⊗ u N , where each u i is sampled independently from CUEðdÞ, d being the local Hilbert space dimension (d ¼ 2 for the spin-1=2 we consider). Here, we extend our diagrammatic approach by introducing a matrix-product-operator (MPO) [67] representation for the operators WðtÞ and V. This extension is shown in Fig. 8: Each (blue) physical index is contracted following the 2-design rule shown in panel (a), while the (green) bond links remain unchanged. To simplify the derivation of the proof, we rewrite O n ðtÞ as O n ðtÞ ≡ hWðtÞi u;n hV † WðtÞVi u;k 0 hWðtÞi u;n hWðtÞi u;k 0 ; ðA1Þ with hWðtÞi u;n ≡ P k s ∈E n ð−1=2Þ d½k s ;k 0 hWðtÞi u;k s ¼ Tr½r n WðtÞ, with the operator r n ¼ r 1;n ⊗ … ⊗ r 1;N written as a tensor product with r i;n ¼ jk The operator r n gathers all the information about the initial states jk s i ∈ E n and the corresponding chosen weights c k s ¼ ð−1=2Þ d½k 0 ;k s . Finally, we also use the tensor product decomposition We can now prove graphically, in Fig. 8, consists of all the 2 N−n permutations of the form This case leads directly to the desired equality In particular, for n ¼ N, we obtain O N ðtÞ ¼ OðtÞ.

APPENDIX B: OTOCs FOR HAAR SCRAMBLING
We now prove Eqs. (9) and (10) of the main text. The case j > LðtÞ is straightforward due to the commutativity of V and W. For j ≤ LðtÞ, the equality written for the OTOC OðtÞ follows directly from the 2-design identities [cf. Fig. 9(a)] and can also be found in Ref. [33].
For the modified OTOCs, we use the mapping to statistical correlations to calculate the values corresponding to Haar scrambling. We show in Figs. 9(b) and 9(c) the evaluation of hWðtÞi u;n¼0;1 hV † WðtÞVi u;k 0 , which can be adapted to derive hWðtÞi u;n hWðtÞi u;k 0 by replacing V by the identity matrix. Since u is a product of local random unitaries, we use the decomposition X ¼ X½L 1 ⊗; …, We also use the notation ρ 0 ½L 1 ¼ Vρ½L 1 V † . This method proves the desired identities To describe analytically the behavior of OTOCs in the MBL phase [24,26], we consider the l-bit model with time evolution UðtÞ ¼ expð−iHtÞ described by the Hamiltonian where h z i are random fields that, as we show below, have no effect on the OTOC; J ij ¼ J R;ij expð−jj − ij=ξÞ, where J R;ij is a random interaction amplitude with probability distribution fðJ R;ij Þ, which we assume to be uniform in ½−J z ; J z , and ξ is the localization length. We study the behavior of the OTOCs for the operators W ¼ σ x j and V ¼ σ x 1 . In the Anderson case J z ¼ 0, each l bit evolves independently, and we have OðtÞ ¼ O n ðtÞ ¼ 1.

Single disorder realization
We now address the general MBL case J z > 0 and first consider a single disorder realization of the interaction matrix ðJ ij Þ, writing the Heisenberg operators as with c i≠1 ¼ 1, c 1 ¼ −1. We then obtain as already shown in Refs. [24,26]. We now evaluate O n¼0;1 ðtÞ using the first line given in Eq. (A3). We first perform the trace operation with respect to the site j, withσ β i (β ¼ x, y, z) the set of Pauli matrices in the "copy" Hilbert space H i associated with site i. We can then calculate the trace over the remaining sites, for instance, All the factors that enter in the numerator and in the denominator in Eq. (A3) are identical, except for the position k ¼ 1 of the V operator. This case leads to

Averaged OTOCs
Considering a distribution of the random realizations of J ij , O n¼0;1 ðtÞ is now obtained from Eq. (A3), where the numerator and the denominator are replaced by their ensemble average. Repeating the above derivation, replacing each cosine contribution by the average we obtain Eq. (13).

Simulations of the disordered XXZ model
To describe the OTOCs in the MBL dynamics [24], we consider the XXZ model with h i z sampled from a uniform distribution ½−Δ; Δ.

APPENDIX D: NUMERICAL SIMULATION WITH DECOHERENCE
To study the competition between scrambling dynamics and decoherence, we consider the kicked Ising model. The dynamics is calculated from the Lindblad master equation with the initial condition ρð0Þ ¼ ρ 0 , and γ the spontaneous emission decay rate (identical for each qubit).

APPENDIX E: NUMERICAL STUDY OF STATISTICAL ERRORS
Here, we present complementary data to Fig. 7, showing the scaling of statistical errors in our protocol with global and local unitaries. The results are shown in Fig. 10. The data confirm that statistical errors for N M → ∞ decrease as 1= ffiffiffiffiffiffi N u p with a growing number of applied random unitaries N u , independently of LðtÞ and N [panels (a)-(c)]. Panel (d) shows the scaling of statistical errors with respect to N M for the first modified OTOC.

APPENDIX F: ROBUSTNESS OF THE PROTOCOLS
In this section, we give three examples showing how random measurements are robust against different kinds of perturbations. In each case, we consider the protocol with local unitaries.

Limited reproducibility of generated random unitaries
First, we analyze the robustness of our protocol with respect to imperfections in the generation of random unitaries. Specifically, we assume that hWðtÞi u;k s is obtained from a random unitary matrix u ¼ u 1 ⊗ Á Á Á u N , with u i ∈ CUEð2Þ, while the second measurement hV † WðtÞVi u 0 ;k 0 is obtained from a slightly different unitary u 0 ¼ u 0 1 ⊗ Á Á Á u 0 N , which we write as Here, R iγ denotes single qubit rotations for spin i along the γ axis, and θ in are assumed to be random angles drawn uniformly in ½−θ; θ, which represent the unwanted mismatch between the unitaries u and u 0 .
a. Analytical understanding of the robustness of the protocol To show that our protocol is robust against such imperfections, we calculate the estimator of the OTOC, obtained from To first order θ in , we can write where the second term vanishes because of θ in ¼ 0. This implies thatÕ n ðtÞ ¼ O n ðtÞ. Our protocol is thus robust to first order in θ against imperfections of the generated unitaries. As shown below, the quadratic contribution scales linearly with the characteristic size LðtÞ of the operator WðtÞ.

b. Numerical example
We now consider a numerical example for the model of scrambling presented in Sec. III B and the estimation of the first modified OTOC n ¼ 0. To assess the robustness of statistical correlations, we consider j > LðtÞ, so that the exact OTOC is O 0 ðtÞ ¼ 1, and imperfect generated unitaries as written above.
The numerical results are shown in Fig. 11 and confirm our analytical treatment: There is no linear correction in θ in the error E of the estimated OTOC. However, the error scales as θ 2 L, showing that the quadratic contribution does not vanish.

Robustness against the depolarizing channel
In the presence of depolarizing noise [39], the realized quantum state can be written as with 0 < p < 1 the depolarization probability, I the identity matrix, and U tot ðtÞ the desired time-evolution unitary operator. This case translates to the measurement outcomes [for the sequence shown in Fig. 1 with U tot ðtÞ ¼ UðtÞu (a) and U tot ðtÞ ¼ UðtÞVu where we use the fact that WðtÞ is traceless. Hence, the two measurements are simply rescaled, and the statistical correlations are not affected by the depolarizing channel. This case is illustrated in Fig. 7

Robustness against readout errors
We now show that our protocols are also robust against readout errors. Here, we consider, for concreteness, W ¼ σ z j . For each random unitary u, and initial state k s (k 0 ), we assume that each measurement of the operator is estimated as hWðtÞi est u;k s ¼ 2P est ðt; ↑; ujk s iÞ − 1; hVWðtÞVi est u;k 0 ¼ 2P est ðt; ↑; uVjk 0 iÞ − 1; ðF5Þ where P est ðt; ↑; jkiÞ is the estimated probability to detect spin j in state ↑ after time evolution from the state jki. We now consider that the estimated probabilities P est ðt; ↑; jkiÞ are built from a sequence of N M measurements, with X i the measurement outcomes obtained with error probability x: ProbðX i ¼ 1Þ ¼ ð1 − xÞPðt; ↑; jkiÞ þ xPðt; ↓; jkiÞ ¼ ð1 − 2xÞPðt; ↑; jkiÞ þ x: In the limit of an infinite number of measurements (N M → ∞), we obtain P est ðt; ↑; jkiÞ ¼ ProbðX i ¼ 1Þ, and thus hWðtÞi est u;k s ¼ ð1 − 2xÞhWðtÞi u;k s ; hVWðtÞVi est u;k 0 ¼ ð1 − 2xÞhVWðtÞVi est u;k 0 : Accordingly, as in the case of the depolarizing channel, readout errors simply rescale the value of the observables; i.e., the estimation of O n ðtÞ is not affected by readout errors.

APPENDIX G: ACCESSING OTOCs FOR NONINFINITE-TEMPERATURE STATES
In this section, we show how to extend our protocol beyond the case of infinite temperature ρ ¼ I=N H .

Thermal OTOCs from global unitaries
First, we show how our protocol can be extended to include finite-temperature corrections to OðtÞ. To this end, we employ a high-temperature expansion of the thermal density matrix ρ β ¼ expð−βHÞ=Z, with H being the many-body Hamiltonian of the system of interest Z ¼ Tr½expð−βHÞ and β the inverse temperature. Specifically, instead of considering the canonical finite-temperature OTOC with ρ β a thermal density matrix, we consider a symmetrized variant, introduced in Ref. [31], Here, and in the following, we assume WðtÞ to be Hermitian, V to be Hermitian and unitary, and all operators W, V, H to be traceless (which is the case for the spin models considered in the main text). We employ a hightemperature expansion of ρ β ∼ ðI − βHÞ=N H þ Oðβ 2 Þ. To first order in β, we find FIG. 11. Robustness of the protocol with respect to a mismatch between generated random unitaries in the estimation of the first modified OTOC. The error E is plotted as a function of the scaling parameter θ 2 L, for L ¼ 4, 6, 8 and 0 < θ < 0.2 Note that for all data points, the error is below 6%.
We now introduce the statistical correlations CðtÞ ¼ hWðtÞi u;k 0 hVWðtÞVi u;k 0 hHi u;k 0 ; with u global random unitaries of the CUE and hAi u;k 0 ¼ hψ u jAjψ u i. Using the 3-design properties [52] of the CUE, we obtaiñ Trðτ½WðtÞ ⊗ VWðtÞV ⊗ HÞ; which shows that O S ½ρ β ðtÞ can be accessed by separately measuring (i) the T ¼ ∞ value OðtÞ as described in the main text and (ii) the additional correlationsCðtÞ. This result shows that thermal OTOCs are experimentally accessible by measuring statistical correlations. The additional requirement compared to the measurement of OðtÞ is the measurement of the energy hHi u;k 0 of random initial states.

OTOCs for arbitrary states from global unitaries
We now generalize the previous approach and show that we can access, for an arbitrary state ρ (thermal or not), another symmetric variant of OTOCs, Note that for thermal states ρ ¼ ρ β , O S 0 ½ρ β ðtÞ and O S ½ρ β ðtÞ coincide to first order in β. Also, when V and ρ commute, this variant is equivalent to the standard OTOC, i.e., O S 0 ½ρðtÞ ¼ O½ρðtÞ. Now, we introduce the statistical correlations C 0 ðtÞ ¼ hWðtÞi u;k 0 hVWðtÞVi u;k 0 hρi u;k 0 ; with the same notations as in Appendix G 1. Note that the last term hρi u;k 0 simply consists in measuring the probability that the state ρ is in the random state jψ u i. We then obtain, using 3-design properties, Trðτ½WðtÞ ⊗ VWðtÞV ⊗ ρÞ; which corresponds to a sum of O S 0 ½ρðtÞ and of lower-order terms that can be measured independently (such as the infinite temperature OTOC). Thus, O S 0 ½ρðtÞ can be measured via statistical correlations.