Probing scrambling and operator size distributions using random mixed states and local measurements

The dynamical spreading of quantum information through a many-body system, typically called scrambling, is a complex process that has proven to be essential to describe many properties of out-of-equilibrium quantum systems. Scrambling can, in principle, be fully characterized via the use of out-of-time-ordered correlation functions, which are notoriously hard to access experimentally. In this work, we put forward an alternative toolbox of measurement protocols to experimentally probe scrambling by accessing properties of the operator size probability distribution, which tracks the size of the support of observables in a many-body system over time. Our measurement protocols require the preparation of separable mixed states together with local operations and measurements, and combine the tools of randomized operations, a modern development of near-term quantum algorithms, with the use of mixed states, a standard tool in NMR experiments. We demonstrate how to efficiently probe the probability-generating function of the operator distribution and discuss the challenges associated with obtaining the moments of the operator distribution. We further show that manipulating the initial state of the protocol allows us to directly obtain the individual elements of the distribution for small system sizes.


I. INTRODUCTION
Understanding how quantum information spreads across the degrees of freedom of a quantum system is a key part of developing a comprehensive picture of nonequilibrium quantum many-body physics.In this context, the notion of scrambling has attracted much attention over the past years due to its relevance in the study of closedsystem thermalization [1], quantum chaos [2], information retrieval in black holes [3,4], and quantum algorithms [5].Scrambling refers to the dynamical delocalization of quantum information [6] and can be diagnosed by the generation of entangled states from initially separable ones, or from the growth of initially local operators [7].
An approach to characterizing scrambling from the unitary evolution of an operator W (t) in a many-body system is to analyze the dynamics of so-called operator size distributions {P k (t)} [8][9][10].These can be obtained by coarse-graining the expansion of W (t) in a complete operator basis, the choice of which depends on the nature of the system.In the case of systems of N spin− 1 2 particles, a natural operator basis is the set of multi-body Pauli operators P N = {1, σ x , σ y , σ z } ⊗N / √ 2 N , which has dimension D = 4 N and forms an orthonormal basis.In the Heisenberg picture an operator W may at time t be * blocher@unm.edu† karthik.chinni@polymtl.ca‡ somanakuttan@unm.edu§ pablo.poggi@strath.ac.uk written as where U(t) is the unitary time evolution operator from the initial time t = 0 to time t and Λ j ∈ P N .Given this expansion of W (t), the operator size distribution is constructed by grouping the elements of the exponentially large Pauli basis according to their size s(Λ).Here the operator size corresponds to the number of non-identity operators in the Pauli string (i.e., its Hamming weight) and thus 1 ≤ s(Λ) ≤ N such that P N = ∪ N k=1 C k , where C k = {Λ | s(Λ) = k}.The resulting operator size distribution reads and measures the size of the support of W (t). It is easy to see that N k=1 P k (t) = 1 for all times t, hence the operator size distribution can be regarded as a coarse-grained probability distribution in the expansion coefficients of W (t).
Of particular interest is the case where W (0) is a sizeone (i.e.single-body) Pauli operator such that P k (0) = δ k,1 .As the operator grows and information becomes scrambled, the distribution shifts to higher values of k and grows in variance.The dynamics of P k (t) have been studied for various many-body models [8,10], and it holds a close connection to the Krylov picture of operator growth [7,11].
While the scrambling dynamics in the system may be described using the operator size distribution, it is not  2), out-of-time-ordered correlation functions (OTOCs), and the NOTOC approach proposed in this article using local expectation values over random mixed states.OTOCs allow (indirect) access of the moments of the operator distribution, while our proposal allows to probe the probability generating function F (x) or the individual probabilities directly, depending on the choice of initial state ρ0.(b) Illustration of the operator evolution (above the dashed line) and the proposed measurement protocol (below the dashed line).The operator size distribution P k (t) is a property of an operator in the Heisenberg picture (and thus independent of a choice of state).Our protocol proposes to access this property by measuring expectation values ⟨W (t)⟩ = Tr[ρ(t)W ] after the time evolution of a suitably-chosen set of initial states ρ0.
obvious how to access this distribution experimentally.A complete characterization of P k (t) for all k requires an exponential amount of resources; this is clearly seen if one reconstructs P k (t) from two-point correlation functions to obtain each of the coefficients f [Λ; W (t)] = Tr (ΛW (t)) in Eq. (1).A partial workaround is given by considering out-of-time-ordered correlators (OTOCs) which take the form It has been shown that moments of P k (t) can be computed by constructing averages of OTOCs over appropriate sets of operators {R i } [8,10,12,13].However, accessing even a single OTOC is often challenging in experiments due to the out-of-time-ordered nature of Eq. (3), and typically requires the use of many-body time-reversal operations [14][15][16][17][18] or auxiliary systems [19].In cases where OTOCs can be accessed without these tools (see for instance the method in [20]), reconstructing the moments of P k (t) requires the measurement of a large number of OTOCs and the task becomes unfeasible for high-order moments [10].
In this article we propose an alternative set of tools to experimentally probe the operator size distribution which circumvents the use of OTOCs completely.We combine the use of ensembles of initially separable mixed states with local random operations and local measurements at the final time in order to access a quantity G(t) which is explicitly not an OTOC, hence we name this quantity a "NOTOC".We show that our measurement protocol, depending on the choice of initial state ensemble, probes either the probability generating function (PGF) F (x, t) = k P k (t)x k of the operator size distribution (method A), or its elements {P k (t)} directly (method B), as illustrated in Fig. 1.We demonstrate our methods numerically and show that the operator size distribution can be accurately probed even when accounting for statistical noise stemming from averaging over random operations and from experimental imperfections.In our analysis of method A, we discuss inherent shortcomings of the problem of inverting the PGF to obtain P k (t) such as its sensitivity to noise, and we show for the case of the Ising model that the PGF can itself be seen as a good indicator of the presence of scrambling in the system.For method B we discuss the efficiency of the method as the system size N increases and show that individual probabilities P k (t) can be reliably obtained as long as k ≪ N .
Our proposed measurement protocol connects to previous works which focused on experimental schemes to diagnose scrambling.In particular, method A recovers a procedure put forward in Ref. [12] to measure operator growth in the special case when the initial states are pure.Likewise, the use of randomized operations to access properties of the operator size distribution makes this proposal complementary to the one in Ref. [20], where similar tools were used to measure OTOCs instead.Finally, the NOTOC proposed here can be seen as a generalization of the fidelity OTOC [1,14,15], which has been widely studied due to its ease of accessibility in experiments [21].Our analysis reveals that this quantity is inherently connected to the operator size distribution, although in a way that is fundamentally different from regular OTOCs, and that had not been revealed up to now to the best of our knowledge.
The structure of this article is as follows: In Sec.II we present the main tools to be used to probe the operator size distribution P k (t), including the choice of initial states, randomized operations and local measurements.Sections III and IV discuss two methods of implementing the necessary state preparation for our measurement protocol, with the method of Sec.III using partially polarized qubit states to obtain the probability-generating function for the operator size distribution, and the method of Sec.IV tailoring separable states to obtain the probability distribution elements directly.In Sec.V we compare our measurement protocol to previous results and proposals in the literature, extend the NOTOC measurement protocol to collective spin systems, and provide a general discussion of the NOTOC toolbox.Finally, we conclude on our work in Sec.VI and discuss possible future extensions.

