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

TL;DR

This paper links Trotter errors in digital quantum simulation of collective spins to quantum chaos in the kicked top, showing how chaos onset causes error proliferation and proposing experimental platforms for studying this phenomenon.

Contribution

It provides a novel interpretation of Trotter errors through quantum chaos theory and demonstrates the feasibility of experimental studies across various quantum platforms.

Findings

Regular kicked top dynamics ensures Trotter error convergence.

Chaos onset correlates with rapid Trotter error growth.

Experimental platforms can probe the chaos-error relationship.

Abstract

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.

Figures

Equations

$$H(t)=H_x+\tau H_z\sum _{n\in \mathbb{Z}}\delta (t-n\tau ),$$

$$U_{\tau }={\mathrm{e}}^{-\mathrm{i}{H}_z\tau }{\mathrm{e}}^{-\mathrm{i}{H}_x\tau },$$

$$U(t)={\mathrm{e}}^{-\mathrm{i}Ht}\approx {\left({\mathrm{e}}^{-\mathrm{i}{H}_{z}t/n}{\mathrm{e}}^{-\mathrm{i}{H}_{x}t/n}\right)}^n=U_{\tau }^n.$$

$$H_x=h_x\sum _{i=1}^N{S}_i^x,H_z=J_{z,\alpha }\sum _{i<j=1}^N\frac{S_i^z{S}_j^z}{{\left|i-j\right|}^{\alpha }},$$

$$J_{z,\alpha}=J_{z}\bigg{/}\bigg{(}\frac{1}{N-1}\sum_{i<j=1}^{N}\frac{1}{\left\lvert i-j\right\rvert^{\alpha}}\bigg{)}.$$

$$\Delta M(t)=\frac{1}{S}\left(⟨S_z(t){⟩}_{\tau }-⟨S_z(t)⟩\right),$$

$$H_{\mu }=h_{\mu }{S}_{\mu }+J_{\mu }{S}_{\mu }^2/(2S+1),$$

$$Q_E(t)=\frac{E_{\tau }(t)-E_0}{E_{T=\infty }-E_0}$$

$$\mathcal{F}(t)=\left|⟨{\psi }_{\tau }(t)⟩\psi (t)\right|.$$

$$\overline{\Delta M}=\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{n^{\prime }=1}^n\Delta M(n^{\prime }\tau )=m_1\tau +m_2{\tau }^2+O({\tau }^3).$$

$$\mathrm{IPR}(|{\psi }_0⟩)=\sum _{m=1}^{\mathcal{D}}{\left|⟨{\psi }_0|{\varphi }_m⟩\right|}^4$$

$$\mathrm{PR}={\left(\sum _{n,m=1}^{\mathcal{D}}|⟨{\psi }_n|{\varphi }_m⟩{|}^4\right)}^{-1}.$$

