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.


INTRODUCTION
In digital quantum simulation (DQS), unitary Hamiltonian evolution is decomposed into a sequence of quantum gates. A common approach to achieve this decomposition utilizes Suzuki-Trotter formulas 1,2 to approximately factorize the time evolution operator. [3][4][5][6][7][8][9][10][11][12] It is a fundamental conceptual question under which conditions this "Trotterization" is a controlled approximation. A recent numerical study 13 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. 14 Motivated by these findings we revisit the kicked top, a paradigmatic model of single-body quantum chaos. 15 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 timedependent Hamiltonian HðtÞ ¼ H x þ τH z X n2Z δðt À nτÞ; (1) which combines precession of the spin of the top S around a fixed axis, H x = h x S x , with nonlinear "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, 15 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 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. 15 Indeed, it is a defining feature of quantum chaotic systems that their spectral statistics are universal and obey predictions from random-matrix theory (RMT). 15 Among the model systems of quantum chaos, the kicked top stands out due to its extraordinary faithfulness to RMT. The discovery of this property initiated a surge of theoretical studies. Recent developments, [16][17][18][19] include proposals 20 to diagnose chaos in the kicked top by measuring out-of-time-ordered correlators, and the discovery of critical quasienergy states. 21 Experimentally, the kicked top was realized as the spin of single Cesium atoms 22 with S = 3, as the collective spin of an ensemble of three superconducting qubits 23 corresponding to S = 3/2, and very recently in nuclear magnetic resonance (NMR) as a spin S = 1 composed of two spin-1/2 nuclei. 24 Novel experimental possibilities include atomic species with high spin becoming available in labs such as Dysprosium 25 with S = 8 and Erbium 26 with S = 19/2, and also the spin S = 7/2 of 123 Sb nuclei; 26 implementations as collective spin in quantum simulators of all-to-all coupled spin-1/2 with "flip-flop" qubits, 27 or trapped ions 28,29 in which S ≳ 50 can be realized; and condensates of ultracold bosonic atoms 30 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 often "Trotterized", 1,2 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 ÀiHx t=n n ¼ U n τ : We note that other ways of decomposing the time evolution operator are being studied in the literature. 31,32 Further, Trotterization is not uniquely defined, and we discuss different schemes in Methods. 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 4,10,11,28,29,33 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, 34 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 nearest-neighbor interactions, α → ∞, and with the addition of a longitudinal field in Eq. (4), some of us addressed in ref. 13 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. 13 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. Ref. 13 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 14 typically have to resort to numerics (see refs 35,36 for exceptions). 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 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 size 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. 13 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 nonzero 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. 13 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.