II. NOTOC MEASUREMENT PROTOCOL
Our measurement protocol draws inspiration from the so-called fidelity OTOCs, which are OTOCs on the form Eq. ( 3) where one chooses the operator R = ρ 0 in Eq. ( 3) to be the projector onto the pure initial state ρ 0 [1,14,15].Fidelity OTOCs are experimentally accessible quantities as [21] ⟨W hence they require only the measurement of a single-time expectation value.However, this identity holds only for pure initial states, and the fidelity OTOCs for mixed initial states takes on a more elaborate form, as discussed in Ref. [21].
Inspired by the simple form of the fidelity OTOC in Eq. ( 4) as a squared expectation value of an operator W for a pure state, we define the quantity of interest for our protocol as with ρ 0 being a generic mixed state.The quantity G(t) defined in Eq. ( 5) is equivalent to an OTOC only for pure states [22], and to emphasize that G(t) is related to fidelity OTOCs while itself being a time-ordered correlation function at all times and for all initial states, we will refer to G(t) as a "NOTOC" throughout this article.In the remainder of this section we will demonstrate the utility of defining G(t) on the form Eq. ( 5), and in Sec.V we will comment further on its relation to fidelity OTOCs.Using Eq. ( 1) the NOTOC G(t) may be written on the form We now consider a generic mixed state on the form for which ⟨Λ i ⟩ = r i .We recall from Sec.I that C k = {Λ | s(Λ) = k} is the set of Pauli operators of size k.If we were to engineer an ensemble of initial states {ρ 0 } on the form Eq. ( 7) such that the state coefficients {r i } were independent, identically distributed random variables with vanishing mean r i = 0, the second term of Eq. ( 6) would vanish under averaging over this ensemble.Furthermore, if we require the state coefficients of the engineered initial state to have finite variance r 2 i = ∆ k for Λ i ∈ C k , the averaged quantity then reads Equation (10) reveals that G(t) is a linear combination of the elements of the probability distribution {P k (t)}, with coefficients proportional to the variance ∆ k times the 2-norm Tr[W † W ] = ∥W ∥ 2 2 of the operator W .The probability distribution {P k (t)} may be extracted from Eq. ( 10) using several methods, and we present two different measurement protocols for systems of spin-1/2 particles in the following sections III and IV.Our measurement protocols provide experimental access to the averaged NOTOC G(t) and the probability distribution {P k } through engineering of the initial state ρ 0 and subsequent measurement of the expectation value ⟨W (t)⟩ = Tr[ρ(t)W ] at the final time t.The main challenge thus lies in the preparation of random initial states ρ 0 whose coefficients {r i } must have appropriate statistics.
The operator size distribution is an operator property independent of the initial state of the system, hence the Heisenberg picture lends itself nicely to the analysis of the operator size distribution's evolution in time.However, in the NOTOC measurement protocol outlined above we probe the operator size distribution using expectation values ⟨W (t)⟩ ρ0 ≡ Tr[ρ 0 W (t)], which depend on the choice of initial state ρ 0 .In this way, the initial state is a control knob used by this protocol to access properties of the operator.As expectation values are quantities independent of the choice of picture, it is more natural to describe our proposed experimental measurement protocol in the Schrödinger picture as the time evolution of an initial state ρ 0 → ρ(t), with which we first evaluate the expectation value ⟨W (t)⟩ ρ0 ≡ Tr[ρ(t)W ], subsequently calculate the NOTOC Eq. ( 5), and finally recover the averaged NOTOC Eq. ( 10) by appropriate averaging over initial states.The NOTOC measurement protocol is thus illustrated in this way in Fig. 1(b), where the Heisenberg operator evolution is shown above the dashed line, and our proposed measurement protocol is illustrated below the dashed line.