$$\mathrm{PR}=\{\begin{array}{ll}{\mathcal{D}}^{-1}& \text{fully localized,}\\ 1& \text{fully delocalized.}\end{array}$$

$$U_{\tau ,\mathrm{s}}={\mathrm{e}}^{-\mathrm{i}{H}_x\tau /2}{\mathrm{e}}^{-\mathrm{i}{H}_z\tau }{\mathrm{e}}^{-\mathrm{i}{H}_x\tau /2},$$

$$U_{\tau }=Q^{†}{U}_{\tau ,\mathrm{s}}Q,Q={\mathrm{e}}^{-\mathrm{i}{H}_x\tau /2}.$$

$$r_n=\frac{\min ({\delta }_n,{\delta }_{n+1})}{\max ({\delta }_n,{\delta }_{n+1})},r=\frac{1}{\mathcal{D}}\sum _{n=1}^{\mathcal{D}}{r}_n$$

$$C(t)=\frac{1}{\mathcal{D}}\mathrm{tr}\left({\left|[W(t),V]\right|}^2\right),$$

$$U(t)={\mathrm{e}}^{-\mathrm{i}Ht}={\left({\mathrm{e}}^{-\mathrm{i}Ht/n}\right)}^n.$$

$${\mathrm{e}}^{-\mathrm{i}H\tau }$$

$${\mathrm{IPR}}_{\mathrm{CUE}}$$

$${\mathrm{IPR}}_{\mathrm{COE}}$$

$$\begin{array}{rl}{E}_{\mathrm{COE}}& ={\left[⟨{\psi }_0⟩Q^{†}{U}^{†}QH{Q}^{†}UQ{\psi }_0\right]}_{\mathrm{COE}}\\ & =\frac{1}{\mathcal{D}+1}\left(\mathrm{tr}(H)+⟨{\psi }_0⟩Q^{†}{\tilde{H}}^{T}Q{\psi }_0\right)\\ & \sim \frac{1}{\mathcal{D}}\mathrm{tr}(H)=E_{T=\infty }\text{for}\mathcal{D}\to \infty ,\end{array}$$

$$F(t)=\frac{1}{\mathcal{D}}\mathrm{tr}\left(W(t{)}^{†}{V}^{†}W(t)V\right),$$

$$F_{\mathrm{COE}}=\frac{1}{\mathcal{N}_{\mathcal{D}}}\left[\left(\mathcal{D}+2\right)\left(\mathop{\mathrm{tr}}(\tilde{V}\tilde{W}^{*}\tilde{V}^{\dagger}\tilde{W}^{T})\right.\right.\\ \left.+\mathop{\mathrm{tr}}(\tilde{V}\tilde{W}^{T}\tilde{V}^{\dagger}\tilde{W}^{*})+\bigl{|}\mathop{\mathrm{tr}}(\tilde{V})\bigr{|}^{2}\mathop{\mathrm{tr}}(\tilde{W}\tilde{W}^{\dagger})\right)\\ -\left(\mathcal{D}+4\right)\mathop{\mathrm{tr}}(\tilde{W}^{2})-\bigl{|}\mathop{\mathrm{tr}}(\tilde{V}\tilde{W}^{T})\bigr{|}^{2}-\bigl{|}\mathop{\mathrm{tr}}(\tilde{V}\tilde{W}^{*})\bigr{|}^{2}\\ \left.-\left(\mathcal{D}+1\right)\left(\mathop{\mathrm{tr}}(\tilde{V})\mathop{\mathrm{tr}}(\tilde{V}^{\dagger}\tilde{W}^{T}\tilde{W}^{*})+\mathop{\mathrm{tr}}(\tilde{V}^{\dagger})\mathop{\mathrm{tr}}(\tilde{V}\tilde{W}^{*}\tilde{W}^{T})\right)\right],$$

$$C(t)=2\left(\frac{S\left(S+1\right)}{3}-\mathrm{Re}(F(t))\right).$$

$$\begin{gathered} \overline{F}=\frac{1}{\mathcal{D}}\sum _{m,n=1}^{\mathcal{D}}\frac{1}{1+{\delta }_{m,n}}(W_{mm}{V}_{mn}^{†}{W}_{nn}{V}_{nm} \\ +W_{mn}{V}_{nn}^{†}{W}_{nm}{V}_{mm}), \end{gathered}$$

$$X=S_x/S,Y=S_y/S,Z=S_z/S,$$

$$S_i(t+\tau )=U_{\tau }^{†}{S}_i(t)U_{\tau }\text{for}i\in \{x,y,z\}.$$

$${\mathrm{e}}^{\mathrm{i}\tau F(S_z)}{S}_{+}{\mathrm{e}}^{-\mathrm{i}\tau F(S_z)}=S_{+}{\mathrm{e}}^{\mathrm{i}\tau f(S_z)},$$

$$\begin{array}{rl}X(t+\tau )& =\mathrm{Re}\left({\mathrm{e}}^{\mathrm{i}\tau h_z}\left[X(t)+\mathrm{i}\mathrm{Re}\left((Y(t)+\mathrm{i}Z(t)){\mathrm{e}}^{\mathrm{i}\tau h_x}{\mathrm{e}}^{\mathrm{i}{J}_{x}X(t)\tau }\right)\right]{\mathrm{e}}^{\frac{J_z\tau }{2}\mathrm{Im}\left((Y(t)+\mathrm{i}Z(t)){\mathrm{e}}^{\mathrm{i}\tau h_x}{\mathrm{e}}^{\mathrm{i}{J}_{x}X(t)\tau }\right)}\right),\\ Y(t+\tau )& =\mathrm{Im}\left({\mathrm{e}}^{\mathrm{i}\tau h_z}\left[X(t)+\mathrm{i}\mathrm{Re}\left((Y(t)+\mathrm{i}Z(t)){\mathrm{e}}^{\mathrm{i}\tau h_x}{\mathrm{e}}^{\mathrm{i}{J}_{x}X(t)\tau }\right)\right]{\mathrm{e}}^{\frac{J_z\tau }{2}\mathrm{Im}\left((Y(t)+\mathrm{i}Z(t)){\mathrm{e}}^{\mathrm{i}\tau h_x}{\mathrm{e}}^{\mathrm{i}{J}_{x}X(t)\tau }\right)}\right),\\ Z(t+\tau )& =\mathrm{Im}\left((Y(t)+\mathrm{i}Z(t)){\mathrm{e}}^{\mathrm{i}\tau h_x}{\mathrm{e}}^{\mathrm{i}{J}_{x}X(t)\tau }\right).\end{array}$$

Full text

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 Trotter1959 ; Suzuki1976 to approximately factorize the time evolution operator Lloyd1073 ; Lanyon57 ; Peng2005 ; Barends2016 ; Langford2016 ; OMalley2016 ; Barends2015 ; Barreiro2011 ; Martinez2016 ; Salathe2015 . It is a fundamental conceptual question under which conditions this “Trotterization” is a controlled approximation. A recent numerical study Heyl2019 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 DAlessio2016 . Motivated by these findings we revisit the kicked top, a paradigmatic model of single-body quantum chaos Haake2010 . 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

[TABLE]

which combines precession of the spin of the top $S$ around a fixed axis, $H_{x}=h_{x}S_{x}$, with non-linear “kicks” given by $H_{z}=J_{z}S_{z}^{2}/(2S+1)$, which are applied periodically at times $t=n\tau$ for all $n\in\mathbb{Z}$. Here, $S_{\mu}$ with $\mu=x,y,z$ are quantum angular momentum operators. The inverse spin size can be regarded as an effective Planck constant, $\hbar_{\mathrm{eff}}=1/S$ Haake2010 , and in the limit $S\to\infty$, 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 Haake2010 . This classical chaotic motion is reflected in the spectrum of the Floquet operator (we set the “true” Planck constant to unity, $\hbar=1$),

[TABLE]

which determines the evolution of a quantum kicked top during a period of duration $\tau$: In the chaotic regime of strong kicking, the spectral statistics of $U_{\tau}$ is described by Dyson’s ensemble of random orthogonal matrices Haake2010 . Indeed, it is a defining feature of quantum chaotic systems that their spectral statistics are universal and obey predictions from random-matrix theory (RMT) Haake2010 . 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 Bandyopadhyay2015 ; Bhosale2017 ; Bhosale2018 ; Kumari2018 include proposals Swingle2016 to diagnose chaos in the kicked top by measuring out-of-time-ordered correlators, and the discovery of critical quasienergy states Bastidas2014 . Experimentally, the kicked top was realized as the spin of single Caesium atoms Chaudhury2009 with $S=3$, as the collective spin of an ensemble of three superconducting qubits Neill2016 corresponding to $S=3/2$, and very recently in NMR as a spin $S=1$ composed of two spin-$\nicefrac{{1}}{{2}}$ nuclei Krithika2019 . Novel experimental possibilities include atomic species with high spin becoming available in labs such as Dysprosium Chalopin2018 with $S=8$ and Erbium Baier2018 with $S=19/2$, and also the spin $S=7/2$ of ${}^{123}\mathrm{Sb}$ nuclei Mourik2018 ; implementations as collective spin in quantum simulators of all-to-all coupled spin-$\nicefrac{{1}}{{2}}$ with “flip-flop” qubits Tosi2017 , or trapped ions Britton2012 ; Garttner2017 in which $S\gtrsim 50$ can be realized; and condensates of ultracold bosonic atoms Strobel2014 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_{\tau}$ in Eq. (2), is equivalent to DQS of a system with Hamiltonian $H=H_{x}+H_{z}$: In DQS, time evolution is often “Trotterized” Trotter1959 ; Suzuki1976 , 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,

[TABLE]

We note that other ways of decomposing the time evolution operator are being studied in the literature Childs2012 ; Berry2015 . 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, $\tau=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 Lanyon57 ; Barreiro2011 ; Martinez2016 ; Zhang2017 ; Britton2012 ; Garttner2017 is the Ising chain described by

[TABLE]

with spin operators $S_{i}^{\mu}$ where $\mu=x,y,z$ and $i=1,\dotsc,N$, and the power-law coupling strength $J_{z,\alpha}$. For such long-range interacting systems, extensivity of the Hamiltonian is ensured by the Kac prescription Kac1963 ,

[TABLE]

In the limit $\alpha\to 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, $\alpha\to\infty$, and with the addition of a longitudinal field in Eq. (4), some of us addressed in Ref. Heyl2019 the crucial question of “Trotter errors,” i.e., errors in DQS due to the approximate factorization $\mathrm{e}^{-\mathrm{i}H\tau}\approx\mathrm{e}^{-\mathrm{i}H_{z}\tau}\mathrm{e}^{-\mathrm{i}H_{x}\tau}$. By performing extensive numerical simulations, Ref. Heyl2019 found that the Trotter error in observables such as the magnetization exhibits threshold behavior as the Trotter step size $\tau$ is changed. Below the threshold, the Trotter error admits a perturbative expansion in $\tau$, while it is uncontrolled above the threshold. Reference Heyl2019 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 DAlessio2016 typically have to resort to numerics (see Refs. Kos2018 ; Bertini2018 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 $\alpha\to 0$ of the Hamiltonian Eq. (4). In this limit, the collective spin $S^{2}=S_{x}^{2}+S_{y}^{2}+S_{z}^{2}$ with components $S_{\mu}=\sum_{i=1}^{N}S_{i}^{\mu}$ where $\mu=x,y,z$ becomes a constant of motion. Accordingly, the $2^{N}$-dimensional Hilbert space of a system of $N$ spin-$\nicefrac{{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}\to 2\left(N+1\right)J_{z}/N$ to accommodate usual conventions). The dimension $\mathcal{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. Heyl2019 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 Heyl2019 . In particular, we explore the influence of deformations of the kicked top to $\alpha>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\leq\alpha\leq 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 Heyl2019 and the fidelity of the quantum stateobtained 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 $\ket{\theta,\phi}=\mathrm{e}^{\mathrm{i}\theta\left(S_{x}\sin(\phi)-S_{y}\cos(\phi)\right)}\ket{S,S_{z}=S}$. Here, $\ket{S,S_{z}}$ is an eigenstate of the spin $z$-component with eigenvalue $S_{z}$. The magnetization error is defined as

[TABLE]

where the expectation value $\langle\cdot\rangle_{\tau}$ is taken under Trotterized dynamics with step size $\tau$ while $\langle\cdot\rangle$ pertains to time evolution with the target Hamiltonian $H$. The accumulated effect of Trotter errors is quantified by the stroboscopic temporal average $\overline{\Delta M}(t)=\frac{1}{n_{t}}\sum_{n=1}^{n_{t}}\Delta M(n\tau)$ where $n_{t}=\left\lfloor t/\tau\right\rfloor$, which is shown in Fig. 1a. For this figure and throughout the paper, we assume that

[TABLE]

for both $\mu=x$ and $\mu=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_{\mathrm{th}}\approx 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 $\tau$ below a threshold value $J_{z}\tau_{*}\approx 2.5$, errors become significant for larger $\tau$. 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_{\mathrm{th}}$, where $t_{c}$ is the coherence time. This is discussed further in Discussion.

As detailed below, $\tau_{*}$ determines the radius of convergence of the Floquet-Magnus (FM) expansion for the effective Hamiltonian $H_{\tau}$ defined by $U_{\tau}=\mathrm{e}^{-\mathrm{i}H_{\tau}\tau}$ (see Ref. Bandyopadhyay2015 for the FM expansion in the context of the kicked top). The FM expansion is guaranteed to converge always if the “driving frequency” $\Omega=2\pi/\tau$ exceeds the time-averaged operator norm of the Hamiltonian, i.e., $\Omega\gtrsim\frac{1}{\tau}\int_{0}^{\tau}dt\left\lVert H(t)\right\rVert$ Casas2001 . 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 $\tau$ in which this criterion applies vanishes for large spin sizes $S$ since $\left\lVert H(t)\right\rVert\sim S$ so that $\tau\lesssim 1/S$. However, we find in DQS of collective spin systems that heating is suppressed up to a Trotter step size $\tau_{*}$ even for large values of $S$. We quantify this absence of heating for $\tau<\tau_{*}$ by the “simulation accuracy” Heyl2019

[TABLE]

where $E_{\tau}(t)=\langle H(t)\rangle_{\tau}$, and $E_{0}=\braket{\psi_{0}}{H}{\psi_{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=\infty}=\mathop{\mathrm{tr}}(H)/\mathcal{D}$ where $\mathcal{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 1b shows the finite-time average, $\overline{Q_{E}}(t)=\frac{1}{n_{t}}\sum_{n=1}^{n_{t}}Q_{E}(n\tau)$, which clearly exhibits threshold behavior characterized by the same timescales $t_{\mathrm{th}}$ and $\tau_{*}$ as the magnetization error.

These results demonstrate that few-body observables are quite robust against Trotter errors for $\tau<\tau_{*}$. This is in stark contrast to the accuracy of the full unitary time-evolution operator: The difference $U(t)-U_{\tau}^{n}$, which quantifies the error made in the Trotterization in Eq. (3), grows at least linearly in both the simulation time $t=n\tau$ and the system size $N$ Berry2007 ; Jordan2012 ; Haah2018 ; Childs2019 . Similarly, the robustness of observables also does not extend to the quantum state $\ket{\psi_{\tau}(t)}=U_{\tau}^{n_{t}}\ket{\psi_{0}}$ obtained under Trotterized dynamics in DQS. In Fig. 1c, we show stroboscopic temporal averages of the fidelity $\mathcal{F}(t)$, i.e., the absolute value of the overlap of $\ket{\psi_{\tau}(t)}$ with the ideal target state $\ket{\psi(t)}=U(t)\ket{\psi_{0}}$,

[TABLE]

As illustrated in the figure, the temporal average $\overline{\mathcal{F}}(t)$ approaches the ideal value of $\overline{\mathcal{F}}(t)=1$ for $\tau\to 0$, but drops sharply already for $J_{z}\tau\ll 1$, and in particular for Trotter step sizes which are much smaller than the threshold value $\tau_{*}$ identified above for the magnetization error and the simulation accuracy. At this threshold value, i.e., for $\tau\approx\tau_{*}$, there is another noticeable drop in $\overline{\mathcal{F}}(t)$. However, also in the region below $\tau_{*}$ the fidelity vanishes with increasing system size as shown in Fig. 1f. For large system sizes, the decay of the fidelity approaches $\mathcal{F}(t)\sim 1/\sqrt{\mathcal{D}}$ set by the Hilbert space dimension $\mathcal{D}=2S+1$ Emerson2002 . 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<\infty$ is obtained by exact diagonalization of the Floquet operator $U_{\tau}$. In the limit $S\to\infty$, the effective Planck constant vanishes, $\hbar_{\mathrm{eff}}=1/S\to 0$ Haake2010 . 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\to\infty$ shown in Fig. 1 corresponds to $n=10^{6}$ iterations. Fluctuations in the data for $\overline{\Delta 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\to\infty$. We note that the different shapes of the curves for the magnetization error in Figs. 1a and d are due to the different choice of the initial state, which is $\ket{\theta,\phi}=\ket{0,0}=\ket{S,S}$ and $\ket{\theta,\phi}=\ket{0.25\pi,0}$ in a and d, respectively. We comment below on these state-dependent variations, also of the threshold value $\tau_{*}$ itself.

For small Trotter steps, the time-averaged magnetization error admits a controlled expansion in powers of $\tau$,

[TABLE]

This expansion can be obtained analytically from the FM series for the effective Hamiltonian by employing time-dependent perturbation theory Heyl2019 , 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 Haah2018 , where an integral representation of Trotter errors for the full time evolution operator Childs2019 proved the correctness of error estimates based on a truncation of the Baker-Campbell-Hausdorff formula Jordan2012 . Further, our findings suggest that the radius of convergence of the expansion Eq. (10) coincides with $\tau_{*}$ and is finite.

Figure 1d shows data for spin sizes ranging form $S=8$ to $S\to\infty$. The smallest value shown, $S=8$, is realized in experiments with Dysprosium atoms  Chalopin2018 , while $S=64$ can be achieved in 1D Lanyon57 ; Zhang2017 or 2D arrays of ions Garttner2017 . A threshold in $\overline{\Delta M}$ is visible over the entire range of values of $S$. The value of $\tau_{*}$ depends only very weakly on system size, and crucially it remains finite even for $S\to\infty$. These observations do not depend qualitatively on the choice of the Hamiltonian parameters or the angles $\theta$ and $\phi$ specifying the initial coherent spin state.

In Fig. 1e, we also show the long-time average $\overline{Q_{E}}$ of the simulation accuracy $Q_{E}(t)$ defined in Eq. (8). Clearly, heating is suppressed to perturbatively small values, $\overline{Q_{E}}=q_{1}\tau+q_{2}\tau^{2}+O(\tau^{3})$, for Trotter steps up to the same threshold value $\tau_{*}$ as for the magnetization error, indicating a finite radius of convergence of the FM expansion even for $S\to\infty$. 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 $\tau$ to the limit $\tau\to 0$. Our results show that a controlled extrapolation is possible from a wide range of values $\tau<\tau^{*}$ with $J_{z}\tau^{*}$ of order one, and is not restricted by the commonly expected requirement $J_{z}\tau\ll 1$ Lloyd1073 . 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 Heyl2019 .

The region of controlled Trotter errors at small Trotter step sizes gives way to a proliferation of errors at $\tau>\tau^{*}$. In particular, the time-averaged simulation accuracy shown in Fig. 1e approaches the maximum value of $\overline{Q_{E}}=1$ with increasing spin size. In Methods we show that this value is compatible with the assumption that the Floquet operator (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\to\infty$, 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 $\tau$, the kicked top is well-known to exhibit quantum chaos Haake2010 . 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 Haake2010 . 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 $\ket{\psi_{0}}$,

[TABLE]

where $\{|\phi_{m}\rangle\}_{m}$ denotes the eigenbasis of the Floquet operator $U_{\tau}$. We average the IPR over the basis of the target Hamiltonian $\{|\psi_{n}\rangle\}_{n}$ and take its inverse to obtain the participation ratio (PR),

[TABLE]

In this definition, we rescaled the PR by a factor of $1/\mathcal{D}$ so that it assumes values between two limiting cases:

[TABLE]

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 $\tau$. However, we would like to emphasize a crucial difference between the PR and Trotter errors of observables: The threshold value $\tau_{*}$ for Trotter errors depends on the initial state of the time evolution as illustrated Figs. 1a,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. 2a, two sets of curves are displayed: Solid lines pertain to the PR calculated with eigenvectors $\ket{\phi_{m}}$ of the Floquet operator $U_{\tau}$ given in Eq. (2); dashed lines correspond to a symmetric Floquet operator,

[TABLE]

which realizes a higher-order Trotter decomposition as explained in Methods, and is related to the Floquet operator $U_{\tau}$ in Eq. (2) by a unitary transformation,

[TABLE]

In the chaotic regime of large Trotter steps, the PR obtained from the Floquet operator $U_{\tau}$ Eq. (2) assumes the value predicted for the CUE of RMT, see Eq. (21). On the other hand, the symmetric Floquet operator $U_{\tau,\mathrm{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_{\tau,\mathrm{s}}$ is symmetric, $U_{\tau,\mathrm{s}}=U_{\tau,\mathrm{s}}^{T}$, it can be regarded as an element of the COE. The PR Eq. (12) probes the statistics of eigenvectors of the Floquet operator $U_{\tau,\mathrm{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_{\tau,\mathrm{s}}$ is related to the non-symmetric Floquet operator $U_{\tau}$ by the unitary transformation given in Eq. (15). Since a unitary transformation does not affect the spectrum of an operator, both $U_{\tau}$ and $U_{\tau,\mathrm{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_{\tau}$ 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 $\sim 1/\mathcal{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 $\ket{\psi_{m}}$ of $H$ are spread out over $\sim\mathcal{D}$ Floquet states $\ket{\phi_{m}}$. This is similar to the many-body case considered in Ref. Heyl2019 and further below, where for small $\tau$ 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 $\sigma(U_{\tau})$ of the Floquet operator. The eigenphases $\{\theta_{n}\}_{n}$ of the unitary operator $U_{\tau}$ are defined through $\sigma(U_{\tau})=\{\mathrm{e}^{\mathrm{i}\theta_{n}}\}_{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 $\tau$, the eigenphases of the Floquet operator $U_{\tau}$ coincide with the phases corresponding to the ideal time-evolution operator over one Trotter step, $\mathrm{e}^{-\mathrm{i}H\tau}$, which are shown as gray lines in the figure. (ii) In the second regime, the phases start to wind around the interval $[-\pi,\pi]$. 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$ Oganesyan2007 . In terms of the phase spacings $\delta_{n}=\theta_{n+1}-\theta_{n}$, the average spacing ratio is defined through

[TABLE]

where $\mathcal{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_{\mathrm{POI}}=2\ln(2)-1\approx 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_{\mathrm{COE}}=0.5996(1)$ Atas2013 . For the Floquet operator $U_{\tau}$ in Eq. (2), the adjacent phase spacing ratio is shown in Fig. 2c. Again, three regimes can be distinguished: In a region of size $\sim 1/S$ in which the spectral width of the target Hamiltonian $H$ is smaller than $\Omega=2\pi/\tau$, $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 Haake2010 ; DAlessio2016 . 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_{x}^{2}+S_{y}^{2}+S_{z}^{2}$, and can be parameterized in terms of the polar and azimuthal angles $\theta$ and $\phi$, respectively, which span the two-dimensional phase space $\left(\theta,\phi\right)\in[0,\pi)\times[0,2\pi)$. In the second regime in Fig. 2c, the eigenphases $\theta_{n}$ of the Floquet operator wind around the unit circle. As shown in the previous section, for values of $\tau$ up to $J_{z}\tau\lesssim 2$, the dynamics is well-captured by the time-independent effective Hamiltonian $H_{\tau}$ obtained in the FM expansion. For the average phase spacing ratio we find $r\approx 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 $\tau$, 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 Heyl2019 , 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 DAlessio2014b . 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_{\tau}$ 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 $\tau_{*}$ 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. 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 $\tau$. This coexistence is reflected in the initial-state dependence of $\tau_{*}$ 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 $\theta\in\{0,\pi\}$ are particularly stable whereas states around $\theta\approx\frac{\pi}{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-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 Larkin1969 ; Maldacena2016a ; Kitaev . The “infinite-temperature” OTOC for two initially commuting operators $V$ and $W$ is defined as

[TABLE]

where $\mathcal{D}$ is the Hilbert space dimension, and the Heisenberg operator $W(t)=U(t)^{\dagger}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 $\mathcal{D}=2S+1$ and we evaluate $C(t)$ stroboscopically at multiples of the Trotter step such that $t=n\tau$ and $W(n\tau)=U_{\tau}^{n\dagger}WU_{\tau}^{n}$. In certain highly chaotic systems Maldacena2016a ; Kitaev ; Hosur2016 , 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 Roberts2017 . 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 Maldacena2016a ; Kitaev ; Hosur2016 ; Pappalardi2018 ; Rozenbaum2017a ; Swingle2016 ; Garttner2017 . 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 Vermersch2018a and has already been demonstrated in NMR Nie2019 . We give details on the experimental requirements in the Supplementary Section 2.

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=\mathrm{e}^{-\mathrm{i}\varphi S_{z}}$ and Hermitian $W=S_{z}$—a particular choice of operators that facilitates the measurement of $C(t)$ Swingle2016 ; Zhu2016 ; Shen2017 ; Yoshida2018 ; Yao2016 ; Li2017 ; Garttner2017 ; Vermersch2018a as also discussed in the Supplementary Section 2. However, we note that for the value $\varphi=10^{-4}$ shown in Fig. 3, the OTOC reduces to the form $C(t)=-\frac{\varphi^{2}}{\mathcal{D}}\mathop{\mathrm{tr}}([S_{z}(t),S_{z}]^{2})+O(\varphi^{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)\sim\mathrm{e}^{\lambda J_{z}t}$ Pappalardi2018 ; Seshadri2018 . This exponential growth extends up to the Ehrenfest time $t_{\mathrm{E}}\sim\left\lvert\ln(\hbar_{\mathrm{eff}})\right\rvert$, where the effective Planck constant is determined by the spin size, $\hbar_{\mathrm{eff}}=1/S$ Haake2010 . For times $t$ beyond the Ehrenfest time $t_{\mathrm{E}}$, the OTOC saturates at a value $C(t)=C_{\mathrm{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 non-zero Trotter step size $J_{z}\tau=0.73$ shown in Fig. 3a, the OTOC is essentially indistinguishable from evolution under $U(t)=\mathrm{e}^{-\mathrm{i}Ht}$ with the target Hamiltonian $H=H_{x}+H_{z}$, which is labeled as $J_{z}\tau=0$ in the figure. From a short-time expansion of the time evolution operator as $\mathrm{e}^{-\mathrm{i}Ht}=1-\mathrm{i}Ht+O(t^{2})$ it follows immediately that at very short times the OTOC grows algebraically, $C(t)\sim 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. Rozenbaum2017a 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_{\mathrm{COE}}$ as shown in Fig. 3b; (ii) Moreover, the growth rate $\lambda$, 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 $\tau$ 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 $\lambda$ very accurately due to the large Ehrenfest time $t_{\mathrm{E}}\sim\ln(S)$. However, we found (see Supplementary Section 2) that the two regimes in the dependence of $\lambda$ on the Trotter step size—$\lambda\approx\mathrm{const.}$ in the regular regime and strong increase of $\lambda$ with $\tau$ in the chaotic regime—can be clearly distinguished also in systems with $S\gtrsim 50$ as in recent measurements of OTOCs with trapped ions Garttner2017 .

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 Lanyon57 ; Barreiro2011 ; Martinez2016 ; Zhang2017 ; Britton2012 ; Garttner2017 . Specifically, we study a kicked spin chain as defined in Eqs. (4) and (5). By tuning the exponent $\alpha$, 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 $\alpha\to 0$ and $\alpha\to\infty$, 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 $\ket{\psi_{0}}=\bigotimes_{i=1}^{N}\ket{\uparrow}$ over up to $2\cdot 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-$\nicefrac{{1}}{{2}}$ can be decomposed into subspaces with fixed total spin $S$. The initial state $\ket{\psi_{0}}$ given above belongs to the subspace of maximal spin $S=N/2$, and for $\alpha\to 0$ the dynamics is restricted to this subspace. That is, for $\alpha\to 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 $\alpha$, 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 $\tau$ and interaction exponent $\alpha$ for a fixed system size $N=14$. To clearly illustrate how the kicked top is approached in the limit $\alpha\to 0$, Fig. 4d shows cuts through the same data for four values of the exponent $\alpha$ form $\alpha=0$ to $\alpha=0.38$. For small $\tau$, we find that the properties observed already with the kicked top extend to $\alpha>0$. Specifically, $\overline{Q_{E}}$ is small with $\overline{Q_{E}}=q_{2}\tau^{2}+O(\tau^{3})$ (the first-order term vanishes for the initial state $\ket{\psi_{0}}=\bigotimes_{i=1}^{N}\ket{\uparrow}$). Upon increasing the Trotter step for $\alpha\lesssim 1/4$, however, $\overline{Q_{E}}$ does not reach the value $\overline{Q_{E}}\to 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 $\alpha=0$, where significant finite-size effects in $\overline{Q_{E}}$ appear (note also the different choice of initial state and Hamiltonian parameters in Fig. 1e). Instead, for $\alpha\gtrsim 1/4$, we find a clear crossover from a perturbative regime at small Trotter steps $\tau$, implying controllable Trotter errors, to a many-body quantum chaotic region with $\overline{Q_{E}}\to 1$. Except for the range $0\leq\alpha\lesssim 1/4$, the influence of the interaction exponent $\alpha$ on the simulation accuracy is only minor in that also the crossover region only shifts slightly upon varying $\alpha$.

A similar picture emerges from the IPR defined in Eq. (11), which measures the localization properties of the state $\ket{\psi_{0}}$ in the eigenbasis $\{\ket{\phi_{m}}\}_{m}$ of the Floquet operator. For $\alpha>0$, the accessible Hilbert space grows exponentially with the number of spins, $\mathcal{D}=2^{N}$. To account for this exponential dependence, we introduce the rate function $\lambda_{\mathrm{IPR}}=\ln(\mathrm{IPR})/N$ Heyl2019 . In Fig. 4b, we show the ratio $\lambda_{\mathrm{IPR}}/\lambda_{\mathcal{D}}$, where $\lambda_{\mathcal{D}}=\ln(2)\left(1-3/N\right)$. This value incorporates leading-order finite-size corrections as follows: On the one hand the system obeys both an inversion as well as Ising $\mathbb{Z}_{2}$ symmetry, such that $\mathcal{D}=2^{N-2}$; On the other hand, RMT allows us to estimate the IPR in an ergodic phase as $\mathrm{IPR}\sim 2/\mathcal{D}$ for $\mathcal{D}\to\infty$, see Methods. Combining these two results we obtain the value of $\lambda_{\mathcal{D}}$ given above. We note that here we restrict ourselves to comparing the IPR obtained for the initial state $\ket{\psi_{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 $\lambda_{\mathrm{IPR}}$ essentially reproduces the observations from the simulation accuracy in the range $0<\alpha\leq 3$ with only slight modifications, which we again attribute to finite-size corrections. The case $\alpha$ close to $\alpha=0$, however, deserves a more detailed discussion. While a system at $\alpha=0$ with total spin $S=N/2$ can access only $\mathcal{D}=2S+1=N+1$ quantum states in the Hilbert space, with any infinitesimal deviation from $\alpha=0$ this number grows to $\mathcal{D}=2^{N}$. For a system with finite $N$ as we consider here, however, we expect that the effective behavior of $\mathcal{D}$ has to interpolate from the polynomial dependence at $\alpha=0$ to the exponential dependence in $N$ upon increasing $\alpha$. This opens up a crossover window at small $\alpha$, whose size diminishes for larger $N$ and where the rate function $\lambda_{\mathrm{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 $\lambda_{\mathrm{IPR}}$ one considers the IPR Eq. (11) itself. This is illustrated in Fig. 4e, which is based on the same data as panel b and shows cuts for several small values of $\alpha$. 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 time-averaged magnetization error Eq. (6) is shown in Fig. 4c. Again, we observe two different regimes depending on the Trotter step $\tau$. 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 $\overline{\Delta M}$ of order one as we find across a wide region in the $\alpha$-$\tau$ plane. For small values of $\tau$ instead, $\overline{\Delta 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\leq\alpha\leq 3$ in combination with the results of Ref. Heyl2019 , which considered a short-range interacting system corresponding to the limit $\alpha\to\infty$, 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 Berry2007 ; Jordan2012 ; Haah2018 ; Childs2019 , 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 Casas2001 . In the chaotic regime, the kicked top is well-known to be faithful to RMT Haake2010 .

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^{N}$-dimensional 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. Heyl2019 , 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 $\alpha>0$.

We note that for small $\alpha>0$ it is not yet fully understoood whether the long-range deformation in Eq. (4) becomes ergodic. In the kicked top, i.e., for $\alpha=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 $\alpha$, recent works argued that in the thermodynamic limit periodic driving will always lead to indefinite heating with the simulation accuracy Eq. (8) $Q_{E}\to 1$ DAlessio2014b ; Lazarides2014 ; Luitz2017 , while heating can be strongly suppressed with $Q_{E}<1$ during a prethermal regime Abanin2015b ; Mori2016 ; Kuwahara2016 ; Machado2017 ; Howell2018 . This regime persists up to a time scale which grows exponentially with the driving frequency, and has recently been observed experimentally Singh2018 . It is an interesting question for future studies how the above scenarios are connected upon tuning $\alpha$. 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 Chalopin2018 ; Baier2018 to quantum simulators of interacting spin-$\nicefrac{{1}}{{2}}$ systems Lanyon57 ; Peng2005 ; Barends2016 ; Langford2016 ; OMalley2016 ; Zhang2017 ; Barends2015 ; Barreiro2011 ; Martinez2016 ; Salathe2015 . 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_{\mathrm{th}}$ it takes to build up the threshold. For the parameters chosen in Fig. 1, this requirement is $J_{z}t_{\mathrm{th}}\approx 25<J_{z}t_{c}$ and can be met in current experiments. (To name an example, Ref. Garttner2017 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 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. Heyl2019 .

**Methods 

** 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 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)=\mathrm{e}^{-\mathrm{i}Ht}$ 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 $\tau=t/n$,

[TABLE]

Second, within a Trotter step of size $\tau$, the exponential $\mathrm{e}^{-\mathrm{i}H\tau}$ 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 $\tau^{m}$, where $m$ can be increased by working with more elaborate factorizations Wiebe2010 . Equation (3) corresponds to the lowest-order decomposition with $m=1$, and a first improvement can be obtained by symmetrization,

[TABLE]

We emphasize that the given error bounds $O(\tau^{m})$ pertain to the full time evolution operator. In contrast, as shown in this work for collective spin systems and in Ref. Heyl2019 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 $\tau$, the Floquet operator Eq. (2) is faithful to the COE of RMT. This is demonstrated in Results, 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 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 $\left\lvert\braket{\psi_{0}}{\phi_{m}}\right\rvert^{4}$ over the distribution of eigenvectors $\ket{\phi_{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 $\mathbb{C}^{\mathcal{D}}$, while for the COE the components of eigenvectors can be assumed real and thus lie on the the unit sphere in $\mathbb{R}^{\mathcal{D}}$ Haake2010 . 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. PhysRev.132.948 . These averages yield the same result for all states $\ket{\phi_{m}}$ with $m=1,\dotsc,\mathcal{D}$ and also do not depend on the reference state $\ket{\psi_{0}}$. We note that for the COE, $\ket{\psi_{0}}$ should be chosen to have real components. This is indeed the case for the eigenstates $\ket{\psi_{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

[TABLE]

The RMT predictions for the PR as defined in Eq. (12) follow from the above results for the IPR as $\mathrm{PR}=1/(\mathcal{D}\,\mathrm{IPR})$, and are shown in Fig. 2a.

With regard to the simulation accuracy $Q_{E}$ defined in Eq. (8), the long-time average value of $\overline{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_{\tau}(t)=\left\langle H(t)\right\rangle_{\tau}$ is equivalent to taking the expectation value of $H$ in an infinite-temperature state, i.e., $E_{T=\infty}=\mathop{\mathrm{tr}}(H)/\mathcal{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_{\tau}^{n_{t}}$ with $n_{t}=\lfloor t/\tau\rfloor$, by $Q^{\dagger}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_{\tau}$. Then, the average over $U\in\mathrm{COE}$, which we denote by $[\dotsm]_{\mathrm{COE}}$ in the following, yields Brouwer1996

[TABLE]

where $\tilde{H}=QHQ^{\dagger}$. Hence, we obtain the RMT prediction for the simulation accuracy in the chaotic regime, $Q_{\mathrm{COE}}=(E_{\mathrm{COE}}-E_{0})/(E_{T=\infty}-E_{0})\sim 1$ for $\mathcal{D}\to\infty$.

Similarly, Figs. 3a,b show that at late times the OTOC $C(t)$ Eq. (17) approaches a value $C_{\mathrm{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

[TABLE]

where $W(t)=U_{\tau}^{n_{t}\dagger}WU_{\tau}^{n_{t}}$ and $n_{t}=\lfloor t/\tau\rfloor$. We assume $\mathop{\mathrm{tr}}(W)=0$ and $V^{\dagger}=V^{-1}$, which is satisfied for $W=S_{z}$ and $V=\mathrm{e}^{-\mathrm{i}\varphi S_{z}}$ shown in Fig. 3. As above, we replace $U_{\tau}^{n_{t}}$ by $Q^{\dagger}UQ$, and by taking the average over $U\in\mathrm{COE}$ we find

[TABLE]

where $\mathcal{N}_{\mathcal{D}}=\mathcal{D}^{2}\left(\mathcal{D}+1\right)\left(\mathcal{D}+3\right),$ $\tilde{W}=QWQ^{\dagger}$, and $\tilde{V}=QVQ^{\dagger}$. To relate $F(t)$ Eq. (24) to the squared commutator in Eq. (17), we note that for $W=S_{z}$ we obtain $\mathop{\mathrm{tr}}(W^{2})/\mathcal{D}=S\left(S+1\right)\!/3$, which yields

[TABLE]

This relation holds at all times $t$ and, in particular, it allows us to express the COE average $C_{\mathrm{COE}}$ in terms of $F_{\mathrm{COE}}$ given in Eq. (25). Moreover, the infinite-time average $\overline{C}$ shown in Fig. 3 can be obtained in terms of the corresponding average $\overline{F}$ of $F(t)$, which reads

[TABLE]

where $W_{nm}=\braket{\phi_{n}}{W}{\phi_{m}}$ etc. and the vectors $\ket{\phi_{n}}$ with $n=1,\dotsc,\mathcal{D}$ form an eigenbasis of the Floquet operator. Figure 3 shows that the RMT expression $C_{\mathrm{COE}}$ of the OTOC agrees with the temporal average obtained from Eq. (27) in the chaotic regime and beyond the Ehrenfest time $t_{\mathrm{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, $\hbar_{\mathrm{eff}}=1/S$. Hence, if the kicked top is realized as the collective spin of a system of $N$ spin-$\nicefrac{{1}}{{2}}$ such that $S=N/2$, the thermodynamic limit $N\to\infty$ coincides with the semiclassical limit $\hbar_{\mathrm{eff}}\to 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,

[TABLE]

and subsequently taking the limit $S\to\infty$. 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

[TABLE]

To evaluate the right-hand side of this equation, we use relations such as Kitagawa1993

[TABLE]

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}+\mathrm{i}S_{y}$ that obey $[S_{z},S_{+}]=S_{+}$. We thus arrive at

[TABLE]

To obtain the simulation accuracy Eq. (8) in the limit $S\to\infty$, we iterate the semiclassical evolution equations (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_{\tau}(t)$. The energy $E_{T=\infty}$ 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, $\left\langle S_{z}(t)\right\rangle_{\tau}$ in Eq. (6) can be replaced by $SZ(t)$, where $Z(t)$ is evolved with Eq. (31); To find $\left\langle S_{z}(t)\right\rangle$, we again interpret $H$ as a classical Hamiltonian and integrate the corresponding Hamiltonian evolution equations up to the time $t$.

**Acknowledgments 

** We thank B. Vermersch, M. Baranov, and M. Bukov for helpful discussions. Work in Innsbruck is supported by the European Research Council (ERC) Synergy Grant UQUAM and the SFB FoQuS (FWF Project No. F4016-N23). MH acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) via the Gottfried Wilhelm Leibniz Prize program. PH acknowledges support by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT), the ERC Advanced Grant EntangleGen (Project-ID 694561), and the ERC Starting Grant StrEnQTh. Numerical simulations were realized with QuTiP Johansson20131234 .