RESULTS
Trotter errors in DQS of collective spin systems We consider two measures of Trotter errors: deviations of the expectation values of physical observables 13 and the fidelity of the quantum state obtained in DQS. Interestingly, these measures exhibit strikingly different behaviors.
We first compare the "magnetization", i.e., the expectation value of S z , under Trotterized and ideal dynamics starting from a spin coherent state θ; ϕ j i ¼ e iθ Sx sinðϕÞÀSy cosðϕÞ ð Þ S; S z ¼ S j i . Here, |S, S z 〉 is an eigenstate of the spin z-component with eigenvalue S. The magnetization error is defined as 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 P nt n¼1 ΔMðnτÞ where n t ¼ t=τ b c, which is shown in Fig. 1a. 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 fine-tuned 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 Discussion.
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. 16 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 τ R τ 0 dt HðtÞ k k. 37 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 timeaveraged 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Þ k k$ 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" 13 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 infinitetemperature state, Hilbert space dimension, and E 0 . A value of Q E (t) = 1 thus indicates heating to infinite temperature. Figure 1b shows the finite-time average, Q E ðtÞ ¼ 1 nt P 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. [38][39][40][41] Similarly, the robustness of observables also does not extend to the quantum state ψ τ ðtÞ j i¼ U nt τ ψ 0 j i obtained under Trotterized dynamics in DQS. In Fig. 1c, 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 〉 F ðtÞ ¼ ψ τ ðtÞjψðtÞ h i j j : 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. 1f. For large system sizes, the decay of the fidelity approaches F ðtÞ $ 1= ffiffiffiffi D p set by the Hilbert space dimension D ¼ 2S þ 1. 42 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  1d 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. 15 Therefore, in the thermodynamic limit, the spin components obey semiclassical stroboscopic evolution equations as detailed in Methods. These evolution equations can be iterated efficiently, and the long-time average for S → ∞ shown in Fig. 1d, e 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. 1d, 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 Fig. 1a, 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 below on these state-dependent variations, also of the threshold value τ * itself.
For small Trotter steps, the time-averaged magnetization error admits a controlled expansion in powers of τ, This expansion can be obtained analytically from the FM series for the effective Hamiltonian by employing time-dependent perturbation theory, 13 which we generalize to arbitrary initial states in the Supplementary Section 1. While we cannot justify the truncation of the FM expansion underlying Eq. (10) with full mathematical rigor, we take the excellent quantitative agreement between Eq. (10) and the numerical data shown in Fig. 1d as evidence that the scaling of the error predicted by truncating the FM expansion is indeed correct. Our numerical findings indicate that a fully rigorous justification of the analytical error estimate can be given, perhaps similarly to the problem of simulating local lattice Hamiltonians, 40 where an integral representation of Trotter errors for the full time evolution operator 41 proved the correctness of error estimates based on a truncation of the Baker-Campbell-Hausdorff formula. 39 Further, our findings suggest that the radius of convergence of the expansion Eq. (10) coincides with τ * and is finite.  28,29 Smaller (larger) values of S result in a more rounded (sharper) threshold. c, f Absence of a sharp threshold in the 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= ffiffi ffi S p 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 Q E . 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 Eq. (31). For a-c, f, the initial state is |θ, ϕ〉 = |0, 0〉 = |S, S〉, while in d, e the initial state is |θ, ϕ〉 = |0.25π, 0〉 Figure 1d shows data for spin sizes ranging form S = 8 to S → ∞. The smallest value shown, S = 8, is realized in experiments with Dysprosium atoms, 24 while S = 64 can be achieved in 1D 4,33 or 2D arrays of ions. 29 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. 1e, 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 the Supplementary Section 1.
The existence of controlled perturbative expansions for Trotter errors of few-body observables up to comparatively large values of the Trotter step size has profound practical implications for DQS. First, to obtain accurate results for few-body observables in DQS, one ideally would like to extrapolate from data at various finite τ to the limit τ → 0. Our results show that a controlled extrapolation is possible from a wide range of values τ < τ* with J z τ* of order one, and is not restricted by the commonly expected requirement J z τ ≪ 1. 3 Second, the ability to retain controlled Trotter errors with relatively large Trotter steps is highly beneficial for current experimental efforts in DQS, as the reduced number of required quantum gates mitigates the effects of imperfect gate operations. 13 The region of controlled Trotter errors at small Trotter step sizes gives way to a proliferation of errors at τ > τ * . In particular, the time-averaged simulation accuracy shown in Fig. 1e approaches the maximum value of Q E ¼ 1 with increasing spin size. In Methods we show that this value is compatible with the assumption that the Floquet operator Eq. (2) is a random unitary matrix. In other words, the Trotterized dynamics becomes completely unrelated to the ideal dynamics generated by the target Hamiltonian H and in this sense universal. As we detail in the following, this breakdown of Trotterization, which becomes sharp for S → ∞, marks the onset of quantum chaos.
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 wellknown to exhibit quantum chaos. 15 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. 15 Thence, the relevant RMT ensemble is the circular orthogonal ensemble (COE) of unitary and symmetric matrices.
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 PR In this definition, we rescaled the PR by a factor of 1=D so that it assumes values between two limiting cases PR ¼ D À1 fully localized; 1 fully delocalized: ( 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. 2a. 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 Fig. 1a, b, 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. 2a, 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 Methods, 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. (21). On the other hand, the symmetric Floquet operator U τ,s yields a value of the PR that matches the COE prediction obtained from Eq. (22). 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 Methods, 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.  2a, 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. 13 and further 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 τ Þ ¼ fe iθn g n and are shown in Fig. 2b 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 timeevolution 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. 43 In terms of the phase spacings δ n = θ n+1 − θ n , the average spacing ratio 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). 44 For the Floquet operator U τ in Eq.
(2), the adjacent phase spacing ratio is shown in Fig. 2c. 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 timeindependent Hamiltonians that have classical counterparts with only one degree of freedom exhibit nonuniversal level statistics. 14,15 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 two-dimensional phase space (θ, ϕ) ∈ [0, π) × [0, 2π). In the second regime in Fig. 2c, 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 few-body observables considered above is exactly analogous to the phenomenology of generic (short-range interacting) many-body systems, 13 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. 45 Therefore, the statistics of the eigenphases of the Floquet operator does not allow one to distinguish between localized and many-body chaotic regimes.
The PR 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 above 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 spincomponents of the top, we obtain in the semiclassical limit the Poincaré sections shown in Fig. 2d-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 L.M. Sieberer et al.
is reflected in the initial-state dependence of τ* in the thermodynamic limit as can be seen in Fig. 2g, where we show the simulation accuracy for trajectories which emanate from the initial states marked by colored dots in Fig. 2d-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 is quantified by out-of-time-ordered correlation functions (OTOCs) which, as we show in the Supplementary Section 2, can be measured in experimental realizations of the kicked top.
Signatures of dynamical localization and quantum chaos in out-oftime-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. [46][47][48] 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Þ y WUð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 ny τ WU n τ . In certain highly chaotic systems, 47-49 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. 50 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 high-energy, condensed matter, and AMO physics. 20,29,[47][48][49]51,52 In experimental realizations of the kicked top, the OTOC could be accessed through a recently introduced scheme which is based on correlations between randomized measurements 53 and has already been demonstrated in NMR. 54 We give details on the experimental requirements in the Supplementary Section 2.
As we show in the following, in the kicked top, both the shorttime 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) 20,29,53,55-59 as also discussed in the Supplementary Section 2. However, we note that for the value φ = 10 −4 shown in Fig. 3, the OTOC reduces to the form CðtÞ ¼ À φ 2 D trð½S z ðtÞ; S z 2 Þ þ Oðφ 4 Þ. The growth of the OTOC is depicted in Fig. 3a. As expected in a chaotic system, for large Trotter step sizes the OTOC grows exponentially, CðtÞ $ e λJzt . 51,60 This exponential growth extends up to the Ehrenfest time t E $ lnð h eff Þ j j , where the effective Planck constant is determined by the spin size, ħ eff = 1/S. 15 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 Methods.
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. 3a, 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. 52 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. 3b; (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. 3c. 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 Supplementary Section 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. 29 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 Fig. 3 Signatures of quantum chaos in OTOCs. a Short-time growth of the infinite-temperature OTOC Eq. (17) and b infinite-time average C. In the chaotic regime, the exponential growth of C(t) saturates to the COE value obtained from Eq. (25). The infinite-time average C follows from Eqs. (26) and (27). c Growth rate λ for the range of times in a during which the OTOC grows exponentially. The rate λ is extracted from the data by a fit to CðtÞ ¼ ae λJz t À b, where a, b, and λ are fit parameters. Error bars correspond to the fitting error. In the regular regime, the growth rate takes a value that is independent of the Trotter step size τ and agrees with the value at J z τ = 0, which we obtained from the ideal evolution with the time-averaged Hamiltonian. For Trotter step sizes in the chaotic region, the growth rate increases with τ 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. 4,10,11,28,29,33 Specifically, we study a kicked spin chain as defined in Eqs. (4) and (5). By tuning the exponent α, the longrange 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 j i ¼ N i¼1 " j i over up to 2 × 10 4 periods. For the quantities shown below we perform a temporal average over all periods.
As explained above, 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 above. On the other hand, for non-zero values of α, the dynamics explores the 2 N -dimensional many-body Hilbert space.
In Fig. 4a, 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. 4d shows cuts through the same data for six values of the exponent α form α = 0 to α = 2. 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 j i ¼ N i¼1 " j i). 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. 1e 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. 1e). Instead, for α ≳ 1/4, we find a clear crossover from a perturbative regime at small Trotter steps τ, implying controllable Trotter errors, to a many-body 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. 13 In Fig. 4b, 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 ! 1, see Methods. 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. (21). As discussed above, we expect that a symmetrized Floquet operator yields an IPR that is consistent with the COE result (22). 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 finite-size corrections, as is also visible in Fig. 4b. 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. 4e, which is based on the same data as panel 4b and shows cuts for several small values of α. Beyond the crossover window, the rate function shown in Fig. 4b 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 timeaveraged magnetization error Eq. (6) is shown in Fig. 4c. Again, we observe two different regimes depending on the Trotter step τ. Compared to the simulation accuracy and the IPR in Fig. 4b, the crossover region appears much sharper, which is also reflected in the cuts through the data shown in Fig. 4f. The overall behavior remains similar: delocalization of the time-evolved state in Hilbert space implies an error ΔM of order one as we find across a wide region in the α-τ plane. For small values of τ instead, ΔM is perturbatively small 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 many-body chaos threshold can exhibit slight modifications. The present study of the regime 0 ≤ α ≤ 3 in combination with the results of ref. 13 , which considered a short-range interacting system corresponding to the limit α → ∞, indicate that the transition to chaos is a generic phenomenon with broad relevance for DQS.

DISCUSSION
In the field of quantum chaos, it is well-known that the kicked top exhibits a transition from regular to chaotic dynamics upon increasing the kicking strength, and the signatures of this transition, e.g., in the spectral statistics of the Floquet operator, are well-understood. Our work connects these results to and highlights their fundamental importance for the seemingly remote field of quantum information science. Namely, we show that the transition to chaos becomes manifest in sharp threshold behavior of Trotter errors of few-body observables in DQS of collective spin systems. Contrary to the difference between the ideal and the Trotterized time-evolution operator, which is often used as a measure for Trotter errors in DQS and which grows at best linearly with simulation time, [38][39][40][41] Trotter errors of observables remain controlled and constant up to arbitrarily long times. A quantitatively accurate description of these Trotter errors is provided in the regular regime by the FM expansion-even up to Trotter step sizes for which one would expect the FM expansion to diverge 37 . In the chaotic regime, the kicked top is well-known to be faithful to RMT. 15 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 singleparticle 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 N -dimensional manybody Hilbert space. In this regard, the kicked top is also a rather peculiar system for DQS, which aims at simulating quantum manybody 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. 13 , can be explained in terms of quantum chaos in the kicked top, and is smoothly connected to "true" quantum many-body 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, 45,61,62 while heating can be strongly suppressed with Q E < 1 during a prethermal regime. [63][64][65][66][67] This regime persists up to a time scale which grows exponentially with the driving frequency, and has recently been observed experimentally. 68 It is an interesting question for future studies how the above scenarios are connected upon tuning α. Irrespective of this question, the fact that the heating time scale is infinite in the kicked top and can be made arbitrarily large in generic short-range interacting systems implies that DQS experiments are not limited by Trotter errors but rather by extrinsic error sources such as qubit decoherence. We emphasize that this conclusion is not restricted to all-to-all interacting systems. Rather, we expect it to hold generically for a large class of model systems which are relevant for DQS.
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 24,25 to quantum simulators of interacting spin-1/2 systems. [4][5][6][7][8][9][10][11][12]33 A first experimental requirement is the ability to engineer nonlinear 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. 29 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). It is an intriguing prospect to experimentally connect the threshold in Trotter errors to quantum chaos by measuring the growth rate and saturation value of the OTOC.
To summarize, in this work, we connect the well-studied quantum chaos of the kicked top with Trotter errors in DQS of L.M. Sieberer et al. 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 many-body systems and many-body quantum chaos presented in ref. 13