III. METHOD A: ACCESSING THE PROBABILITY GENERATING FUNCTION OF THE PROBABILITY DISTRIBUTION
In this section we present an experimentally relevant measurement protocol for obtaining the squared expec-tation value Eq. ( 10) using mixed states similar to those used in NMR together with random local operations [23].The starting point of our measurement protocol is the preparation of the product state where each qubit is in a statistical mixture of being maximally mixed 1 and polarized along the z-axis, with the parameter ε controlling the amount of polarization.In Eq. ( 13) we have defined the subset C z k ⊂ C k of sizek Pauli operators that consist of only σ z terms (e.g., 1 ⊗ σ z ⊗ 1).The polarization parameter ε takes values |ε| ≤ 1, and we note that for ε = ±1 the initial state ρ ini is a pure state.
To create a random state ρ 0 for the experiment whose expansion coefficients r i in Eq. ( 7) satisfy the appropriate statistics for Eq. ( 10), we apply random local rotations to the initial state Eq. ( 13) via the unitary where rot is a random rotation of the ith qubit that we will discuss momentarily.In the last equality of Eq. ( 16) we expanded U rot Λ∈C z k Λ U † rot on the Pauli operators Q ∈ C k , as the random local rotations do not change the operator size.Comparing the form of Eq. ( 16) to that of Eq. ( 7), we make the identification We choose the random local rotations U rot such that the coefficients {q Q } are independent, identically distributed random variables with vanishing mean q Q = 0 and finite variance , which is consistent with the assumptions made in Sec.II.Substituting this back into Eq.( 10) -and choosing W to be a non-identity observable with trace Tr[W 2 ] = d [24] -we thus find averaged squared expectation value Equation ( 17) is an (at most) N th degree polynomial in ε 2 with coefficients proportional to elements P k (t) of the probability distribution of interest.This form is reminiscent of a probability-generating function (PGF) of the probability distribution {P k (t)} [25], and we now show that Eq. ( 17) is indeed a PGF by introducing the explicit form of ∆ k .
To obtain the correct statistics for the coefficients r Q , as well as to ensure that all operators Q ∈ C k are sampled for all k, we propose to take the single-qubit rotation operator U (i) rot to be sampled from a uniform distribution over SU (2).In each random instance the local rotation transforms the ith site Pauli-Z as with α = x, y, z, with The polar angle θ i and azimuthal angle ϕ i are thus random variables taking values in the interval [0, π) and [0, 2π), respectively.
The coefficients q Q of Eq. ( 16) are then expressible as products of these random numbers n α , with all factors being independent of each other thanks to the local rotations being uncorrelated.Taking each pair of angles (θ i , ϕ i ) to be uniformly distributed over the sphere, one readily obtains that n α = 0 leading to r Q = 0 as required by our protocol.Due to symmetry the variance ∆ k is expected to be independent of α = x, y, z, hence we can compute it for any component.We find that and thus we get a factor of a 1/3 for each non-identity operator in a given multi-body Pauli operator Q.This leads to the variance ∆ k = 1/3 k which in turn implies that Using the result Eq. ( 21) and letting x := ϵ 2 /3 for notational convenience, Eq. ( 17) may be rewritten as which is the probability-generating function (PGF) for the probability distribution {P k (t)} [25].From the PGF one may extract information about the corresponding probability distribution, including the elements and moments of the probability distribution.We point out that the uniform sampling of the continuous group SU(2) is not strictly necessary, as it suffices to sample over a finite set of rotations given the correct first and second moments.This is equivalent to constructing a y .Solid black lines are the exact PGF, while colored dots denote the PGF obtained using a subset of 1000 random rotations.Results are shown for several choices of the state parameter ε (see legend in panel b), and for both the integrable case θ = 0 (a) and the chaotic case θ = π/6 (b).For the case ε = 1 the blue crosses show the PGF obtained using a subset of 100 random rotations.
unitary 2-design and sampling the local operations from it, and can be done by choosing U rot such that each qubit takes on one of the values ±X, ±Y , and ±Z in each shot, with a total of 6 N unique rotation unitaries U rot needed to obtain Eq. ( 17) without approximation.This shortcut to uniform averaging also results in ∆ k = 1/3 k and r i = 0, and thus Eq. ( 22) is unchanged when using this method of averaging over a discrete set of rotations.

