Digital Quantum Simulation, Trotter Errors, and Quantum Chaos of the Kicked Top

This work aims at giving Trotter errors in digital quantum simulation (DQS) of collective spin systems an interpretation in terms of quantum chaos of the kicked top. In particular, for DQS of such systems, regular dynamics of the kicked top ensures convergence of the Trotterized time evolution, while chaos in the top, which sets in above a sharp threshold value of the Trotter step size, corresponds to the proliferation of Trotter errors. We show the possibility to analyze this phenomenology in a wide variety of experimental realizations of the kicked top, ranging from single atomic spins to trapped-ion quantum simulators which implement DQS of all-to-all interacting spin-1/2 systems. These platforms thus enable in-depth studies of Trotter errors and their relation to signatures of quantum chaos, including the growth of out-of-time-ordered correlators.


I. INTRODUCTION
Digital quantum simulation (DQS) is based on the "Trotterization" of quantum dynamics, i.e., on approximating unitary Hamiltonian evolution by a sequence of quantum gates [1][2][3][4][5][6][7][8][9][10][11]. It is a fundamental conceptual question under which conditions this approximation is controlled. A recent numerical study [12] by some of us found that Trotter errors in DQS of generic many-body systems remain bounded below and proliferate above a dynamical transition to many-body quantum chaos [13]. Motivated by these findings we revisit the kicked top, a paradigmatic model of single-body quantum chaos [14]. Resorting to this well-studied model system allows us to gain insights into Trotter errors in DQS of collective spin systems. Moreover, the kicked top connects smoothly to a paradigmatic DQS of an Ising chain with long-ranged power-law interactions.
The dynamics of the kicked top is described by the time-dependent Hamiltonian which combines precession of the top's spin S around a fixed axis, H x = h x S x , with non-linear "kicks" given by H z = J z S 2 z /(2S + 1), which are applied periodically at times t = nτ for all n ∈ Z. Here, S µ with µ = x, y, z are quantum angular momentum operators. The inverse spin size can be regarded as an effective Planck constant, eff = 1/S [14], and in the limit S → ∞, a semiclassical description of the dynamics applies. Then, the precession of a classical top is only slightly perturbed by weak * lukas.sieberer@uibk.ac.at kicks, whereas strong kicks cause the top to tumble chaotically [15]. This classical chaotic motion is reflected in the spectrum of the Floquet operator (we set the "true" Planck constant to unity, = 1), which determines the evolution of a quantum kicked top during a period of duration τ : In the chaotic regime of strong kicking, the spectral statistics of U τ is described by Dyson's ensemble of random orthogonal matrices [16]. Indeed, it is a defining feature of quantum chaotic systems that their spectral statistics are universal and obey predictions from random-matrix theory (RMT) [14]. Among the model systems of quantum chaos, the kicked top stands out due to its extraordinary faithfulness to RMT. The discovery of this property [15] initiated a surge of theoretical studies [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33]. Recent developments include proposals [34] to diagnose chaos in the kicked top by measuring out-of-time-ordered correlators, and the discovery of critical quasienergy states [35]. Experimentally, the kicked top was realized as the spin of single Caesium atoms [36] with S = 3 and as the collective spin of an ensemble of three superconducting qubits [37] corresponding to S = 3/2. Novel experimental possibilities include atomic species with high spin becoming available in labs such as Dysprosium [38] with S = 8 and Erbium [39] with S = 19/2; implementations as collective spin in quantum simulators of all-to-all coupled spin-1 /2 systems [40,41] in which S 50 can be realized; and condensates of ultracold bosonic atoms [42] corresponding to even larger collective spins S on the order of several hundreds.
Formally, the dynamics of the kicked top, described by repeated application of the Floquet operator U τ in Eq. (2), is equivalent to DQS of a system with Hamiltonian H = H x + H z : In DQS, time evolution is "Trotterized" [43,44], i.e., the run-time t of a simulation is split into n steps of duration t/n, and within each step the time evolution operator is approximately factorized, U (t) = e −iHt ≈ e −iHzt/n e −iHxt/n n = U n τ . (3) We note that Trotterization is not uniquely defined, and we discuss different schemes in Appendix A. To establish the formal equivalence with a periodically kicked system, in Eq. (3) we identified the Trotter step size with the kicking period, τ = t/n. One of the main goals of DQS is to enable studying the dynamics of quantum many-body systems in regimes where the growth of entanglement prohibits classical simulations. A typical model system which is implemented in trapped-ion quantum simulators [2,8,10,40,41,[45][46][47][48][49][50] is the Ising chain described by with spin operators S µ i where µ = x, y, z and i = 1, . . . , N , and the power-law coupling strength J z,α . For such long-range interacting systems, extensivity of the Hamiltonian is ensured by the Kac prescription [51], In the limit α → 0, this model describes all-to-all interactions between spins and becomes equivalent to the quantum kicked top. In the opposite limit of nearestneighbor interactions, α → ∞, and with the addition of a longitudinal field in Eq. (4), some of us addressed in Ref. [12] the crucial question of "Trotter errors," i.e., errors in DQS due to the approximate factorization e −iHτ ≈ e −iHzτ e −iHxτ . By performing extensive numerical simulations, Ref. [12] found that the Trotter error in observables such as the magnetization exhibits threshold behavior as the Trotter step size τ is changed. Below the threshold, the Trotter error admits a perturbative expansion in τ , while it is uncontrolled above the threshold. Reference [12] attributed this behavior to a transition from dynamical localization to a regime showing features of many-body quantum chaos. Studies of many-body quantum chaos [52][53][54][55][56] typically have to resort to numerics (see Ref. [57] for an exception). In single-body quantum chaos, on the other hand, sophisticated techniques have been developed to establish deep connections to RMT even analytically. To harness this knowledge, we first focus here on the all-to-all interacting limit α → 0 of the Hamiltonian Eq. (4). In this limit, the collective spin S 2 = S 2 x + S 2 y + S 2 z with components S µ = N i=1 S µ i where µ = x, y, z becomes a constant of motion. Accordingly, the 2 N -dimensional Hilbert space of a system of N spin-1 /2 can be decomposed into decoupled subspaces of fixed total spin S. In each of these subspaces, the Trotterized dynamics (3) reproduces the dynamics of a kicked top of sizes S. We focus specifically on the subspace with maximal spin S = N/2, in which the Hamiltonians given in Eq. (4) reduce to those in Eq. (1) (up to an inconsequential additive shift of the total energy and a rescaling J z → 2 (N + 1) J z /N to accommodate usual conventions). The dimension D = 2S + 1 = N + 1 of this subspace scales linearly with the number of spins N . Chaotic motion of the kicked top is ergodic within this subspace, and it does not explore the full many-body Hilbert space of dimension 2 N . In this regard, the kicked top does not display many-body quantum chaos.
We show that in this setting nevertheless observables exhibit the same threshold behavior as described in Ref. [12] for a generic many-body system. Moreover, we show that the threshold is determined by the onset of quantum chaos in the kicked top. While this implies that physical observables in DQS can remain close to their "ideal" values, we find that the fidelity of the simulated many-body quantum states drops sharply with time for any non-zero Trotter step in large systems. In view of the variety of possible experimental implementations of the kicked top listed above, we determine requirements on system sizes and decoherence rates to observe the threshold behavior. It is an interesting conceptual question whether in an ideal DQS Trotter errors could be controlled up to arbitrarily long times and in the thermodynamic limit in which N and hence S tend to infinity. For collective spin models, we find that this is indeed the case.
Finally, we connect the results for collective spin systems to the case of generic many-body systems [12]. In particular, we explore the influence of deformations of the kicked top to α > 0, where the problem can no longer be solved on the basis of collective spin operators. Instead, we resort to exact diagonalization of systems up to N = 14, where we recover the essential features of the kicked top. The threshold in Trotter errors persists over the entire studied range 0 ≤ α ≤ 3 from all-to-all to dipolar interactions.

II. TROTTER ERRORS IN DQS OF COLLECTIVE SPIN SYSTEMS
We consider two measures of Trotter errors: deviations of the expectation values of physical observables [12] and the fidelity of the quantum state [26,58] obtained in DQS. Interestingly, these measures exhibit strikingly different behaviors.
Here, |S, S z is an eigenstate of the spin z-component with eigenvalue S z . The magnetization error is defined as  The discontinuities in the curves pertaining to short times are due to the rounding t/τ to integer numbers of Trotter steps in the definition of the temporal averages. We chose S = 50 corresponding to recent experiments with trapped ions [40,41]. Smaller (larger) values of S result in a more rounded (sharper) threshold. (c,f) Absence of a sharp threshold in fidelity. In contrast to physical observables, the (c) time-averaged fidelity F(t) does not exhibit persistent threshold behavior as a function of the Trotter step size τ . Instead, (f) the fidelity drops with increasing system size as ∼ 1/ √ S even for values of τ within the regular regime. (d,e) Long-time averages of the (d) magnetization error ∆M and (e) heating as quantified by QE. Dotted lines represent results from perturbation theory as explained in the main text. The data for S → ∞ is obtained from the semiclassical equations of motion (C4). For (a-c) and (f), the initial state is |θ, φ = |0, 0 = |S, S , while in (d) and (e) the initial state is |θ, φ = |0.25π, 0 .
where the expectation value · τ is taken under Trotterized dynamics with step size τ while · pertains to time evolution with the target Hamiltonian H. The accumulated effect of Trotter errors is quantified by the stroboscopic temporal average ∆M (t) = 1 nt nt n=1 ∆M (nτ ) where n t = t/τ , which is shown in Fig. 1(a). For this figure and throughout the paper, we assume that for both µ = x and µ = z, with h x = 0.1J z , J x = 0.7J z , and h z = 0.3J z . However, our results do not depend on the precise values of these parameters (except for finetuned points of higher symmetry such as h x = h z , J x = J z ). As shown in the figure, after a characteristic time of J z t th ≈ 25 for the spin size S = 50 and the above parameters, the time-averaged magnetization error exhibits clear threshold behavior: While Trotter errors remain small for Trotter step sizes τ below a threshold value J z τ * ≈ 2.5, errors become significant for larger τ . A concrete experimental realization will also suffer from other imperfections. However, the regimes of bounded and uncontrolled Trotter errors can still be distinguished if t c > t th , where t c is the coherence time. This is discussed further in Sec. VI. As detailed below, τ * determines the radius of convergence of the Floquet-Magnus (FM) expansion for the effective Hamiltonian H τ defined by U τ = e −iHτ τ (see Ref. [32] for the FM expansion in the context of the kicked top). The FM expansion is guaranteed to converge always if the "driving frequency" Ω = 2π/τ exceeds the time-averaged operator norm of the Hamiltonian, i.e., Ω 1 τ τ 0 dt H(t) [59]. Intuitively, for such high driving frequencies, the system cannot absorb energy from the drive and thus the energy as measured by the expectation value of the time-averaged Hamiltonian H = H x + H z is conserved. The range of values of the Trotter step τ in which this criterion applies vanishes for large spin sizes S since H(t) ∼ S so that τ 1/S. However, we find in DQS of collective spin systems, that heating is suppressed up to a Trotter step size τ * even for large values of S. We quantify this absence of heating for τ < τ * by the "simulation accuracy" [12,55] where E τ (t) = H(t) τ , and E 0 = ψ 0 |H|ψ 0 is the energy of the initial state at t = 0. In Q E (t), the energy increase is normalized with respect to the difference between the energy of an infinite-temperature state, E T =∞ = tr(H)/D where D = 2S + 1 is the Hilbert space dimension, and E 0 . A value of Q E (t) = 1 thus indicates heating to infinite temperature. Figure 1(b) shows the finite-time average, Q E (t) = 1 nt nt n=1 Q E (nτ ), which clearly exhibits threshold behavior characterized by the same timescales t th and τ * as the magnetization error.
These results demonstrate that few-body observables are quite robust against Trotter errors for τ < τ * . This is in stark contrast to the accuracy of the full unitary time-evolution operator: The difference U (t)−U n τ , which quantifies the error made in the Trotterization in Eq. (3), grows at least linearly in both the simulation time t = nτ and the system size N [60]. Similarly, the robustness of observables does also not extend to the quantum state |ψ τ (t) = U nt τ |ψ 0 obtained under Trotterized dynamics in DQS. In Fig. 1(c), we show stroboscopic temporal averages of the fidelity F(t), i.e., the absolute value of the overlap of |ψ τ (t) with the ideal target state |ψ(t) = U (t) |ψ 0 , As illustrated in the figure, the temporal average F(t) approaches the ideal value of F(t) = 1 for τ → 0, but drops sharply already for J z τ 1, and in particular for Trotter step sizes which are much smaller than the threshold value τ * identified above for the magnetization error and the simulation accuracy. At this threshold value, i.e., for τ ≈ τ * , there is another noticeable drop in F(t). However, also in the region below τ * the fidelity vanishes with increasing system size as shown in Fig. 1(f). For large system sizes, the decay of the fidelity approaches F(t) ∼ 1/ √ D set by the Hilbert space dimension D = 2S + 1 [26]. Thus, there is no persistent threshold behavior in the fidelity of the simulated quantum state. In light of this fragility of the quantum state, the robustness of local observables seems rather remarkable.
Thus far, we have focused on the temporal evolution of Trotter errors and, in particular, on the dynamical buildup of a threshold in physical observables. Notably, in collective spin systems, both the infinite-time and thermodynamic limits are accessible. Figure 1(d) shows the infinite-time average of the magnetization error. The data for S < ∞ is obtained by exact diagonalization of the Floquet operator U τ . In the limit S → ∞, the effective Planck constant vanishes, eff = 1/S → 0 [14]. Therefore, in the thermodynamic limit, the spin components obey semiclassical stroboscopic evolution equations as detailed in Appendix C. These evolution equations can be iterated efficiently, and the long-time average for S → ∞ shown in Fig. 1 corresponds to n = 10 6 iterations. Fluctuations in the data for ∆M decrease upon further increasing the number of iterations.
As can be seen in Fig. 1(d), the threshold behavior identified above at finite times persists for t → ∞. We note that the different shapes of the curves for the magnetization error in Figs. 1(a) and (d) are due to the different choice of the initial state, which is |θ, φ = |0, 0 = |S, S and |θ, φ = |0.25π, 0 in (a) and (d), respectively. We comment on these state-dependent variations, also of the threshold value τ * itself, in Sec. III below.
For small Trotter steps, the time-averaged magnetiza-tion error admits a controlled expansion in powers of τ , (10) This expansion can be obtained analytically from the FM expansion for the effective Hamiltonian by employing time-dependent perturbation theory [12] which we generalize to arbitrary initial states in Appendix D. As illustrated in Fig. 1(d), the perturbative expansion up to second order in the Trotter step size captures the dependence of ∆M on τ up to the threshold τ * of uncontrolled Trotter errors. This indicates that the radius of convergence of the expansion Eq. (10) coincides with τ * . Figure 1(d) shows data for spin sizes ranging form S = 8 to S → ∞. The smallest value shown, S = 8, is realized in experiments with Dysprosium atoms [38], while S = 64 can be achieved in 1D [2,49] or 2D arrays of ions [41]. A threshold in ∆M is visible over the entire range of values of S. The value of τ * depends only very weakly on system size, and crucially it remains finite even for S → ∞. These observations do not depend qualitatively on the choice of the Hamiltonian parameters or the angles θ and φ specifying the initial coherent spin state.
In Fig. 1(e), we also show the long-time average Q E of the simulation accuracy Q E (t) defined in Eq. (8). Clearly, heating is suppressed to perturbatively small values, Q E = q 1 τ + q 2 τ 2 + O(τ 3 ), for Trotter steps up to the same threshold value τ * as for the magnetization error, indicating a finite radius of convergence of the FM expansion even for S → ∞. As in the expansion for the magnetization error in Eq. (10), the coefficients q 1 and q 2 can be obtained analytically by time-dependent perturbation theory, see Appendix D. On the other hand, for τ > τ * and with increasing spin size, the time-averaged simulation accuracy approaches the maximum value of Q E = 1. This is to be expected for chaotic dynamics, as we detail in the following and in Appendix B below.

III. INTERPRETATION OF THE PROLIFERATION OF TROTTER ERRORS AS QUANTUM CHAOS
As we show in the following, the threshold behavior in the magnetization error and the simulation accuracy can be traced back to the transition from dynamical localization to quantum chaos in the kicked top. Indeed, for strong kicking, which corresponds to large Trotter steps τ , the kicked top is well-known to exhibit quantum chaos [14]. Then, statistical properties of the eigenvectors and eigenvalues of the Floquet operator Eq. (2) are faithful to RMT predictions. In the case of Hamiltonians H x,z as given in Eq. (7) with all coefficients h x,z and J x,z different from zero, the Floquet operator does not have any symmetries other than time reversal [14]. Thence, the relevant RMT ensemble is the circular orthogonal ensemble (COE) of unitary and symmetric matrices.  Depending on the Trotter step size, trajectories which emanate from initial states marked by the colored dots are regular or chaotic. This is reflected in the simulation accuracy QE shown in (g). The simulation accuracy is averaged over the first 10 6 recursion steps.
Our first indicator of quantum chaos is connected to the inverse participation ratio (IPR) for a state |ψ 0 , where {|φ m } m denotes the eigenbasis of the Floquet operator U τ . We average the IPR over the basis of the target Hamiltonian {|ψ n } n and take its inverse to obtain the participation ratio (PR), In this definition, we rescaled the PR by a factor of 1/D so that it assumes values between two limiting cases: A low PR thus indicates that the two eigenbases are very similar (localization) whereas a high PR indicates equal absolute overlaps between all the eigenstates (delocalization). We show the PR for the kicked top in Fig. 2(a).
Similar to the Trotter errors discussed above, the PR exhibits a steep increase at a critical value of the Trotter step size τ . However, we would like to emphasize a crucial difference between the PR and Trotter errors of observables: The threshold value τ * for Trotter errors depends on the initial state of the time evolution as illustrated Figs. 1(a,b) and (d,e). On the other hand, the PR as defined in Eq. (12) contains an average over the entire eigenbasis of the target Hamiltonian, and consequently it is a "global" rather than a state-dependent and thus "local" measure of the onset of chaos. We discuss this distinction in more detail below. In Fig. 2(a), two sets of curves are displayed: Solid lines pertain to the PR calculated with eigenvectors |φ m of the Floquet operator U τ given in Eq. (2); dashed lines correspond to a symmetric Floquet operator, which realizes a higher-order Trotter decomposition as explained in Appendix A, and is related to the Floquet operator U τ in Eq. (2) by a unitary transformation, In the chaotic regime of large Trotter steps, the PR obtained from the Floquet operator U τ Eq. (2) assumes the value predicted for the CUE of RMT, see Eq. (B1). On the other hand, the symmetric Floquet operator U τ,s yields a value of the PR that matches the COE prediction obtained from Eq. (B2). This difference can be understood by noting that since the Floquet opreator U τ,s is symmetric, U τ,s = U T τ,s , it can be regarded as an element of the COE. The PR Eq. (12) probes the statistics of eigenvectors of the Floquet operator U τ,s , and in the chaotic regime it approaches the RMT prediction. As we explain in Appendix B, the RMT prediction is obtained by taking an average over the ensemble of eigenvectors of random matrices of the COE, i.e., unit-norm vectors with real components. U τ,s is related to the non-symmetric Floquet operator U τ by the unitary transformation given in Eq. (15). Since a unitary transformation does not affect the spectrum of an operator, both U τ and U τ,s exhibit the same level statistics in the chaotic regime which we discuss in more detail below. However, as shown in Fig. 2(a), the unitary transformation Q strongly affects the PR. Indeed, our results for the PR corresponding to U τ indicate that applying Q to the ensemble of eigenvectors corresponding to the COE yields an ensemble of vectors that resembles the one of the CUE, which is comprised of unit-norm vectors with complex components.
At small Trotter steps, the PR does not assume the lowest possible value ∼ 1/D = 1/(2S + 1) = 1/(N + 1). Instead, it appears to converge to a finite value, which indicates that even in this localized regime the basis states |ψ m of H are spread out over ∼ D Floquet states |φ m . This is similar to the many-body case considered in Ref. [12] and in Sec. V below, where for small τ an exponential number of Floquet states is required to represent generic states, but with an exponent that is smaller than the maximal possible one.
Signatures of the transition to chaos can also be found in the spectrum σ(U τ ) of the Floquet operator. The eigenphases {θ n } n of the unitary operator U τ are defined through σ(U τ ) = {e iθn } n and are shown in Fig. 2(b) for a system with S = 8. Already by visual inspection of the figure, three regimes can be distinguished: (i) In the region of very small τ , the eigenphases of the Floquet operator U τ coincide with the phases corresponding to the ideal time-evolution operator over one Trotter step, e −iHτ , which are shown as gray lines in the figure. (ii) In the second regime, the phases start to wind around the interval [−π, π]. This leads to level crossings, which, however, are avoided ones with only very narrow gaps. (iii) A dramatic change can be observed in the third regime, where phase repulsion becomes strongly pronounced, which indicates significant overlaps between the respective eigenstates.
A quantitative one-parameter measure of the degree of phase repulsion is given by the average adjacent phase spacing ratio r [61]. In terms of the phase spacings δ n = θ n+1 − θ n , the average spacing ration is defined through where D = 2S + 1 denotes the dimension of the Hilbert space. The phase spacing ratio takes universal values for distinct ensembles of random matrices: In the absence of correlations between the eigenphases, which corresponds to Poissonian statistics of the adjacent phase spacings, the phase spacing ratio is given by r POI = 2 ln(2) − 1 ≈ 0.39; In contrast, the repulsion between eigenphases of random matrices sampled from the COE, is reflected in a higher universal value of the adjacent phase spacing ratio, r COE = 0.5996(1) [62]. For the Floquet operator U τ in Eq. (2), the adjacent phase spacing ratio is shown in Fig. 2(c). Again, three regimes can be distinguished: In a region of size ∼ 1/S in which the spectral width of the target Hamiltonian H is smaller than Ω = 2π/τ , r is determined by the spectral statistics of H. The corresponding value of r depends on the spin size S and does not assume any of the universal ones given above. This is not surprising: It is well-known that quantum systems with time-independent Hamiltonians that have classical counterparts with only one degree of freedom exhibit non-universal level statistics [13,14]. This is the case, in particular, for the considered spin systems, where the expectation values of the collective spin operators span a two-dimensional phase space corresponding to a single degree of freedom in the classical limit. Indeed, for a classical top S with spin components S x,y,z , the dynamics is restricted to the surface of a sphere, S 2 = S 2 x + S 2 y + S 2 z , and can be parameterized in terms of the polar and azimuthal angles θ and φ, respectively, which span the twodimensional phase space (θ, φ) ∈ [0, π) × [0, 2π). In the second regime in Fig. 2(c), the eigenphases θ n of the Floquet operator wind around the unit circle. As shown in the previous section, for values of τ up to J z τ 2, the dynamics is well-captured by the time-independent effective Hamiltonian H τ obtained in the FM expansion. For the average phase spacing ratio we find r ≈ 0.26 in this regime. This non-universal value depends on the choice of parameters h x,z and J x,z . Finally, for large values of τ , the adjacent phase spacing ratio assumes the universal value predicted by RMT for the COE.
We note that while the behavior of observables considered in Sec. II is exactly analogous to the phenomenology of generic (short-range interacting) many-body systems [12], here we find a crucial difference: Even in the dynamically localized regime, the time-evolution operator of an interacting many-body system displays COE phase statistics [53]. Therefore, the statistics of the eigenphases of the Floquet operator does not allow one to distinguish between localized and many-body chaotic regimes.
The participation ratio and the adjacent gap ratio include averages over the entire set of eigenstates and eigenphases of the Floquet operator U τ in Eq. (2). Hence, they pertain to "global" properties of the dynamics, and show that at large Trotter step sizes the Floquet operator is faithful to RMT-a hallmark of single-body quantum chaos. In contrast, the magnetization error and the simulation accuracy discussed in the previous Sec. II measure Trotter errors in the time evolution starting from a single initial state. Therefore, they yield a more "local" measure of irregularity in the dynamics, and indeed we found that the value τ * of the threshold in Trotter errors depends on the initial state. To explain this initial-state dependence, we consider again the semiclassical limit of the top. From the stroboscopic Heisenberg equations for the spin-components of the top, we obtain in the semiclassical limit the Poincaré sections shown in Fig. 2(d-f). Regimes of regular and chaotic dynamics coexist in the classical phase space, and their relative weight is tuned by the Trotter step size τ . This coexistence is reflected in the initial-state dependence of τ * in the thermodynamic limit as can be seen in Fig. 2(g), where we show the simulation accuracy for trajectories which emanate from the initial states marked by colored dots in Fig. 2(d-f). Note that for our choice of parameters h x,z and J x,z , states pertaining to θ ∈ {0, π} are particularly stable whereas states around θ ≈ π 2 are generally more sensitive. The signatures of chaos discussed so far show clearly that the proliferation of Trotter errors is a concomitant effect of the transition from regular to chaotic dynamics in the kicked top. However, in contrast to Trotter errors, in experiments these signatures are not directly accessible. We proceed to discuss how chaos can be quantified experimentally by measuring out-of-time-ordered correlation functions (OTOCs).

IV. SIGNATURES OF DYNAMICAL LOCALIZATION AND QUANTUM CHAOS IN OUT-OF-TIME-ORDERED CORRELATORS
A hallmark of chaotic dynamics in classical systems is the butterfly effect, which is caused by extreme sensitivity to initial conditions: In a chaotic system, two trajectories that emanate from nearby points in phase space deviate from each other exponentially fast. The rate entering this exponential, which is called the Lyapunov exponent, is a quantitative measure of chaos.
A generalization of these concepts to the quantum domain can be given in terms of OTOCs [63][64][65]. The "infinite-temperature" OTOC for two initially commuting operators V and W is defined as where D is the Hilbert space dimension, and the Heisenberg operator W (t) = U (t) † W U (t) is evolved unitarily with U (t). Here, we consider a single spin S with Trotterized dynamics, for which the Hilbert space dimension is D = 2S + 1 and we evaluate C(t) stroboscopically at multiples of the Trotter step such that t = nτ and W (nτ ) = U n † τ W U n τ . In certain highly chaotic systems [64][65][66], the OTOC exhibits an extended period of exponential growth, and the growth rate serves as a quantum analog of the classical Lyapunov exponent. Moreover, the infinite-time limiting value of the OTOC is maximal if the time-evolution operator U (t) is faithful to an ensemble of RMT [67]. Due to its properties as a diagnostic of chaos, the OTOC has gained a lot of attention recently in a wide range of physical contexts including highenergy, condensed matter, and AMO physics [34,41,[64][65][66][68][69][70][71][72][73][74][75][76].
As we show in the following, in the kicked top, both the short-time growth as well as the long-time saturation of the OTOC are indicative of the transition from dynamical localization to chaos. This is illustrated in Fig. 3, which shows the infinite-temperature OTOC (17) with unitary V = e −iϕSz and Hermitian W = S z -a particular choice of operators that facilitates the measurement of C(t) [34,41,[77][78][79][80][81][82] as also discussed in Appendix E. However, we note that for the value ϕ = 10 −4 shown in Fig. 3, the OTOC reduces to the form The growth of the OTOC is depicted in Fig. 3(a). As expected in a chaotic system, for large Trotter step sizes the OTOC grows exponentially, C(t) ∼ e λJzt [73,83]. This exponential growth extends up to the Ehrenfest time t E ∼ |ln( eff )| [84,85], where the effective Planck constant is determined by the spin size, eff = 1/S [14]. For times t beyond the Ehrenfest time t E , the OTOC saturates at a value C(t) = C COE that is compatible with replacing the time-evolution operator in Eq. (17) by a random unitary matrix drawn from the COE as shown in Appendix B.
The behavior of the OTOC is markedly different in the regular regime. In particular, for the smallest nonzero Trotter step size J z τ = 0.73 shown in Fig. 3(a), the OTOC is essentially indistinguishable from evolution under U (t) = e −iHt with the target Hamiltonian H = H x + H z , which is labeled as J z τ = 0 in the figure. From a short-time expansion of the time evolution operator as e −iHt = 1 − iHt + O(t 2 ) it follows immediately that at very short times the OTOC grows algebraically, C(t) ∼ t 2 . The algebraic behavior persists approximately up to a time J z t of order one. Then, the OTOC grows exponentially also in the regular regime. Similar behavior, i.e., exponential growth of the OTOC in a (classically) regular parameter regime, was recently discussed in Ref. [74] for the kicked rotor. However, there are two key differences between the regular and chaotic regimes: (i) In the regular regime, the exponential growth gives way to much slower growth at late times, and eventually saturates at a value below C COE as shown in Fig. 3(b); (ii) Moreover, the growth rate λ, which we extract from the numerical data in the region of exponential growth, assumes a constant but weakly system-size dependent value in the regular regime as shown in Fig. 3(c). In contrast, it rises sharply upon increasing the Trotter step size τ beyond the threshold to chaotic dynamics.
In Fig. 3, we show data for spin sizes S = 256, 512, 1024. These rather large system sizes allow us to determine the growth rate λ very accurately due to the large Ehrenfest time t E ∼ ln(S). However, we found (see Appendix E 2) that the two regimes in the dependence of λ on the Trotter step size-λ ≈ const. in the regular regime and strong increase of λ with τ in the chaotic regime-can be clearly distinguished also in systems with S 50 as in recent measurements of OTOCs with trapped ions [41].

V. THE KICKED TOP AS A LIMIT OF LONG-RANGE INTERACTING SPIN CHAINS
In the previous sections, we discussed extensively the kicked top as an ideal model system for the study of fundamental questions of both DQS and the associated quantum chaos threshold. Now we explore to which extent our observations persist under natural deformations of the Hamiltonian in experimentally relevant settings. For this purpose, we consider Ising chains with algebraically decaying instead of all-to-all interactions, which exhibit a direct realization in trapped ion quantum computers [2,8,10,40,41,[45][46][47][48][49][50]. Specifically, we study a kicked spin chain as defined in Eqs. (4) and (5). By tuning the exponent α, the long-range Ising model allows us to smoothly interpolate between the two paradigmatic model systems of the infinite-range kicked top studied above and the nearest-neighbor Ising chain, reached at α → 0 and α → ∞, respectively. Motivated by realistic trapped-ion experiments, we consider open boundary conditions. For the numerical data shown in what follows, we choose h x /J z = 1/4, but the main picture does not change for different parameter sets. We use a Lanczos algorithm with full reorthogonalization to propagate the initial state |ψ 0 = N i=1 |↑ over up to 2·10 4 periods. For the quantities shown below we perform a temporal average over all periods.
As explained in Sec. I, the many-body Hilbert space of N spin-1 /2 can be decomposed into subspaces with fixed total spin S. The initial state |ψ 0 given above belongs to the subspace of maximal spin S = N/2, and for α → 0 the dynamics is restricted to this subspace. That is, for α → 0 the Trotterized dynamics of the spin chain as defined in Eqs. (4) and (5) becomes equivalent to the kicked top Eq. (1) with S = N/2 and a rescaling of the parameter J z as pointed out in Sec. I. On the other hand, for non-zero values of α, the dynamics explores the 2 Ndimensional many-body Hilbert space.
In Fig. 4(a), we plot the numerically obtained data for the simulation accuracy Q E defined in Eq. (8) in the long-time limit as a function of Trotter step τ and interaction exponent α for a fixed system size N = 14. To clearly illustrate how the kicked top is approached in the limit α → 0, Fig. 4(d) shows cuts through the same data for four values of the exponent α form α = 0 to α = 0.38. For small τ , we find that the properties observed already with the kicked top extend to α > 0. Specifically, Q E is small with Q E = q 2 τ 2 + O(τ 3 ) (the first-order term vanishes for the initial state |ψ 0 = N i=1 |↑ ). Upon increasing the Trotter step for α 1/4, however, Q E does not reach the value Q E → 1 expected in the quantum chaotic regime. We attribute this to the small number N = 14 of simulated spins and consequently to large finite-size corrections. This agrees also with our observations in Fig. 1(e) for the kicked top at α = 0, where significant finite-size effects in Q E appear (note also the different choice of initial state and Hamiltonian parameters in Fig. 1(e)). Instead, for α 1/4, we find a clear crossover from a perturbative regime at small Trotter steps τ , implying controllable Trotter errors, to a manybody quantum chaotic region with Q E → 1. Except for the range 0 ≤ α 1/4, the influence of the interaction exponent α on the simulation accuracy is only minor in that also the crossover region only shifts slightly upon varying α.
A similar picture emerges from the IPR defined in Eq. (11), which measures the localization properties of the state |ψ 0 in the eigenbasis {|φ m } m of the Floquet operator. For α > 0, the accessible Hilbert space grows exponentially with the number of spins, D = 2 N . To account for this exponential dependence, we introduce the rate function λ IPR = ln(IPR)/N [12]. In Fig. 4(b), we  show the ratio λ IPR /λ D , where λ D = ln(2) (1 − 3/N ). This value incorporates leading-order finite-size corrections as follows: On the one hand the system obeys both an inversion as well as Ising Z 2 symmetry, such that D = 2 N −2 ; On the other hand, RMT allows us to estimate the IPR in an ergodic phase as IPR ∼ 2/D for D → ∞, see Appendix B. Combining these two results we obtain the value of λ D given above. We note that here we restrict ourselves to comparing the IPR obtained for the initial state |ψ 0 to the RMT prediction for the CUE Eq. (B1). As discussed in Sec. II, we expect that a symmetrized Floquet operator yields an IPR that is consistent with the COE result (B2).
The rate function λ IPR essentially reproduces the observations from the simulation accuracy in the range 0 < α ≤ 3 with only slight modifications, which we again attribute to finite-size corrections. The case α close to α = 0, however, deserves a more detailed discussion. While a system at α = 0 with total spin S = N/2 can access only D = 2S + 1 = N + 1 quantum states in the Hilbert space, with any infinitesimal deviation from α = 0 this number grows to D = 2 N . For a system with finite N as we consider here, however, we expect that the effective behavior of D has to interpolate from the polynomial dependence at α = 0 to the exponential dependence in N upon increasing α. This opens up a crossover window at small α, whose size diminishes for larger N and where the rate function λ IPR exhibits strong finitesize corrections, as is also visible in Fig. 4(b). Within the crossover window, the threshold to chaos is nevertheless clearly visible if instead of the rate function λ IPR one considers the IPR Eq. (11) itself. This is illustrated in Fig. 4(e), which is based on the same data as panel (b) and shows cuts for several small values of α. Beyond the crossover window, the rate function shown in Fig. 4(b) reproduces the observations from the simulation accuracy, while the crossover region from quantum localized to chaotic appears slightly broader.
The previous considerations provide numerical evidence for a quantum many-body chaos transition in the Trotterized time evolution. For completeness, we now aim to confirm our general arguments on controllable Trotter errors by studying also the long-time limit of a local observable, namely the magnetization M as done also for the kicked top. Our numerical data for the time-averaged magnetization error Eq. (6) is shown in Fig. 4(c). Again, we observe two different regimes depending on the Trotter step τ . Compared to the simulation accuracy and the IPR in Fig. 4(b), the crossover region appears much sharper, which is also reflected in the cuts through the data shown in Fig. 4(f). The overall behavior remains similar: Delocalization of the timeevolved state in Hilbert space implies for ∆M a vanishing value as we find across a wide region in the α-τ plane.
For small values of τ instead, ∆M is nonzero and Trotter errors remain controllable.
Overall, this analysis suggests that a continuous deformation of the kicked top Hamiltonian to algebraically decaying interactions does not influence our general observations in the previous parts of the manuscript, while only details such as the location of the quantum manybody chaos threshold can exhibit slight modifications.

VI. DISCUSSION
In this work, we connected the well-studied quantum chaos of the kicked top with Trotter errors in DQS of collective spin systems. Our approach of harnessing the available knowledge on the kicked top, which is a paradigmatic model of single-body quantum chaos, is complementary to the numerical study of the connection between Trotter errors in DQS of generic manybody systems and many-body quantum chaos presented in Ref. [12].
We obtained the following key results: (i) Trotter errors of observables exhibit the same behavior as was found numerically in "true" many-body systems in Ref. [12], with a regime of controlled Trotter errors at small Trotter steps separated by a sharp threshold from a regime of large errors. Hence, completely analogous physical phenomena can occur in a system with an exponentially large Hilbert space dimension of 2 N for N spin-1 /2 and in a collective system of effective dimension 2S + 1 = N + 1. (ii) There is no analogous threshold behavior in the fidelity of the quantum states in DQS. The fragility of the full quantum state contrasts the robustness of local observables. (iii) In the small-τ regime, Trotter errors are perturbative in τ , which indicates that the Floquet-Magnus expansion converges, and truncating the expansion at low orders yields an accurate description of the system dynamics. While DQS is in general limited by experimental noise and decoherence, it is an interesting conceptual question whether in an "ideal" DQS Trotter errors could be controlled up to arbitrarily long times and in the thermodynamic limit. For collective spin models, we find that this is indeed the case. (iv) The threshold in Trotter errors, which marks the diver-gence of the FM expansion, is determined by the onset of quantum chaos in the kicked top. (v) Quantum chaos in the kicked top can be demonstrated experimentally by measuring the initial growth rate as well as the longtime steady state value of OTOCs. Accessing this rather novel diagnostic of quantum chaos is an intriguing possibility for the various experimental implementations of the kicked top. As we show in Appendix E, this is possible with the experimental resources that are required to generate the dynamics of the kicked top in the first place. (vi) Finally, a sharp threshold to quantum chaos persists across a wide range of exponents α in spin chains with power-law interactions.
It is worthwhile to emphasize again the distinction between single-particle and many-body quantum chaos, which is central to the results obtained in this work. On one hand, the quantum kicked top has emerged as a paradigmatic model for single-particle quantum chaos upon increasing the kicking strength. In the semiclassical limit, the problem becomes chaotic in the classical sense of exponentially diverging trajectories upon slightly perturbing initial conditions. From a quantum many-body perspective, however, the kicked top is an integrable and therefore non-ergodic system, because its dynamics is restricted to a N + 1-dimensional subspace of the 2 Ndimensional many-body Hilbert space. In this regard, the kicked top is also a rather peculiar system for DQS, which aims at simulating quantum many-body chaotic dynamics in regimes in which the growth of entanglement prohibits efficient classical simulations. Nevertheless, in this work we show that the phenomenology of Trotter errors in this particular DQS is analogous to the behavior found in a generic many-body system in Ref. [12], can be explained in terms of quantum chaos in the kicked top, and is smoothly connected to "true" quantum manybody systems described by Eq. (4) with α > 0.
We note that for small α > 0 it is not yet fully understoood whether the long-range deformation in Eq. (4) becomes ergodic. In the kicked top, i.e., for α = 0, threshold behavior persists even in the thermodynamic limit and up to arbitrarily long times. On the other hand, in generic many-body systems such as the Ising chain Eq. (4) with sufficiently large α, recent works argued that in the thermodynamic limit periodic driving will always lead to indefinite heating with the simulation accuracy Eq. (8) Q E → 1 [53,54,56], while heating can be strongly suppressed with Q E < 1 during a prethermal regime [86][87][88][89][90][91][92][93]. This regime persists up to a time scale which grows exponentially with the driving frequency, and has recently been observed experimentally [94]. It is an interesting question for future studies how the above scenarios are connected upon tuning α.
Our theoretical predictions concerning Trotter errors and the onset of quantum chaos can be tested experimentally in a variety of different systems ranging from single atomic spins [38,39] to quantum simulators of interacting spin-1 /2 systems [2-11, 47, 49, 50, 95-99]. A first experimental requirement is the ability to engineer non-linear spin Hamiltonians with time-dependent coefficients. Moreover, in any experimental realization, decoherence will eventually wash out the threshold shown in Fig. 1: Experimental imperfections such as technical noise or coupling of the system to its environment will lead to further heating and drive the system towards a featureless infinite-temperature state. In such a state, expectation values of observables become independent of the Trotter step size and the sharp threshold is destroyed. Therefore, the decoherence time t c should be longer than the time scale t th it takes to build up the threshold. For the parameters chosen in Fig. 1, this requirement is J z t th ≈ 25 < J z t c and can be met in current experiments. (To name an example, Ref. [41] achieved J z on the order of several kHz and thus nearly two orders of magnitude larger than the decoherence rate 1/t c of several ten Hz.) In this work, we exploit the formal equivalence of the dynamics of the kicked top which is determined by the Floquet operator given in Eq. (2), and the Trotterized dynamics in DQS of collective spin systems as described by Eq. (3). For a given target Hamiltonian H = H x + H z which is composed of two non-commuting pieces H x,z , Trotterization is not unique. In the following, we briefly discuss different Trotterization schemes, and how they affect our conclusions concerning Trotter errors and quantum chaos.
Trotterization of the unitary time evolution operator U (t) = e −iHt with H = H x + H z is carried out in two steps: First, the simulation time t is divided into n Trotter steps of duration τ = t/n, Second, within a Trotter step of size τ , the exponential e −iHτ of the Hamiltonian H is approximated by a product of exponentials of H x and H z alone. This approximation introduces an error of order τ m , where m can be increased by working with more elaborate factorizations [101,102]. Equation (3) corresponds to the lowest-order decomposition with m = 1, and a first improvement can be obtained by symmetrization, We emphasize that the given error bounds O(τ m ) pertain to the full time evolution operator. In contrast, as shown in this work for collective spin systems and in Ref. [12] for a generic many-body system, Trotter errors of fewbody observables remain bounded even at long times, i.e., after repeated application of the evolution operators in Eqs. (A2) and (A3) corresponding to elementary Trotter steps. However, it is an interesting question for future studies how the threshold behavior in Trotter errors of few-body observables and the quantum chaos properties discussed in Secs. II and III, respectively, are affected by higher-order Trotterization. Notably, already going from the first-order formula Eq. (A2) to the second-order formula Eq. (A3) leads to a significant qualitative difference in the PR Eq. (12) shown in Fig. 2(a). This is discussed in detail in Sec. III.

Appendix B: RMT predictions for the IPR, observables and OTOCs
For large Trotter step sizes τ , the Floquet operator Eq. (2) is faithful to the COE of RMT. This is demonstrated in Sec. III, where we consider the participation ratio PR and eigenphase-spacing statistics, cf. Fig. 2. However, also the results for the long-time averages of magnetization error and the simulation accuracy presented in Sec. II, as well as steady-state saturation of OTOCs discussed in Sec. IV, are indicative of the RMT properties of the Floquet operator. Here, we derive RMT predictions for these quantities.
RMT predictions for the IPR defined in Eq. (11) can be obtained by averaging the overlap | ψ 0 |φ m | 4 over the distribution of eigenvectors |φ m of random matrices drawn from the appropriate ensemble. In particular, the components of eigenvectors of random matrices sampled from the CUE are distributed uniformly on the unit sphere in C D , while for the COE the components of eigenvectors can be assumed real and thus lie on the the unit sphere in R D [14]. Averages over these eigenvector distributions can be performed by explicitly parameterizing unit-norm eigenvectors in terms of spherical coordinates and carrying out the resulting integrals as discussed for the CUE in Ref. [103]. These averages yield the same result for all states |φ m with m = 1, . . . , D and also does not depend on the reference state |ψ 0 . We note that for the COE, |ψ 0 should be chosen to have real components. This is indeed the case for the eigenstates |ψ n of the target Hamiltonian which we consider in Eq. (12), if they are expanded in the basis of eigenstates of J z . We find The RMT predictions for the PR as defined in Eq. (12) follow from the above results for the IPR as PR = 1/(D IPR), and are shown in Fig. 2(a). With regard to the simulation accuracy Q E defined in Eq. (8), the long-time average value of Q E = 1 which is reached in the chaotic regime as shown in Fig. 1(b,e), indicates that the temporal average of the energy E τ (t) = H(t) τ is equivalent to taking the expectation value of H in an infinite-temperature state, i.e., E T =∞ = tr(H)/D. As we show now, the same result is predicted by RMT. To this end, we replace in Eq. (8) the n t -th power of Floquet operator, U nt τ with n t = t/τ , by Q † U Q where Q is given in Eq. (15) and U is a random unitary matrix from the COE. As discussed in Sec. III, the COE is comprised of symmetric matrices, and the unitary transformation Q acts to symmetrize the Floquet operator U τ . Then, the average over U ∈ COE, which we denote by [· · · ] COE in the following, yields [104] whereH = QHQ † . Hence, we obtain the RMT prediction for the simulation accuracy in the chaotic regime, Similarly, Figs. 3(a,b) show that at late times the OTOC C(t) Eq. (17) approaches a value C COE which as we explain in the following can also be obtained by taking an average over the COE. First, we consider the COE average of the OTOC F (t) defined as where W (t) = U nt † τ W U nt τ and n t = t/τ . We assume tr(W ) = 0 and V † = V −1 , which is satisfied for W = S z and V = e −iϕSz shown in Fig. 3. As above, we replace U nt τ by Q † U Q, and by taking the average over U ∈ COE we find where N D = D 2 (D + 1) (D + 3) ,W = QW Q † , and V = QV Q † . To relate F (t) Eq. (B4) to the squared commutator in Eq. (17), we note that for W = S z we obtain tr(W 2 )/D = S (S + 1)/3, which yields This relation holds at all times t and, in particular, it allows us to express the COE average C COE in terms of F COE given in Eq. (B5). Moreover, the infinite-time average C shown in Fig. 3 can be obtained in terms of the corresponding average F of F (t), which reads where W nm = φ n |W |φ m etc. and the vectors |φ n with n = 1, . . . , D form an eigenbasis of the Floquet operator. Figure 3 shows that the RMT expression C COE of the OTOC agrees with the temporal average obtained from Eq. (B7) in the chaotic regime and beyond the Ehrenfest time t E .
Appendix C: The semiclassical limit of the kicked top As pointed out above, in the kicked top, the role of an effective Planck constant is played by the inverse spin size, eff = 1/S. Hence, if the kicked top is realized as the collective spin of a system of N spin-1 /2 such that S = N/2, the thermodynamic limit N → ∞ coincides with the semiclassical limit eff → 0. In this limit, the spin S obeys stroboscopic evolution equations that can be obtained from the Heisenberg equations of motion for the spin components S x,y,z by introducing rescaled variables, and subsequently taking the limit S → ∞. In this limit, the commutators [X, Y ] etc. vanish, and X, Y , and Z become effectively classical variables. The stroboscopic evolution equations of the spin operators in the Heisenberg picture are determined by To evaluate the right-hand side of this equation, we use relations such as [105] e iτ F (Sz) S + e −iτ F (Sz) = S + e iτ f (Sz) , where f (S z ) = F (S z + 1) − F (S z ), and which holds for any analytic function F and pairs of operators S z and S + = S x + iS y which obey [S z , S + ] = S + . We thus arrive at X(t + τ ) = Re e iτ hz X(t) + i Re (Y (t) + iZ(t))e iτ hx e iJxX(t)τ exp J z τ 2 Im (Y (t) + iZ(t))e iτ hx e iJxX(t)τ , Y (t + τ ) = Im e iτ hz X(t) + i Re (Y (t) + iZ(t))e iτ hx e iJxX(t)τ exp J z τ 2 Im (Y (t) + iZ(t))e iτ hx e iJxX(t)τ , Z(t + τ ) = Im (Y (t) + iZ(t))e iτ hx e iJxX(t)τ .

(C4)
To obtain the simulation accuracy Eq. (8) in the limit S → ∞, we iterate the semiclassical evolution equations (C4) and insert the resulting values of the spin components as per Eq. (C1) in the target Hamiltonian H = H x + H z , which we interpret as a function of the classical variables S x,y,z . This yields E τ (t). The energy E T =∞ at infinite temperature is obtained by averaging H over the classical phase space, i.e., a sphere of radius S.
Finally, E 0 is the energy of the initial spin configuration. The magnetization error Eq. (6) can be obtained similarly: In the semiclassical limit, S z (t) τ in Eq. (6) can be replaced by SZ(t), where Z(t) is evolved with Eq. (C4); To find S z (t) , we again interpret H as a classical Hamiltonian and integrate the corresponding Hamiltonian evolution equations up to the time t.

Appendix D: Perturbation Theory
In the regular regime of Trotter step sizes below τ * , corrections to the ideal long-time averages of observables due to Trotterization can be obtained by treating the difference between the effective Hamiltonian H τ obtained from the FM expansion and the target Hamiltonian H in time-dependent perturbation theory as described in Ref. [12]. Here, we generalize the discussion of Ref. [12], which assumed a fixed initial state |ψ 0 polarized along the z-axis, to an arbitrary choice of |ψ 0 .
The starting point is the FM expansion of the effective Hamiltonian H τ up to second order in τ (see Ref. [32] in the context of the kicked top) where the operator V = τ C 1 + τ 2 C 2 acts as a time-independent correction to the target Hamiltonian H = H x + H z . The effective Hamiltonian gives rise to a time evolution according to the Schrödinger equation which we transform into an interaction picture according to We obtain the time-evolution operator from time t 0 to t as where T denotes a time-ordered exponential. We determine W (t) = W (t, 0) (t 0 = 0) up to second order in the Trotter step τ .

Corrections to the energy
We are interested in an expansion of the quantity where we used Eq. (D1), in the Trotter step τ for an arbitrary initial state |ψ 0 . To this end, we replace e −iHτ t by W (t) which yields the time dependent correction up to second order in the Trotter step τ : The long-time averages of the above expectation values are determined in the diagonal ensemble [106], where the states |λ form a basis of eigenstates of the target Hamiltonian, H |λ = λ |λ . This simplifies the correction to This expression yields the dotted lines in Fig. 1(e).

Corrections to arbitrary observables
In the same way as before we derive an expansion up to second order in the Trotter step τ for the correction to an arbitrary observable A starting from an initial state |ψ 0 . After taking long-time averages in the diagonal ensemble we obtain the result Appendix E: Measurement of OTOCs using statistical correlations of randomized measurements As discussed in Sec. IV, the onset of chaos and thus the breakdown of DQS is diagnosed by the OTOC Eq. (17). In the following, we describe here how the OTOC can be measured in experiments.
The measurement of the OTOC defined in Eq. (17) is very challenging due to the alternation between operators at time zero and a finite time t. Protocols to measure C(t) (or F (t) defined in Eq. (B4)) that have been proposed in the literature thus involve backward time evolution [34,[77][78][79]107] or the coupling to auxiliary quantum systems [80], and the first measurements have been implemented only recently [41,48,81]. Here we propose to measure the OTOC by adapting a recently developed measurement protocol [82] which does not require the above ingredients. Instead, it is an example of measurement protocols which are based on statistical correlations between randomized measurements. This concept was first devised for Rényi entropies [108][109][110] and recently demonstrated in a trapped-ion quantum simulator [45]. The crucial ingredient in this class of measurement protocols are random unitary transformations which form (approximate) unitary two-designs. Thus, we first comment on the generation of such random unitaries in collective spin models and subsequently discuss how to use them to measure the OTOC.

Generation of random unitaries
In systems realizing collective spin models, the same tools that are used to engineer Hamiltonians of the form given in Eq. (7) can be applied to generate random unitaries from approximate unitary two-designs. Unitary two-designs are ensembles that approximate the circular unitary ensemble (CUE) up to 2nd order correlations [67]. We adapt ideas from Ref. [109,110] to create such random unitaries from a series of random quenches as a random time evolution operator u, where the Hamlitonians H α µ take the same form as in Eq. (7), For each quench, parameters h α µ and J α µ are chosen independently from uniform distributions h α µ ∈ [0.5/τ, 1.5/τ ], J α µ ∈ [5.5/τ, 6.5/τ ] such that the kicked top which would be realized by periodic repetition of any of the u α is deep in the chaotic regime [15].
To test whether the ensemble of thusly created random unitaries indeed forms a unitary two-design, we consider first their phase spacing statistics as characterized by the ratios of adjacent phase spacing [52]. As displayed in Fig. 5, we observe a very rapid convergence of the mean value r with the number of quenches η towards the CUE value r CUE = 0.5996 (1). Second, to test not only eigenvalues but also the eigenvectors of the created unitaries, we consider the second frame potential [67], which is defined as where E is an ensemble of unitary matrices. As shown in Ref. [67], for E being a unitary two-design, the frame potential takes a minimal value F CUE = 2!, and F represents a fine-grained measure of distance between an ensemble E and a unitary two-design. However, we note  Figure 5. Generation of approximate two-designs. Convergence of NU = 500 generated random unitaries to an approximate two-design as a function of the number of random quenches η for different spin sizes J. (a) Average adjacent phase spacing ratio. The black line corresponds toc the CUE value. (b) Difference of the frame potential FE g of the generated unitaries to the frame potential FE CUE = k! + d 4 /NU of an ensemble of same size |Eg| = |ECUE| = 500 which is drawn from the CUE. The black dashed line is the statistical error threshold ∼ 1/|Eg|. that for small ensembles |E| D 4 , the diagonal contribution D 4 /|E| in the sum Eq. (E3) dominates, irrespective of the properties of E. Thus, to test numerically using a finite sample size whether the generated random unitaries are sampled from an approximate unitary twodesign, we consider the quantity F ECUE = 2! + D 4 /|E CUE | is the average frame potential of a finite number |E CUE | of unitaries sampled from the CUE (see also Ref. [111]). In Fig. 5, F ECUE is displayed as function of number of quenches η, for a sample size |E CUE | = |E g | = 500. We observe rapid convergence of F (2) Eg − F (2) CUE towards the statistical error threshold ∼ 1/|E g |, signaling that the generated random unitaries are indeed sampled from an (approximate) unitary twodesign.

Protocol to measure the OTOC
The generated random unitaries, sampled from an approximate unitary two-design, can be used to measure 10 −2 C(t = τ n) J z τ 0.7 1.7 2.7 3.5 Figure 6. Measurement of the OTOC. Time evolution of the exact (solid lines) and "measured" (dots) squared commutator C(t) = 2 (S (S + 1)/3 − Re(F (t)) for a system with collective spin S = 64. Parameters are the same as in Fig. 2 of the main text. The "measurements" are simulated using NU = 100 random unitaries generated with η = 5 random quenches (same parameters as in Fig. 5). Projection noise is not included. The black line is the COE value.
the OTOC according to the following recipe [82]: (i) To an arbitrary pure state |ψ 0 , apply a random unitary u generated via random quenches (see Appendix E 1), to prepare the state |ψ u = u |ψ 0 . This randomized state serves as the input for two independent experiments.
(ii.a) Evolve the state in time with U (t) and measure the expectation value W (t) u = ψ u |U † (t)W U (t)|ψ u of the operator W . (ii.b) The state ψ u is perturbed by the operator V , evolved in time and W is measured, to obtain V † W (t)V u = ψ u |V † U † (t)W U (t)V |ψ u . (iii) Repeat steps (i) and (ii) for a sample of random unitaries u to estimate the statistical correlations of the two expectation values obtained in (ii.a) and (ii.b) across the unitary ensemble. As shown in Ref. [82], the OTOC can be expressed in terms of these statistical correlations as where · · · denotes the ensemble average over random unitaries u. This expression underpins the interpretation of F (t) as a measure of the butterfly effect in a quantum system: In analogy to rapidly diverging classical chaotic trajectories, the expectation values in the numerator can become decorrelated due to the perturbation V of the initial state. This is diagnosed by the decay of F (t). In Fig. 6, data points correspond to simulated experiments using 100 generated random unitaries. The data demonstrates that the protocol outlined above can serve as another experimental probe to delimit the region of applicability of DQS and the onset of quantum chaos.