Higher-order Trotterization
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 noncommuting 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. 69 Eq. (3) corresponds to the lowest-order decomposition with m = 1, and a first improvement can be obtained by symmetrization e ÀiHτ ¼ e ÀiHx τ e ÀiHz τ þ OðτÞ; e ÀiHτ ¼ e ÀiHx τ=2 e ÀiHz τ e ÀiHx τ=2 þ Oðτ 2 Þ: 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. 13 for a generic many-body system, Trotter errors of few-body observables remain bounded even at long times, i.e., after repeated application of the evolution operators in Eqs. (19) and (20) 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 Results are affected by higher-order Trotterization. Notably, already going from the first-order formula Eq. (19) to the second-order formula Eq. (20) leads to a significant qualitative difference in the PR Eq. (12) shown in Fig.  2a. This is discussed in detail in Results.

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 Results, where we consider the 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 Results, as well as steady-state saturation of OTOCs, 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 . 15 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. 70 . These averages yield the same result for all states |ϕ m 〉 with m ¼ 1; ; D and also do 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. 2a.
With regard to the simulation accuracy Q E defined in Eq. (8), the longtime average value of Q E ¼ 1 which is reached in the chaotic regime as shown in Fig. 1b, 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¼1 ¼ 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=τ b c, by Q y UQ, where Q is given in Eq. (15) and U is a random unitary matrix from the COE. As discussed in Results, 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 71 whereH ¼ QHQ y . Hence, we obtain the RMT prediction for the simulation accuracy in the chaotic regime, Q COE = (E COE − E 0 )/(E T=∞ − E 0 )~1 for D ! 1.
Similarly, Fig. 3a, b shows 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 y τ WU nt τ and n t ¼ t=τ b c. We assume trðWÞ ¼ 0 and V y ¼ 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 y UQ, and by taking the average over U ∈ COE we find F COE ¼ 1 N D D þ 2 ð Þ trṼW ÃṼ yW T À Á À Â þtrṼW TṼ yWÃ À Á þ jtrṼ À Á j 2 trWW y À Á À D þ 4 ð Þ trW 2 À Á À jtrṼW T À Á j 2 À jtrṼW Ã À Á j 2 À D þ 1 ð Þ trṼ À Á trṼ yWTW Ã À Á þ trṼ y À Á trṼW ÃWT À Á À Á Ã ; (25) where N D ¼ D 2 D þ 1 ð ÞDþ3 ð Þ ;W ¼ QWQ y , andṼ ¼ QVQ y . To relate F(t) Eq. (24) 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. (25). Moreover, the infinitetime 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. (27) in the chaotic regime and beyond the Ehrenfest time t E .
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 L.M. Sieberer et al.
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 S i ðt þ τÞ ¼ U y τ S i ðtÞU τ for i 2 fx; y; zg: To evaluate the right-hand side of this equation, we use relations such as 72 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 that obey [S z , S + ] = S + . We thus arrive at Xðt þ τÞ ¼ Re e iτhz XðtÞ þ iRe ðYðtÞ þ iZðtÞÞe iτhx e iJx XðtÞτ À Á Â Ã e Jz τ 2 Im ðYðtÞþiZðtÞÞe iτhx e iJx XðtÞτ ð Þ ; Yðt þ τÞ ¼ Im e iτhz XðtÞ þ iRe ðYðtÞ þ iZðtÞÞe iτhx e iJx XðtÞτ À Á Â Ã e Jz τ 2 Im ðYðtÞþiZðtÞÞe iτhx e iJx XðtÞτ ð Þ ; Zðt þ τÞ ¼ Im ðYðtÞ þ iZðtÞÞe iτhx e iJx XðtÞτ À Á : (31) To obtain the simulation accuracy Eq. (8) in the limit S → ∞, we iterate the semiclassical evolution Eq. (31) and insert the resulting values of the spin components as per Eq. (28) 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. (31); 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. 73

DATA AVAILABILITY
The data sets generated and analyzed during the current study are available from the corresponding author on reasonable request.