A. Application of the measurement protocol
In the remainder of this section we demonstrate via numerical simulations our measurement protocol for the case of the 1D tilted field Ising model.We present results on how to obtain the PGF approximately using a subset of random rotations, and show that the PGF may be used as an indicator of quantum information scrambling.Finally, we discuss an experimentally relevant method of extracting the elements of the probability distribution from the PGF and discuss its sensitivity to noise.
The Hamiltonian for the 1D tilted field Ising model is given by where the operators σ i α are the usual Pauli operators on site i with α = x, y, z.The model describes N spin-1 2 particles interacting in one dimension via nearest-neighbor interactions in the presence of an external magnetic field with both a transverse and longitudinal component which are parameterized by the angle 0 ≤ θ ≤ π/2.This model is a well-known platform for studies of many-body quantum chaos [10,26] since it is nonintegrable for generic choices of θ and presents two integrable limits at θ = 0 (where the model reduces to the usual transverse field Ising model) and θ = π/2 (where the Hamiltonian is diagonal in the z-basis).The scrambling properties of the dynamics generated by the Hamiltonian Eq. ( 23) have also been studied in relation to their quantum chaos characteristics.For instance, it has been established that even the integrable limit θ = 0 can lead to scrambling, and in this case the mean operator size typically presents long-lived oscillations for finite system sizes.The system typically shows the highest degree of chaoticity for θ ≃ π/6, where the operator sizes grow quickly and then equilibrate, and temporal fluctuations are suppressed [10].Throughout the following we will consider the case J = B with N = 6 qubits.
In Fig. 2(a) and (b) the solid lines behind the colored dots show the exact PGF F (x, t) for the edge site operator W (0) = σ (1) y ≡ σ y ⊗ 1 ⊗ . . .⊗ 1, for several choices of the state parameter ε.We consider the integrable case θ = 0 (a) and the chaotic case θ = π/6 (b) of the tilted field Ising model Eq.(23).For all t the two limiting cases F (x = 1, t) = 1 and F (x = 0, t) = 0, with x = ε 2 /3, follow from the definition of the PGF Eq. ( 22).We observe in Fig. 2 that the primary effect of varying the state parameter ε in the considered parameter range is a change of amplitude of the PGF.
The colored dots (blue crosses) in the two panels of Fig. 2 show the PGF extracted using a random subset of 1000 (100) rotations out of the 6 N = 46656 rotations needed for the exact result.Figure 2 demonstrates that we may obtain the PGF to very good accuracy using a heavily reduced number of rotations compared to the exact result, both in the integrable case θ = 0 and the chaotic case θ = π/6.This is encouraging for the experimental feasibility of implementing the present protocol.
In Ref. [10] moments of the operator size probability distribution were used as signatures of quantum information scrambling in the tilted field Ising model, with the integrable case θ = 0 leading to oscillatory dynamics of the mean operator size whereas in the chaotic case θ = π/6 the mean operator size grew and saturated only after initial oscillations.When comparing the curves for θ = 0 with the corresponding curves in θ = π/6 in Fig. 2, we see a clear difference in behavior for Jt ≥ 4, with θ = 0 curves at later times exhibiting oscillations that are not present for θ = π/6.In the following we thus explore whether the PGF may be used as an indicator of scrambling, similar to the mean operator size used in Ref. [10].
Figure 3 illustrates the PGF F (x, t) for different B-field angles θ and times t as a function of the PGF parameter x.For the initial time t = 0 the PGF is linear in x independent of the choice of θ, however at finite times t the behavior of the PGF is significantly different across different θs.In the chaotic case θ = π/6 the PGF F (x, t) approximately coincides with the result obtained from a Haar random state (dashed red line) for all displayed times greater than Jt = 2.52, indicating a fast equilibration.For θ = π/3 the solid curves trend toward the Haar random curve for increasing time Jt, but have yet to equilibrate at Jt = 10 (light blue curve).The two integrable cases θ = 0 and θ = π/2 display oscillatory behavior in the PGF F (x, t) for a given x, similar to what was observed in Fig. 2 -e.g., the dark blue Jt = 2.52 curve is below both the black Jt = 0 curve and the light blue Jt = 10 curve for all x in both cases.We also do not observe an equilibration of the PGF to the result for the Haar random state.
While the PGF F (x, t) primarily serves as a quantity from which one may extract information about the corresponding probability distribution, we now illustrate how the PGF itself may be used to characterize quantum information scrambling.We propose to do this by analyzing the time-average of the PGF for a fixed argument x, and the time-averaged temporal fluctuations In Ref. [10] we studied analog constructions for the mean operator size (i.e., the first moment of the operator size distribution) and found that they allowed to distinguish different scrambling and quantum chaos regimes of this model.In Fig. 4(a) we show the time-average F (x) as a function of the magnetic field angle θ.For different accessible values of x we see that that F (x) dips in the highly chaotic regime and grows near the integrable limits.This behavior originates in the fact that the chaotic case shows quick scrambling and subsequent equilibration to the Haar-random behavior, where the PGF is closer to 0, while the integrable cases show the largest typical values of the PGF, as was noted for θ = 0 in the discussion of Fig. 2. A similar functional form is observed for the timeaveraged temporal fluctuations ∆F (x) 2 which we show in Fig. 4(b).This indicates that temporal fluctuations of the PGF are suppressed in the chaotic and enhanced in the integrable cases, a behavior also observed also for the mean operator size in Fig. 3 of Ref. [10].Our present findings thus show that it is not necessary to extract the probability distribution {P k (t)} from the PGF F (x, t) in order to characterize the scrambling of quantum information for the considered Ising model.Rather, the time average and time fluctuations of the PGF are themselves good indicators of the presence of scrambling.As demonstrated previously in Fig. 2 the PGF can be approximated to good accuracy using a heavily reduced number of random rotations, and as the PGF is not itself sensitive to the noise unlike quantities extracted from the PGF (which we will comment on in the following Sec.III B), the PGF provides an experimentally relevant quantity for characterizing quantum information scrambling.The results of Figs. 3 and 4 further emphasize the utility of the PGF itself as a quantifier of quantum information scrambling.

B. Extracting probability distribution elements and moments from the PGF
In the absence of noise on the PGF F (x, t) -e.g., in a numerical study of the exact PGF -the corresponding probability distribution {P k (t)} can be extracted from the PGF using Eq. ( 22) in the following way.One first calculates the PGF F (x, t) for N different values of x, and then inverts the system of equations Eq. ( 22) to obtain the elements P k (t) of the probability distribution.However, in the presence of a noisy PGF due to, e.g., a finite number of sampled rotations or experimental imperfections, we find that the procedure of inverting Eq. ( 22) leads to high sensitivity to noise, which we associate with a poorlyconditioned linear system of equations.
For a more systematic approach we can use well-known properties of the PGF to recover the full probability distribution.From the definition of the PGF Eq. ( 22), one finds that [25] i.e., the elements {P k (t)} of the probability distribution are accessible through the evaluation of derivatives of the PGF with respect to the parameter x at x = 0.
Likewise the moments of the probability distribution may be accessed through derivatives of the PGF at x = 1: In the following we will focus on the extraction of the elements of the probability distribution via Eq.( 26).We leave the discussion of the moment extraction for later in this section.
The nth derivative of the PGF F (x, t) at x = 0 may be implemented for both numerical studies and in experiments using a forward-only finite difference method on the form [27] where a is the accuracy of the finite difference method.We leave the details of our implementation of the forwardonly finite difference method to appendix A. In practice, applying this method implies preparing several states of the form Eq. ( 13) with various levels of purity (controlled by the parameter ε).To explore the experimental feasibility of extracting the probability distribution elements from the PGF, in the following we characterize the sensitivity of the method to noise.The noise may have origin in an approximative PGF due to choosing a reduced subset of rotations, as discussed in Sec.III A, or stem from experimental imperfections.We simulate the effect of the noise in the following way.First the exact PGF F (x, t) is calculated, on top of which we add noise.The noise is assumed to be Gaussian and multiplicative in nature [28], hence we write the noisy PGF as where δ(η) ∼ N(0, η) is a normal-distributed number with vanishing mean and variance η 2 .From the noisy PGF we then use Eq. ( 28) to extract the elements {P η k (t)} of the probability distribution in the presence of noise.
We now consider the N = 6 transverse field Ising model with J = B and θ = 0.In Fig. 5 we show for the operator W (0) = σ (1) y the time averaged error of the extracted probability distribution elements as a function of the finite difference step size ∆x.The error is calculated by comparing the elements P η k (t) extracted from the noisy PGF with the exact elements P 0 k (t), and has been averaged over 100 realizations of the Gaussian distributed noise δ(η).An accuracy of a = 1 was used for the finite difference method; see details in appendix A.
When using finite-difference methods, there are two major sources of error: truncation errors associated with the accuracy a of the method, and rounding errors due to uncertainty on the input function (in the present case the rounding error is due to the noise added to the PGF) [27].The former error prefers ∆x as small as possible, while the latter prefers larger ∆x.Hence, in choosing the step size ∆x, one has to balance the contributions from the two errors.We observe for the case η = 10 −2 in Fig. 5 that there exists step sizes ∆x for which the first two elements P 1 and P 2 can be extracted with reasonable errors -for P 1 the time-averaged error is much smaller than unity for all choices of distance ∆x.For higher k the errors exceed 10 −1 for all choices of ∆x, and is dominated by rounding errors [27].We note here that the step size ∆x is upper bounded by (1/3)/(n + a − 1), as we need to sample (n + a − 1) points between x = 0 and x = 1/3.
By decreasing the noise amplitude by two orders of magnitude to η = 10 −4 , we observe in the lower panel of Fig. 5 that choosing ∆x ≥ 0.03 yields time averaged errors smaller than 10 −1 for the k ≤ 3 elements, while the k ≥ 4 elements of the probability distribution still carry too large  30)] for the elements of the probability distribution {P k }, averaged over 0 ≤ Jt ≤ 10.We have used multiplicative noise with strength η = 10 −2 (upper panel) and η = 10 −4 (lower panel), added to the PGF F (x, t) before using finite difference derivatives with step size ∆x and accuracy a = 2 to extract the noisy probability distribution elements {P η k (t)}.In both panels the curves from top to bottom at ∆x = 0.02 are the time averaged errors ∆P η k for P6 (black dashed), P5 (blue dashed), P4 (red dashed), P3 (black solid), P2 (blue solid), and P1 (red solid), averaged over 100 noise realizations.
errors to be usable for the analysis of scrambling.We thus see that the above derivative-based method Eq. ( 26) of extracting the elements of the probability distribution becomes increasingly sensitive to noise with the degree of the derivative one needs to evaluate, and even for highly precise measurements of the PGF F (x, t) the method Eq. ( 26) of extracting the probability distribution from the PGF is too sensitive to noise to allow us access to P k for k ≥ 4.
As a consequence of the high sensitivity to noisewhether the noise originates in experimental imperfections or an approximate PGF due to a finite number of sampled rotations -it does not seem experimentally feasible to extract the probability distribution elements from the PGF.Instead, one should use a protocol dedicated to the extraction of the probability distribution elements.To this end we present in Sec.IV an experimentallyrelevant protocol that is able to access the elements of the probability distribution without using the PGF.
Finally, we comment on the extraction of the moments of the probability distribution using Eq.(27).Accessing the moments of the probability distribution using the previously described finite difference method requires one to create states with x > 1/3 to approximate the derivative.However, with the constraint of |ε| ≤ 1, it is impossible to access x > 1/3 through the procedure outlined in this section.If one could engineer an ensemble of correlated initial states, the variance ∆ k could possibly be manipulated sufficiently to allow access to values of x > 1/3.However, in the present work we do not investigate this further and leave it for possible future extensions of this work.

IV. METHOD B: ACCESSING ELEMENTS OF THE PROBABILITY DISTRIBUTION FOR SMALL OPERATOR SIZES
In the previous Sec.III we presented a measurement protocol that can access the operator size probability distribution {P k (t)} by measuring the probability-generating function (PGF) F (x, t) in Eq. ( 22).Here we return to the general result in Eq. (10) and develop an alternate NOTOC measurement protocol to access the individual probabilities.This is achieved by employing a different choice of initial state that replaces the one in Eq. ( 13).The protocol will provide a direct way of obtaining the probabilities and bypasses the sensitivity issues encountered when trying to invert the PGF.
In order to directly access P k for a given 1 ≤ k ≤ N , we choose a subset M k of k particles in the system and prepare them the classically correlated state while the other N − k particles are left in a maximally mixed state.The resulting state is clearly unentangled as it can be produced as a statistical mixture of product states in the computational basis.Crucially, it has the desired property of being written solely in terms of weightk operators, as discussed in Sec.II.Thus we may proceed in a similar fashion to method A: we apply a round of uniformly-chosen random rotations on each of the particles in M k .The resulting full state takes the form where the modified set C k (M k ) is composed of all kbody Pauli operators acting on the particles in M k .The application of random, local rotations ensures that the coefficients q Q are independent, identically distributed random variables with mean q Q = 0 and variance q 2 Q = 1/3 k .It follows from Eq. ( 9) that and thus we find In order to obtain P k (t) exactly one needs to obtain G(t) for all possible subsets of M k of k particles, thus requiring N k repetitions of the protocol.While this is in general inefficient, it can be feasible as long as k and N are not too large.In practice, the number of initial states can be reduced in at least two ways.One option is to approximate the sum in Eq. ( 34) by randomly sampling a reduced number of sets M k .Alternatively one can exploit symmetries of the system.For instance, the Hamiltonian in Eq. ( 23) has a reflection symmetry with respect to the middle of the chain which can be easily leveraged to reduce the state count by a factor of 2. If the system has a full translational symmetry, then the state count can be reduced by a factor of N .The procedure to achieve this is described in the Appendix B.
In the following we present numerical results that demonstrate this method's ability to reconstruct P k (t) for various instances of the tilted field Ising model introduced in Sec.III A with N = 6 particles.We choose the initial operator to be W (0) = σ (2) x .Results are obtained by using all possible choices of M k for each k, and by sampling over M rot realizations of randomized rotations.
In Fig. 6 it can be seen that the dynamics of each size probability P k (t) is faithfully reproduced by the present method, in some cases using as little as M rot = 100 rotations.For larger operator sizes k ≥ 4 the effects of finite sampling are more pronounced, however this is to be expected as the typical probabilities are also smaller.The cases displayed in the figure correspond to θ = 0 (integrable transverse field model, blue curves and symbols) and θ = π/6 (chaotic model, red curves and symbols), and it can be readily seen that the protocol accesses the typical features expected in both cases, namely a fast spread of the operators (as seen in the decay of P 1 (t) at short times), a subsequent oscillatory behavior for θ = 0, and equilibration for θ = π/6.
We note that the method presented in this section will not be feasible to obtain the full probability distribution for large system sizes N ≫ 1, since at some point an exponential number of the subsets M k 's would be required to exactly recover them from Eq. (34).However, this method can be used to verify whether a system a system has scrambled to k−body operators as a long as k ∼ O( 1).An alternative benchmark is to compare the values of the obtained P k 's with the ones corresponding to Haarrandom evolution, i.e.
We show these values as gray dashed lines in each of the plots of Fig. 6, where we observe that the observable tends to equilibrate to these values in the chaotic regime (θ = π/6) of the model.

A. Extensions to other systems
In this section we present an extension of the NOTOC measurement protocol introduced in Sec.II to collective spin systems, where every spin interacts with every other spin in the system.This class of systems preserves is the collective spin operator.Some well known Hamiltonians of this type include p-spin models (p = 2 case is referred to as the Lipkin-Meshkov-Glick model) and the quantum kicked top model [29][30][31].The dynamics of these models is often studied in the symmetric subspace, the subspace of the whole Hilbert space associated with J = N/2 quantum number where N is the number of spins.The dimension of this subspace increases linearly with the number of spins in the system, It is natural to study scrambling in this type of system by decomposing the time-evolved operators in a basis associated with the polynomials of collective spin operators, referred to as the spherical tensor operator basis [10,32].A spherical tensor operator T (k) q is an operator that transforms under rotations in the same manner as spherical harmonics [33] where D(α, β, γ) = e −iJzα e −iJyβ e −iJzγ and therefore D q ′ q (α, β, γ) ≡ ⟨k, q ′ |D(α, β, γ)|k, q⟩.The explicit form of the spherical tensor operators is given by [32] where C J m ′ J m;k q = ⟨J, m ′ |J, m; k, q⟩ is a Clebsch-Gordan coefficient, and q = {−k, −k + 1, ..., k} for a given rank of the tensor operator k = {0, 1, ..., N }.These operators form an orthonormal basis Tr (T spans the Hilbert space associated with the symmetric subspace.Note that a kth rank tensor operator is simply a kth order polynomial of collective spin operators, hence the basis consists of polynomials of collective spin operators ranging from order 0 to order N .The rank of these operators can be used to construct the operator size distribution for this class of systems, similar to the operator size (Hamming weight) used for the Pauli basis.
It is then natural to analyze the dynamics of the system by considering an operator of a particular rank k at the initial time, and characterizing the scrambling dynamics of the system by considering for each tensor rank k the time evolution of the element P k (t) of the operator size distribution.
In the following we describe the details of a NOTOC measurement protocol for collective spin systems that can access the elements P k (t) of the operator size distribution for small operator sizes.As in Sec.II we require the initial state to be prepared from an ensemble that satisfies r i = 0 and r 2 i = ∆ k in Eq. ( 7).The nonzero coefficients r i are here associated with a particular rank k in the expansion of the density operator in the spherical tensor operator basis {T (k) q }.To obtain states that satisfy these properties for the expansion coefficients, we start with a mixed state given by where k > 0. e 0 is the absolute value of the smallest eigenvalue of T (k) 0 , which is used to ensure the semidefiniteness of the density operator.Note that T (k) 0 is a kth order polynomial in J z .For instance, T 0 and c (2) 0 are normalization coefficients.For k = 1 the state ρ 0 is a thermal equilibrium state of a system placed in an external magnetic field at high temperatures, ρ = 1 d (1 + ϵJ z ) where ϵ ≪ 1 [34].Higher-order states could potentially be prepared as equilibrium states of a system with Hamiltonian consisting of higher-order J z terms.Applying global rotations of the form D(ϕ, θ) = D(α = ϕ, β = θ, γ = 0) = e −iJzϕ e −iJyθ on the state ρ 0 leads to where r k,q ′ = (e 0 d ss ) −1 ⟨k, q ′ |D(ϕ, θ)|k, 0⟩, and the angles θ and ϕ are sampled randomly from a uniform distribution on the surface of a sphere.This state has zero mean and nonzero variance [35] as expected, where dΩ = sin θdθdϕ.Equation ( 41) is analogous to Eq. (21) in that the expression for the variance depends only on a coarse-grained property of the operator basis element, namely, the operator size for the Pauli basis and the operator rank for the spherical tensor basis.Hence we have where W is the operator of interest.This result shows that method B of Sec.IV can be naturally adapted to collective spin systems.In general, the procedure shown in this section illustrates how to choose a combination of initial states and randomized operations tailored to the choice of operator basis such that the average NOTOC connects to operator size distributions.Notice that N + 1 repetitions of the above protocol for different k will provide us probabilities associated with all operator sizes.However, the state in Eq. ( 38) is not easily accessible for higher values of k, so this protocol might only be suitable for small system sizes.

B. Relation to previous proposals
The toolbox presented in Sec.II presents some noteworthy connections with previous works which have studied how to diagnose complex many-body dynamics in different settings.For instance, Qi et al. [12] proposed a method to probe the growth of an operator O in quantum quench experiments using pure product states of qudits (of local dimension d L ) followed by random local operations.The authors showed that the variance δO(t) 2 of the expectation value ⟨O(t)⟩ over the random realizations yields where F (x, t) is the probability-generating function associated with the operator size distribution of O(t).For qubits, d L = 2 and so the method probes the PGF at x = 1/3.This is exactly the case for the NOTOC when choosing states of the form proposed in method A (Sec.III) using pure states (ε = 1), see Eq. (22).Therefore, our proposed method A recovers the protocol of Ref. [12] for the case of pure states and generalizes it by showing that the PGF can be sampled in a continuum of values x ≤ 1/3 by using mixed states of qubits.It remains to be studied whether this generalization carries over to the case of qudits with d L > 2. Additionally, these methods do not allow to access x > 1/3 directly and it is unclear whether the toolbox of Sec.II provides a way around this by using a clever choice of initial states.Regarding this aspect, we point out that a recent work proposes an alternative method to access the PGF F (x, t) (in principle for any x) using a single-particle mixed states (similar in form to Eq. ( 31) when k = 1) and resorting to forward U (t) and backward U † (t) evolution.Importantly, the analysis we presented in Sec.III B concerning the large sensitivity to noise of the process of obtaining the operator distribution {P k (t)} from its PGF F (x, t) applies to all methods that aim at obtaining the PGF.Our findings indicate that obtaining the distribution from the PGF might be unfeasible in experiments, but also show that properties of the PGF itself could be used a probe for scrambling directly.More detailed work should be carried out to explore this further.The idea of using mixed states to probe properties of operator evolution was also used recently by Peng et al. in Ref. [23], where the goal was to measure single-site two-time correlation functions on the form i Tr[σ ] in an NMR experiment.This measurement was carried out by first preparing the weakly polarized initial state ρ 0 ∝ (1 + ε i σ (i) z ), then allowing this initial state to become locally randomized by the effect of on-site disorder, and ultimately measuring a tunable observable using inductive measurements, rotations, and on-site disorder.An analogous method for measuring two-site two-time correlation functions was also proposed by the authors.The measurement protocol proposed in Sec.III of the present work shares significant overlap in methodology with that of Ref. [23], however the goals of the two measurement protocols differ and thus the prepared random mixed states and ultimate measurements are also different.For the purpose of extracting the operator size distribution, we note that by using the two-time correlation function studied in Ref. [23] one will have to extend the method of Ref. [23] to all m-site correlation functions, with m ≤ N .This will likely not be feasible due to the non-local nature of the operators to be measured, as well as the issue of the exponentially growing number of operators one needs to measure which was discussed in Sec.IV for our proposed measurement protocol.
Finally, we comment on the connection between our proposal and that of Vermersch et al. in [20], where the authors propose a way to measure OTOCs without using time-reversal operations or auxiliary systems.Instead, their proposal is based on performing randomized unitaries on a set of initial states and extracting the OTOCs from the statistical correlations between the measurement results.In principle, this method allows one to reconstruct the operator size distribution if one repeats the procedure for (exponentially many) choices of the operator R in Eq. ( 3).This can be achieved by using averages of OTOCs to obtain the moments of the {P k (t)} distribution, as outlined in Sect.VI of [10].In contrast, our method can be seen as a way of using similar tools (i.e.preparation of product states, randomized local operations, and local measurements) to probe the operator size distribution directly, thus bypassing the calculation of OTOCs.

C. Relation to fidelity OTOCs
Finally, we discuss the physical interpretation of the different tools and quantities used to study quantum information scrambling.The first aspect is related to the fidelity OTOCs, introduced in Sec.II, which are a class of correlation functions obtained from the usual OTOCs of the form in Eq. ( 3) by choosing the early-time operator R to be the projector onto the initial state R = ρ 0 .The use of fidelity OTOCs attracted widespread attention in the community because they are technically an OTOC but can be measured as a single expectation value of a (often) local operator W when the initial state is pure, ρ 0 = |ψ 0 ⟩ ⟨ψ 0 | [21].Fidelity OTOCs have interesting connections to quantities like the quantum Fisher information [1] and the Loschmidt echo [36].
The fact that typically ρ 0 is a non-local operator makes the fidelity OTOC relinquish the usual interpretation of OTOCs as measures of information scrambling.In particular, the relation between OTOCs involving Pauli operators and moments of the operator size distributions [6,10] does not apply to fidelity OTOCs.However, our present work shows that fidelity OTOCs are indeed connected fundamentally to operator size distributions if one takes the initial pure state |ψ 0 ⟩ to be a product state and then considers the average fidelity OTOC over many realizations of local random rotations on the initial state.

VI. CONCLUSION
In this article we have proposed a measurement protocol for probing quantum information scrambling by measuring operator size distributions.Our measurement protocol requires the preparation of separable mixed states followed by local operations and final-time measurements of local operators, and circumvents the typical use of outof-time-ordered correlation functions (OTOCs) to probe scrambling properties [37,38].We have demonstrated that the choice of initial separable mixed states in our measurement protocol provides multiple ways to access the operator size distribution, and we comment on the experimental feasibility for two particular methods, one based on first extracting the probability-generating function for the operator size distribution, and the other focused on obtaining the elements of the operator size distribution directly.
The application of our measurement protocol is illustrated in detail for the 1D tilted field Ising model, a well-known platform for studying many-body quantum chaos [26,39,40].We numerically demonstrate the characterization of quantum information scrambling for this model using our proposed measurement protocol.We have related our proposed protocol with well-established methods of characterizing scrambling, including those based on OTOCs, and were able to establish connections between the well-established fidelity OTOC and operator size distributions using our results presented in this article.Finally, we have exemplified the extension of our measurement protocol to other types of quantum systems by considering our protocol for the case of collective spin systems.The collective spin case emphasizes further the role that state preparation plays in our measurement protocol.
In the discussion of our proposed measurement protocol's connection to the probability-generating function, we found that the extraction of moments of the operator size distribution was not possible due to constraints on the prepared initial states that prevent us from accessing x > 1/3 in the probability-generating function F (x). From preliminary numerical results we expect that the extraction of the first and second moments of the operator size distribution would be both resilient to noise and reasonable to implement in experiment if one could obtain values of x > 1/3, hence finding a way to prepare initial states allowing the extraction of moments of the operator size distribution would be an obvious extension of the present work.This potential extension should be compared to the proposal in Ref. [9] where the probabilitygenerating function can be accessed for any x, but at the expense of requiring the implementation of time-reversal of the many-body evolution.
In Sec.VA we presented the extension of our NOTOC measurement protocol to collective spin systems.The extension of the measurement protocol to many-body qudit systems and other systems of interest for studying the nature of scrambling would be an interesting task that we leave for future work.We point out that studies of operator size distributions for many-body systems beyond qubits are also scarce, with some exceptions [6,41].In particular, we note the case where the system of interest interacts with the environment, thus forcing one to probe scrambling and the operator size distribution in the presence of decoherence [9].The NOTOC measurement protocol may be extended to open quantum systems using the analysis presented in Sec.II B of Ref. [21], from which a detailed analysis of the effect of decoherence may be carried out.Crucially one should revisit the definition of the coarse-grained operator size distribution Eq. ( 2) and consider the effect of decoherence on, e.g., the normalization of this distribution.Additionally, the relationship between the fidelity OTOC Eq. ( 4) and the NOTOC Eq. ( 5) becomes more complicated for open quantum systems [21].
While in theory one can increase the accuracy a arbitrarily to minimize the error O(∆x a ) of the approximation for any ∆x ≪ 1, in practice the method becomes increasingly sensitive to noise as we increase the accuracy a.This increased sensitivity to noise is due to the number of points sampled n + a limiting the values that ∆x can take, and in the vicinity of vanishing ∆x the dominant error is not the approximation of the derivative in Eq. (A1), but rather rounding errors due to the presence of the noise [27].
Figure 7 displays a comparison between the extraction of the elements of the probability distribution for accuracies a = 1 (solid curves) and a = 2 (dashed curves).Plotted is the time averaged error ∆P η k as defined in Eq. (30), averaged over times Jt ∈ [0, 10] for each realization of the noise Eq. ( 29) and subsequently averaged over 100 noise realizations.We are considering the same system parameters and operator as in Sec.III B. We observe for the higher noise amplitude η = 10 −2 (upper panel) that the lower accuracy method generally outperforms the higher accuracy method for the considered step sizes ∆x, the exception being for P 1 when considering ∆x ≥ 0.02.The crossover between the a = 1 and a = 2 curves at ∆x ≈ 0.02 for P 1 indicates a transition from being dominated by the sensitivity of the method to the added noise Eq. ( 29) for smaller ∆x values, to being dominated by the error due to the finite difference method's approximations for larger ∆x values.
For the case of η = 10 −4 (lower panel in Fig. 7) we clearly see that the a = 2 method outperforms the a = 1 method for P 1 , and a clear crossover point is also observed for P 2 .Therefore, when working with a sufficient small noise amplitude such as η = 10 −4 , it may be beneficial to use a higher accuracy for the finite difference method.that T † l U (t)T l = U (t).Then we have when we have defined Q l ≡ T † l QT l and used the translation invariance property.By definition, the LHS of the last equation equals and so we have that Suppose we start from a fixed initial state ρ 1 = , where the sum is over only a subset of operators of weight k, starting at position 1 up to k.If we measure ⟨W 1+l (t)⟩ for l = 0, . . ., N − 1, then In conclusion, that means that measuring ⟨W 1+l (t)⟩ yields the same result as having done the protocol starting from T l ρ 1 T † l , thus reducing the state count by a factor of N provided one can measure expectation values in all sites in the case of full translational invariance.The procedure works similarly if one only has a reflection symmetry (i.e.open boundary conditions).

Figure 1 .
Figure 1.(a) Relationship between the operator size distribution P k (t) of Eq. (2), out-of-time-ordered correlation functions (OTOCs), and the NOTOC approach proposed in this article using local expectation values over random mixed states.OTOCs allow (indirect) access of the moments of the operator distribution, while our proposal allows to probe the probability generating function F (x) or the individual probabilities directly, depending on the choice of initial state ρ0.(b) Illustration of the operator evolution (above the dashed line) and the proposed measurement protocol (below the dashed line).The operator size distribution P k (t) is a property of an operator in the Heisenberg picture (and thus independent of a choice of state).Our protocol proposes to access this property by measuring expectation values ⟨W (t)⟩ = Tr[ρ(t)W ] after the time evolution of a suitably-chosen set of initial states ρ0.

Figure 3 .
Figure 3. Exact PGF F (x, t) illustrated for several times t and B-field angles θ as a function of the parameter x.Solid lines are the PGF for N = 6 qubits and, for comparison, the red dashed line is the PGF for a Haar random probability distribution.

Figure 4 .
Figure 4. Time average F (x) and fluctuations ∆F (x) 2 of the PGF F (x, t) as a function of the B-field angle θ, plotted for three values of x corresponding to ε = 1 (red), ε = 2/3 (blue), and ε = 1/3 (black curve).The time average and fluctuations are calculated for 5 ≤ Jt ≤ 80 to exclude the initial transient dynamics from the results.

Figure 5 .
Figure 5.Time averaged error ∆P ηk [see Eq.(30)] for the elements of the probability distribution {P k }, averaged over 0 ≤ Jt ≤ 10.We have used multiplicative noise with strength η = 10 −2 (upper panel) and η = 10 −4 (lower panel), added to the PGF F (x, t) before using finite difference derivatives with step size ∆x and accuracy a = 2 to extract the noisy probability distribution elements {P η k (t)}.In both panels the curves from top to bottom at ∆x = 0.02 are the time averaged errors ∆P η k for P6 (black dashed), P5 (blue dashed), P4 (red dashed), P3 (black solid), P2 (blue solid), and P1 (red solid), averaged over 100 noise realizations.

Figure 6 .
Figure 6.Application of Method B to obtaining the operator size distribution {P k (t)} for the case of the Ising model of Eq. (23) with N = 6 and J/B = 1.Each panel shows a different value of k = 1, . . ., 6. Full lines correspond exact numerical results obtained by solving the Heisenberg evolution of the initial operator, here chosen to be W (0) = σ (2) x .Symbols correspond to numerical simulations of Method B (Sect.IV) with different choices of the number of sampled rotations Mrot = 100 (circles) and Mrot = 500 (crosses).Different colors denote different regimes of the model: integrable case θ = 0 (blue) and chaotic case (θ = π/6) (red).