Trotter errors from dynamical structural instabilities of Floquet maps in quantum simulation

We study the behavior of errors in the quantum simulation of spin systems with long-range multi-body interactions resulting from the Trotter-Suzuki decomposition of the time-evolution operator. We identify a regime where the Floquet operator underlying the Trotter decomposition undergoes sharp changes even for small variations in the simulation step size. This results in a time evolution operator that is very different from the dynamics generated by the targeted Hamiltonian, which leads to a proliferation of errors in the quantum simulation. These regions of sharp change in the Floquet operator, referred to as structural instability regions, appear typically at intermediate Trotter step sizes and in the weakly-interacting regime, and are thus complementary to recently revealed quantum chaotic regimes of the Trotterized evolution (Sieberer et al., npj Quantum Information 5, 1 (2019)). We characterize these structural instability regimes in $p$-spin models, transverse-field Ising models with all-to-all $p$-body interactions, and analytically predict their occurrence based on unitary perturbation theory. We further show that the effective Hamiltonian associated with the Trotter decomposition of the unitary time-evolution operator, when the Trotter-step size is chosen to be in the structural instability region, is very different from the target Hamiltonian, which explains the large errors that can occur in the simulation in the regions of instability. These results have implications for the reliability of near-term gate-based quantum simulators, and reveal an important interplay between errors and the physical properties of the system being simulated.


I. INTRODUCTION
A primary application of quantum computers is simulation of quantum many-body systems that are classically intractable. Quantum simulation has applications in a wide variety of fields such as quantum chemistry [3,4], high-energy physics [5,6], condensed-matter physics [7,8] and quantum machine learning [9,10]. The devices of the current era, characterized by intermediate system size (∼ 100 qubits) and lack of full fault-tolerant error correction, are referred to as noisy intermediate scalequantum (NISQ) devices. The goal of NISQ-era devices is to perform less-demanding tasks than required for universal quantum computer, but ones that can still surpass the capability of the classical computers [11,12]. Recently, there have been a number of studies about achieving quantum advantage in NISQ devices in the context of quantum simulation [7,[13][14][15][16], optimization [17] and sampling [18].
To implement quantum simulation in a gate-based architecture, a common approach is to approximate the time-evolution operator generated by target simulation Hamiltonian using the Trotter-Suzuki decomposition [19][20][21][22][23]. Consider the target time-independent * kchinni@unm.edu Hamiltonian, H tar , given by and suppose that the time evolution operator associated with each of the individual terms in the target Hamiltonian, {e −iH1t , e −iH2t }, can be implemented. Then, the target time-evolution operator U tar = exp (−iH tar t) can be implemented through the first-order Trotter-Suzuki decomposition given by with the Trotter-step size given by τ = t/n. In the limit n → ∞, the Trotterized unitary in Eq. (2) becomes identical to U tar . In practice, n is a finite number and leads to errors in the overall simulation, which are bounded by where ||.|| is the spectral norm.
Here, n = O [H 1 , H 2 ] t 2 / is chosen so that the simulation has an overall accuracy [21]. Henceforth, we will refer to the errors resulting from the Trotter-Suzuki decomposition as Trotter errors.
In recent works, the errors resulting from Trotter approximation have been analyzed by associating the unitary resulting from Trotter-Suzuki decomposition U trot with the Floquet operator of a time-dependent arXiv:2110.03568v2 [quant-ph] 8 Apr 2022 periodically-"kicked" Hamiltonian H δ (t) [1,2]. This analysis revealed the existence of a regime in which dynamics of low-order observables given by the kicked Hamiltonian yields an accurate approximation to that of the target time-independent Hamiltonian. However, this regime breaks down for large enough Trotter step size, where errors are not controlled anymore, and the Floquet dynamics significantly deviate from the target dynamics. This uncontrolled error regime was attributed to the fact that kicked system develops quantum-chaotic features that are absent in the target Hamiltonian.
The existence of such nontrivial crossover in the Trotter error behavior reveals that the physical properties of the simulated Hamiltonian obtained via mapping to an effective time-dependent Floquet system can determine behavior of quantum-simulator errors. A related recent study, involving the simulation of adiabatic dynamics, characterized the emergence of quantum simulation errors to the closing of the spectral gap in the simulated Hamiltonian [22,23]. These works demonstrate that the physics associated with Floquet dynamics can affect how strongly errors saturate the Trotter-error bounds.
In this work we identify a new physical mechanism whereby the Floquet dynamics associated with Trotter-Suzuki decomposition lead to large errors in quantum simulation. Specifically, we study quantum simulation of p−spin models [24][25][26][27][28], which describe long-range interacting systems with multiple-body interactions, or equivalently, the dynamics of a large collective spin with nonlinear evolution to the p th power. Such models are meanfield models, where the mean-field dynamics exactly describes the thermodynamic/classical limit [27]. The classical dynamics on the phase space thus informs us about the potential mechanisms for errors in a quantum simulator. While the target mean-field collective spin models we will study here are classically simulable, they represent useful models for benchmarking quantum simulators. That is, these models can be used to analyze the reliability of the output of a quantum simulator by comparing its output with the classically obtained result.
We show that, even in the absence of chaos, there are parameter regimes where bifurcations arise in the meanfield dynamics associated with U trot that can cause a qualitative change in the structure of the classical phase space. These bifurcations appear only in the classical limit of the Floquet map corresponding to U trot and are absent in the ideal target unitary U tar , signaling the existence of regimes where the target and simulated dynamics strongly diverge. These dynamical instabilities can lead to large errors in the expected value of various observables. We name these parameter regions "structural instabilities" of the simulated unitary U trot , and characterize them in detail for p-spin models.
As the mean-field dynamics are equivalent to the classical dynamics in the thermodynamic limit, we study the structure of the classical phase space associated with U trot (τ ). We show that the structure of phase space undergoes significant changes accommodated by multiple bifurcations in all the instability regions as τ is varied by a small amount. A precursor of this phenomenon was previously noted in Ref. [29], where the route from integrability to chaos was studied for periodically kicked pspin systems. Here, we show that that this phenomenon does not require us to be in the thermodynamic limit, and use unitary perturbation theory to show that in the regime of structural instability, the eigenstates of the Trotterized unitary substantially deviate from those of target-unitary eigenstates, leading to an effective Hamiltonian H eff associated with U trot = e −iH eff τ , that is very different from the target Hamiltonian. Crucially, our analysis shows that structural instabilities appear at Trotter step sizes which are typically smaller than those required for the system to transition to chaos, making this novel regime particularly relevant to understanding the behavior of errors in near-term quantum simulators.
The remainder of the manuscript is organized as follows. In Sec. II, we discuss the mapping between quantum simulation via Trotter-Suzuki decomposition and Floquet dynamics, and discuss the intuition behind the emergence of quantum simulation errors from a classical perspective. We then introduce p-spin models and analyze two different regions of the parameter space that result in large errors as a result of Trotterization for a general Hamiltonian: the chaotic region and structural instability region caused by bifurcations. Then we specifically analyze how the structural instability regions lead to large errors in the Trotter approximation of the timeevolution operator for p-spin models. Following this, in Sec. III, we study the regions of structural instabilities through the use of unitary perturbation theory, which allows us to analytically predict the behavior of errors in long-time averaged magnetization. In Sec. IV, we construct the effective Hamiltonian associated with a general structural instability region in the p-spin Hamiltonian, which explains the presence of large errors in various observables and the significant changes taking place on the classical phase space in these regions. Finally, we show that the effective Hamiltonian construction accurately describes the appearance of unstable fixed points, and in particular predicts the growth-rate of the out-oftime-order correlator (OTOC).

II. TROTTER ERRORS IN QUANTUM SIMULATION OF p-SPIN MODELS
In this section we describe the connection between the Trotter-Suzuki decomposition of the target unitary evolution operator and Floquet dynamics, and discuss the implications of Trotterization in the classical limit in terms of Hamiltonian flows and area-preserving maps. We then introduce p-spin models and their kicked counterparts, and use these models to illustrate the different possible scenarios that could lead to the proliferation of Trotter errors.
A. Trotterized evolution and kicked systems: quantum and classical In a quantum simulator, the time evolution map generated by the desired target Hamiltonian H tar can be implemented through the use of the Trotterized unitary given in Eq. (2). This unitary map can be analyzed by identifying U δ (τ ) as the time evolution operator generated by a time-dependent periodic Hamiltonian H δ (t), which takes the the form The corresponding unitary evolution for one time period is given by the Floquet operator, where T is the time-ordering operator. Due to the impulse-like driving present in Eq. (4), these are sometimes called "kicked" systems. Since H tar and H δ (τ ) are different Hamiltonians except in the limit τ → 0, the Trotterized simulation is expected to be different from the ideal Hamiltonian evolution, resulting in simulation errors for a finite-sized Trotter step size τ , as shown in Eq. (3).
A way of studying the physical mechanisms behind these errors is to rewrite Eq. (4) as where g τ (t) = τ f τ (t) − 1. Since H tar = H 1 + H 2 , the evolution of H δ (t) corresponds to that of H tar under the action of an additional time-dependent perturbation, whose action is weak for small τ , as can be deduced from Eq.
. If H tar corresponds to an integrable Hamiltonian, the inclusion of a time-dependent perturbation is expected to break such integrability. The proliferation of errors in the quantum simulation of an integrable system that becomes chaotic as a result of Trotterization was discussed in [1] for the case of the quantum kicked top. Furthermore, a similar behavior was found even when considering quantum Hamiltonians with many-body quantum chaos, where the classical limit is not clear cut [2]. From this picture, we can expect the transition from regularity to chaos in the Trotterized unitary to be a fairly general phenomenon.
In contrast, here we will focus on a regime where the perturbation is not strong enough to make the system chaotic, but nonetheless has the potential to make the perturbed dynamics significantly different from the target dynamics. In classical Hamiltonian systems, the Kolmogorov-Arnold-Moser (KAM) theorem guarantees that for sufficiently small τ , the regularity of the original Hamiltonian is preserved [30], and thus one expects that the perturbed evolution to be still a good approximation for the ideal dynamics. However, there are situations like the presence of resonant regular orbits, which fall outside the validity of the KAM theorem [31] and in which even small perturbations can have a big effect in the dynamics of the system. These changes are signaled by the emergence of new fixed points through bifurcations, among other mechanisms. As this happens, the emergence of instabilities and the development of significant changes in the phase space structures also have the potential to substantially modify the dynamics of the effective system H δ from that of the ideal target Hamiltonian H tar .
As we will show, for some models these features can persist in quantum dynamics far from the classical limit, and these instabilities determine parameter regimes of high Trotter errors. Moreover, since these features appear at smaller perturbation strengths (before the transition to chaos), they would affect the Trotterized evolution at smaller values of the Trotter-step size τ , making them particularly relevant for quantum simulation. A way to identify these high Trotter-error regions in the quantum regime is to determine the regions where the eigenstates of U δ are very different from U tar . This follows from the correspondence principle. The quasi-probability distribution (e.g., Husimi distribution) of the eigenstates associated with the time-evolution operator is expected to have a significant overlap with the corresponding phasespace trajectories, when both of them are plotted on the classical phase space [32]. Thus, the eigenstates of U δ and U tar will be very different whenever the corresponding classical trajectories are different. In the following we will explore this phenomenon in the quantum simulation of long-range interacting spin models, whose mean-field limit is equivalent to the thermodynamic/classical limit, , with α = π, π/2, respectively. These represent regions of structural instability that are not associated with chaos as the largest Lyapunov exponent corresponding to these regions is zero.
by analyzing the eigenstates of U δ and U tar . In addition, as we will discuss, some of the underlying physical mechanisms behind these phenomena are much more general, and are expected to appear in other types of quantum many-body systems.

B. p-spin models
The family of magnetic models usually referred to as p-spin models describes a collection of N spin-1/2 particles on a fully connected graph, interacting through pbody Ising-like coupling in the presence of an external transverse magnetic field. This collection of particles experience two different orderings, a paramagnetic phase induced by an external homogeneous magnetic field, and a ferromagnetic phase induced by a p-body Ising-like interaction. The Hamiltonian for this family of models is given by where h is the strength of the external field, γ is the strength of the p-body Ising-like interaction, and J µ = µ with µ = x, y, z, are collective spin operators. The interaction term has been normalized with this particular choice of p and N to make the equations of motion have a universal form for all p in the mean-field limit and energy extensive [28]. In this Hamiltonian, the total angular momentum is conserved, [H, J 2 ] = 0, constraining the dynamics to the symmetric subspace, and the system is also invariant under the action of the parity operator, Π = e iπJz , for even p spin models. For the remainder of this manuscript, we consider the following single-parameter version of the p−spin Hamiltonian, which is equivalent to Eq. (7) upon rescaling of the energy, and where s is constrained to be in the range 0 ≤ s ≤ 1, interpolating between pure paramagnetic ordering and a pure ferromagnetic ordering. For p = 2, the above Hamiltonian reduces to the Curie-Weiss model [33], which is a special instance of the Lipkin-Meshov-Glick (LMG) model [34][35][36]. This family of models has been extensively studied in the context of quantum annealing [24,25], where a classification based on the properties of the ground state phase transition (GSQPT) was constructed [26,27]. Such classification splits the family of models in two classes, for p = 2 the GSQPT is continuous and second order, for p > 2 the GSQPT is first order and discontinuous. Furthermore, in the context of dynamical criticality, this family of models exhibit dynamical quantum phase transitions [28].
The dynamics of the p-spin Hamiltonians in the meanfield limit can be obtained by neglecting the fluctuations AB ≈ A B in the Heisenberg equations of motion resulting in the coupled differential equations of the form [28,37] with {X, Y, Z} = lim J→∞ 1 J { J x , J y , J z }, describing the motion of a classical "top". The mean-field limit coincides with the classical limit, eff = N −1 → 0 or equivalently N → ∞, in the case of p-spin models [27]. Due to the conservation of angular momentum, the resulting dynamics of the top is constrained to a unit sphere X 2 +Y 2 +Z 2 = 1. This implies that all the p-spin models are integrable in the classical limit as they correspond to autonomous systems with one degree of freedom.
The first-order Trotterized unitary map U trot = (U δ (τ )) n of the Hamiltonian evolution in Eq. (8) is generated by the time-evolution of the corresponding kicked model, We refer to this family of models as the kicked p-spin models [29], whose Floquet operator is given by The equations of motion are then obtained using the map Note that these models also conserve the angular momentum, [H δ , J 2 ] = 0 similar to the case of H tar implying that the dynamics of H δ is constrained to a unit sphere defined by X 2 i + Y 2 i + Z 2 i = 1 in the classical limit. The classical equations of motion are given by the following map, where α = −(1 − s)τ and k = −sτ . A comparison of targeted Hamiltonian flow, Eq. (9), at periodic intervals with the area-preserving map generated by the Trotterized kicked Hamiltonian, Eq. (10), indicates where errors may occur in the quantum simulation.
The kicked models can exhibit chaos since the energy of the top is no longer conserved. For small values of τ , however, where the kicked-systems in Eq. (11) very closely approximate the evolution of corresponding Hamiltonians in Eq. (8), the dynamics is regular. As the value of τ is increased, all kicked p-spin models become chaotic provided s is not close to zero or one; the models with larger value of s develop chaos at smaller values of τ . The kicked p = 2 model (Haake's kicked-top [38]) develops chaos due to the period-doubling cascade and reaches the regime of strong chaotic trajectories faster (as a function of the coefficient of non-linear term in the Hamiltonian) than other p-spin models [29]. Higher-order kicked p-spin models, on the other hand, develop chaos due to instability of higher-period orbits present on the XY plane [29].
In order to analyze the emergence of quantum simulation errors in the different parameter regimes of the Trotterized evolution, we perform a systematic comparison between the basis of eigenstates of U tar and U trot .
In Figs. 1a and 1b, we plot the average dissimilarity between both sets of eigenstates for the p-spin model with p = 2 and p = 4, correspondingly. This quantity measures the difference between two sets of eigenvectors and is defined as where IPR is the average inverse participation ratio (IPR) of the eigenstates of U tar , denoted by {|φ where d is the dimension of the Hilbert space. Also, is the average IPR in the circular orthogonal ensemble (COE) [1]. We expect the average IPR between the eigenstates of U tar and U trot to be equal to IPR COE when the dynamics of U trot becomes fully chaotic, as COE is the appropriate ensemble due to the symmetries present in the kicked p-spin Hamiltonian [29]. The dissimilarity defined in Eq. (13) ranges between 0, when the eigenstates are identical to the reference basis, and d−1 d (up to normalization) for the case when the eigenstates are completely delocalized in the reference basis.
The dissimilarity of the eigenvectors is shown in the heat map in Fig. 1a and 1b. At larger values of s, the dissimilarity arises mainly due to the chaos present in the kicked p-spin models. This is indicated by the region of parameter space with a positive classical Lyapunov exponent (see in Fig. 1c and 1d for p = 2 and p = 4, correspondingly. For more details on calculation of Lyapunov exponents, see Appendix A and [29,39]). In contrast, the large dissimilarity present at smaller values of s on the heat map cannot be attributed to chaos since the associated Lyapunov exponents are zero. We identify these parameter regimes as structural instabilities of the operator U trot , since small changes in the parameter space (s, τ ) induce substantial changes in the nature of the eigenstates of the operator. In the following section, we will derive the precise location of these structural instabilities, and relate them to the proliferation of errors in quantum simulation.

III. TROTTER ERRORS IN OBSERVABLES DUE TO STRUCTURAL INSTABILITIES
In this section we predict the location of regions of structural instability identified in Fig. 1 and analyze their effect on quantum simulation errors. In order to x ≡ U (0) U . Hence the eigenvectors of the Floquet operator U δ are expected to be close to the eigenstates of J z with a small correction first order in s. Employing unitary perturbation theory [40], we find Here m is the first-order eigenphase correction and |φ (1) m is the corresponding first-order correction to the eigenstate. This can be expressed as Note that the ideal unitary map, U (0) , has degenerate eigenvalues when e i(1−s)τ (m−m ) = 1 (equivalently at τ = r 2π (1−s)(m−m ) for some positive nonzero integer, r). We denote the degenerate points by τ * p,m−m corresponding to the kicked p-spin model with a particular solution for m − m . The neighborhood of these degenerate points can be divided into two subregions. First, in the zone surrounding the degenerate point, the firstorder correction to the eigenphase is on the order of gap between the eigenphases, and the nondegenerate perturbation theory is not valid in this subregion. We refer to this subregion as the "immediate vicinity" of the degenerate point. In the subregion beyond this immediate vicinity ("outer vicinity"), the non-degenerate perturbation theory is valid, and the expression in Eq. (16) predicts a large correction to the eigenstates whenever (J p x ) m ,m = 0. We refer to the whole region consisting of immediate and outer vicinity subregions surrounding the degenerate points as the structural instability regions whenever (J p x ) m ,m = 0. We will refer to τ * p,m−m as the structural instability points, which are the central points of these regions, and denote width of the regions by w.
For concreteness, consider the area around τ = τ * 2,2 = r π (1−s) , that is, corresponding to p = 2 and m − m = 2. In this case, we have that Thus, by virtue of Eq. (16), we obtain a substantial change in the eigenstates of the unitary evolution operator, indicating the existence of an instability region. At τ = τ * 2,2 , the Floquet operator is given by U δ (τ * 2,2 ) = e iπJz e i sτ * 2J J 2 x whose eigenstates are . This can be understood from the fact that the two terms in U δ (τ * 2,2 ) commute and therefore have a common set of eigenvectors given by the parity respecting eigenstates of J 2 x . Hence, as the Trotter-step size is varied in the vicinity of τ * 2,2 = r π (1−s) , the eigenstates of U δ change rapidly from those of J z eigenstates at the left edge of the structural instability region (τ ∼ τ * 2,2 − w 2 ) to those of parity respecting eigenstates of J 2 x at the central point of the instability region before changing back to those of J z eigenstates as τ is increased further towards the other edge of the instability region (τ ∼ τ * 2,2 + w 2 ). The above argument holds for all even p-spin models, so all these models have structural instability regions in the vicinity of the curve (1 − s)τ * p,2 = rπ.
In general, for a given p, J p x has nonzero matrix elements in the J z -basis only in alternating diagonal bands up to offset p. This results in structural instabilities in the vicinity of the curves given by .., 2(1)} for even (odd) values of p. We corroborate this in Fig. 1a and 1c, where we show the curves (1 − s)τ * = π (black-dashed line) and (1 − s)τ * = π 2 (black-dotted line) overlap with the dissimilarity region present on the heat map. Note that the number of structural instability regions increases with p in a given range of τ as the number of choices for m − m increase with p.
The significant change in the eigenstates of U δ that we find in the structural instability regions implies that U δ becomes very different from U tar , which can lead to large errors in a Trotterized quantum simulation algorithm. For concreteness, we focus on the simulation of the long-time average of the collective spin observables J i , defined by where J i (lτ ) = (U † ) l J i (U ) l , with U being the map associated with the time-evolution operator for t = τ and i = {x, y, z}. We analyze the error in J z given by where J z tar and J z trot are the time-averaged magnetizations obtained under the target unitary and the Trotterized unitary correspondingly. As studied in [1,41], quantum simulation is expected to be robust to imperfections in the nonchaotic regime of H δ for expectation values of macroscopic observables that are not sensitively dependent on the full state of the system compared to quantities such as the fidelity of preparing a target state. However, the Trotterization of p-spin models leads to a large region of error, even in the simulation of a macroscopic quantity such as J z . In Fig. 2 we have plotted the error E ∞ z (τ ) for an initial spin coherent state |Θ 0 = π 2 , Φ 0 = 0 at s = 0.1 as a function of the Trotterstep size for p = 2, 3 and 4 in parts (a),(b) and (c) respectively. As expected, the errors increase in structural instability regions in the vicinity of τ * 2,2 = π 1−s for p = 2, cases (all cases shown by the vertical dashed lines). Note that at every structural instability region shown in Fig. 2, the errors first increase rapidly, then decrease in some intermediate region before increasing again resulting in seemingly two separate error peaks. The presence and location of the error dip between the error peaks is dependent on the initial condition, and will be analyzed in the next section.
An important consequence of this analysis is that the structural instability regions associated with higher values of p occur at smaller values of τ for a given value of s, since τ * p,m−m = r The spectral gap given by the denominator term in Eq. (16) determines the width of the error regions. From the Taylor expansion of the inverse spectral gap around the structural instability points τ = τ * , we have, The parabolic form of the denominator centered at τ * with the width given by The error in the long-time averaged magnetization along the z-axis to first-order in s for initial spin-coherent states |Ψ (0) = |Θ, Φ can be expressed analytically. We find, where A r1,r2 = −J + r 1 |A| − J + r 2 for an operator A. The error expression for an arbitrary initial state is shown in the Appendix B. Equation (20) predicts an error peak, captured by the cotangent term in the outer vicinity region, for each value of q in the summation, which corresponds to having an error peak at every structural instability region. This result from perturbation theory also agrees very well with the numerically obtained curves as shown by the dotted lines in Fig. 2 for the spin-coherent states centered at |Θ = π 2 , Φ = 0 except in the immediate vicinity of the degenerate point, where the errors have inverted-triangular shape (analytic prediction in this region is not shown in Fig. 2 because these predictions diverge here, and we expect this because the non-degenerate perturbation theory is not valid in this region). This behavior holds true for most of the other spin-coherent states at small values of s.

IV. EFFECTIVE HAMILTONIAN AND EMERGENT SYMMETRIES
As seen in the previous section, Trotterization of the p-spin models leads to large errors in the vicinity of certain parameter regimes corresponding to so-called structural instabilities. This was understood in the classical limit to be a consequence of bifuractions occurring in the area-preserving map of the Trotterized model that radically shift the structure of phase space, and which manifests at the quantum level in a Floquet operator whose eigenstates are very different from those of the target p-spin Hamiltonian. This difference in the structure of eigenstates can be further elucidated through the construction of an effective Hamiltonian associated with the Trotterized unitary. For small values of s, we have U tar ∼ e iκJz , and the evolution is essentially precession of states around the z-axis. These precessions are well approximated by the Trotterized evolution U trot = U δ (τ ) n away from the structural stabilities. However, near these instabilities, the phase space of the Trotterized evolution undergoes major structural changes leading to a evolution very different from precessions of the state around the z axis.
For instance, consider p = 2, the LMG Currie model, whose phase space in the classical limit associated with U tar is shown in Fig. 3(a). U δ (τ ) is shown in Fig. 3(b-d) at s = 0.1 in the neighborhood of the instability that is present at τ * 2,2 = π 1−s ≈ 3.49. We understand this structure of phase space by analyzing the form of U 2 δ , For s, ∆τ τ * 2,2 +∆τ 1, the unitary map can be written as meaning that the dynamics for p = 2 every two time steps can be described by the effective Hamiltonian H (2,2) eff . For ∆τ > 0, the effective Hamiltonian is in fact the LMG Hamiltonian of the form shown in Eq. (8) with an additional overall multiplicative factor and an effective s parameter given by leading to s eff is always greater than 0.5 in the region of the structural instability. This implies that the Trotter approximation of the unitary map with the original Hamiltonian having a small s value (i.e., being in the paramagnetic phase) leads to simulation of the dynamics of the same model but with a large value of s, (i.e., corresponding to the ferromagnetic phase, up to every alternate step). This paradoxical effect can also be seen in the classical phase space shown in 3(b-d) for ∆τ > 0, where the trajectories change from precessions around the z-axis for τ > π 1−s (1 + s 1−2s ) (s eff > 0.5) to precessions around x-axis at τ = τ * 2,2 = π 1−s (s eff = 1) as τ is decreased (from right to left in Fig. 3).Notice that the time interval required to trace out the ferromagnetic phase dynamics using the Trotterized unitary (effective Hamiltonian) is 1 s eff slower compared to the time required for ideal LMG Hamiltonian (p = 2 in Eq. (8)) with s = s eff to follow the same effective dynamics. This can be traced back to the difference in the overall multiplicative factor, H In the mean-field picture, this process is accommodated by a period-doubling bifurcation at τ = π 1−s (1 +  dynamics trace out the two parity-broken trajectories, simultaneously as shown in Fig. 3e, where the phase-space trajectories are plotted as a function of angular coordinates θ and φ. For a given initial condition the Trotterized dynamics trace out one lobe (red-colored trajectories) in all the odd steps of the evolution and the other lobe in all the even steps (purple-colored trajectories) of the evolution. In this way, the Trotterized trajectory jumps between two separate parity-broken trajectories of the LMG Hamiltonian tracing out the ideal LMG dynamics with s = s eff only every alternate step. On the other hand, for the initial conditions associated with the parity conserving trajectory of the LMG Hamiltonian, every step of the Trotterized unitary traces out the ideal LMG trajectory with s = s eff (black-colored trajectories). Similar phenomenon takes place in the instability region for ∆τ < 0 except that the bifurcation now takes place in the Z < 0 hemisphere. Performing a similar analysis of the structural instability region present at τ * p,2 = π 1−s for the even p-spin models shows that the Trotterized evolution also results in simulation of the ferromagnetic phase of the corresponding p-spin Hamiltonian (or the Hamiltonian with a relative negative sign for ∆τ < 0) even though the target evolution is associated with the paramagnetic-phase dynamics.
The Trotterized unitary dynamics in the vicinity of other instabilities can differ even more substantially from the target, ideal dynamics. For example, consider the instability at τ * 4,4 = π 2(1−s) for p = 4. The effective Hamiltonian can be derived in a similar manner as described above, yielding whose phase-space trajectories undergo two different 1to-4 bifurcations, as can be seen in Fig. 3(g-i). The change in phase space structure results in trajectories that are very different from those associated with the target dynamics shown in Fig. 3f. Similar to the case of p = 2, the parity-broken trajectories of H (4,4) eff located on the phase-space are traced out by the Trotterized dynamics only every fourth step. In the intermediate steps, the Trotterized dynamics leads to jumps between various parity-broken lobes present on the phase space as shown in Fig. 3j, where the red, purple, blue and green colored trajectories represent every first, second, third and fourth step respectively.
More generally, the dynamics in the vicinity of instability at (1−s)τ * p,q = r 2π q has a 1-to-q bifurcation present on the classical phase of the Trotterized unitary and can be understood by analyzing U δ (τ + ∆τ ) q . We see where The Hamiltonian in Eq. (29) simplifies further when the Hamiltonian has parity symmetry, which is the case for all even p-spin models, as the ( q 2 + k) th term in the summation becomes identical to the k th term, reducing the number of terms in the summation from q to q 2 . This effective Hamiltonian captures the dynamics of every q th step of the Trotterized unitary in the vicinity of 1-to-q bifurcation. The associated phase-space of the effective Hamiltonian is invariant around the z-axis under 2π q rotation since it commutes with e i 2π q Jz : [H (p,q) eff , e i 2π q Jz ] = 0. This is an emergent symmetry that appears in the structural instability region. We also want to point out that even though bifurcations facilitate the structural changes in the instability regions, not all bifurcations lead to such sharp changes in phase space. Only the subset of bifurcation points, identified here as "significant" bifurcations, lead to extensive changes in the structure of the eigenstates and result in large Trotter errors. These significant bifurcations appear only in the structural instability regions.
In summary, the effective Hamiltonian formulation explains the presence of large error peaks in the simulated time-averaged magnetization. At certain Trotter step sizes that correspond to structural instabilities, the Trotterized unitary evolution operator simulates a very different Hamiltonian from the target Hamiltonian. From a mean-field perspective, the major structural changes that take place inside the regions of structural instability always correspond with significant bifurcations on the classical phase space, which lead to creation of periodic hyperbolic points (hyperbolic fixed points of higher period), which are absent on the phase space associated with the target Hamiltonian H tar ∼ κJ z .

V. INFORMATION SCRAMBLING INSIDE STRUCTURAL INSTABILITY REGIONS
As mentioned in the previous section, one of the signatures of the structural instability regions in the Trotterized unitary is the emergence of or unstable (hyperbolic) periodic fixed points. Thus, in the classical limit, the trajectories in the vicinity of this point are expected to show exponential divergence when the dynamics is observed at the appropriate stroboscopic times, i.e., every q steps in the region around τ * p,q . Recently, it has been shown that the presence of hyperbolic points can lead to information scrambling deep inside the quantum regime, a phenomenon that was dubbed saddle point scrambling [42,43]. Note that this saddlepoint scrambling is not indicative of chaotic behavior as mentioned in [42,43], rather its origin lies in the exponential divergence of the trajectories that happens only in the localized region around the separatrix, which includes the hyperbolic fixed point. Here, we show that the effective Hamiltonian constructed in the previous section correctly identifies the presence of unstable fixed points and allows us to characterize the saddle point scrambling emerging in the simulator. Notice that such scrambling exists naturally in the dynamics of the ferromagnetic phase of the ideal p-spin models, as the existence of a critical (for p = 2) or a bifurcation (for p > 2) point is always accompanied by the emergence of an saddle point. This is easily understood as a pitchfork or a saddle-node bifurcation of the corresponding classical equations of motion. However, this type of scrambling is absent in the the paramagnetic phase, which is the target dynamics being simulated.
The smoking gun of scrambling behavior is the exponential operator growth of certain correlation functions [44]. In particular, we characterize the saddle-point scrambling in the quantum simulation using the short time growth of the "infinite temperature" square commutator where d is the dimension of the Hilbert space, the operators V (t) and W (0) are chosen so that they commute at the initial time, and V (t) is the Heisenberg evolution of V (0). In this work, we choose V (0) = W = J z , and study the short time growth of the square commutator In the presence of an instability, be it a saddle point or a hyperbolic periodic point, the above quantity is expected to grow like c(t) ∼ e λ saddle t , where λ saddle is the associated growth rate at the saddle point. The exponential growth of the OTOC is seen in Fig. 4a and 4c for the p = 2 system with Trotter step size in the vicinity of τ *

2,2
and p = 4 system with Trotter step size around τ * 4,4 , respectively. The exponents associated with the saddle point derived from these numerics are plotted in Fig. 4b  and 4d as a function of the Trotter step size with the red dotted line.
The analytic expressions for λ (p,q) saddle associated with a given p and 1-to-q bifurcation can be obtained by solving for the eigenvalues of the Jacobian matrix associated with the linearized classical flow around the appropriate unstable fixed point of the effective Hamiltonian, which are labelled by blue solid lines in Figs. 4b and 4d. In the system with p = 2 around τ * 2,2 , i.e the 1-to-2 bifurcation, one finds where M (4,4) + is the largest eigenvalue of where and X 2 sd = Y 2 sd are the Cartesian coordinates of the unstable point of the classical flow associated with the effective Hamiltonian at this instability for the p = 4 system. We give the explicit expression of Z sd along with details on the derivation of the exponents in appendix C. In in Fig. 4(b) and Fig. 4(d) we compare the exponents extracted from the numerical calculation of the OTOC (red stars), for a system size of N = 128, with the analytical expressions obtained with the mean-filed limit of the effective Hamiltonian (solid blue). Notice the good agreement despite the small system size used in the simulation. Thus, as first observed in [42] the short time growth of the OTOC has the form e λ saddle t .

VI. CONCLUSIONS AND OUTLOOK
In this work we have identified a new mechanism leading to the proliferation of errors in a quantum simulator when the algorithm employs the Trotter-Suzuki decomposition. In the mean-field limit, these regions of structural instability are characterized by multiple bifurcations leading to rapid global changes in the structure of the phase space as the Trotter-step size is varied slightly.
The effects of these bifurcations can be seen in smaller systems, with N far from the thermodynamic limit. A method to identify the structural instability regions in these smaller systems is to seek the regions where the eigenstates of the Trotterized unitary differ significantly from the eigenstates of the target unitary. From the correspondence principle, the quasiprobability distribution is expected to overlap with the trajectories of the phase space in the classical limit, so rapid changes in the structure of the phase space are reflected in the modifications of the eigenstates associated with the Trotterized unitary. The sudden changes in the structure of the eigenstates of the Trotterized unitary result in Floquet dynamics that is very different from the targeted evolution of a given state, resulting in large errors in various observables. In this work, we provide analytic expressions for these high-error regions computed using unitary perturbation theory, for the case of p-spin models. We showed that inside the structural instability regions the effective Hamiltonian, which is the generator of the Trotterized evolution, is very different from the target Hamiltonian, providing further justification for the presence of large errors in these regions. The effective Hamiltonian reveals the emergence of new unstable fixed points in the structural instability regions, indicating the presence of saddle point scrambling in the simulator, as manifested by the exponential growth of the OTOCs in these regions.
An important conclusion of the perturbation theory analysis of Sec. III is that structural instability regions will appear at smaller values of the Trotter step size τ as the value of p increases. This indicates that, in general, gate-based quantum simulations of p-body interactions (beyond the usual two-body case of p = 2) is likely to lead to more parameter regimes where Trotter errors proliferate. While this is certainly true in the case of all-to-all interactions analyzed here, the study of its manifestation in systems with different interactions (i.e., with long, but finite, interaction range) will require further study. We conjecture that the structural instabilities studied here in the mean-field case are more universal, and carry over to Floquet states with more general many-body interactions. Note that models with all-to-all interaction graphs are often seen to correctly capture the physics of sys-tems with more complex but still long-range interactions. For example, for the case of p = 2, the mean-field phenomenology associated with the LMG model informs the DQPT behavior of finite-range interacting systems described by x for α 2 (α = 0 corresponds to the LMG model) [14,45]. These results strongly suggest that mean-field models can give insight into many-body behavior (particularly, the nonequilibrium dynamics of physical observables), and the notion of structural instabilities will very likely be present in finite-range models with nonzero values of α for the p = 2 case. Further study of the structural instability regions in the finite-range interaction case, and the extension of these models for p ≥ 3 is left for future work.
Beyond their implication in the proliferation of errors in quantum simulation, the regions of structural instabilities are a manifestation of fundamental effects in the nonequilibrium dynamics of the Floquet system. Particularly, we find that the dynamics of the driven system in Eq. (11), whose Hamiltonian is periodic with period T , H(t + T ) = H(t), shows signatures of Floquet time crystal behavior in a region of structural instability. A Floquet time crystal is an out-of-equilibrium phase of matter that breaks discrete time-translation symmetry [46][47][48][49]. These time-crystal phases in the kicked p-spin models can be studied with the help of area-preserving maps associated with the mean-field limit of the Floquet system. In particular, the existence of periodic elliptic points, which can be observed clearly in the phase space portraits of Fig. 3, will lead to the dynamics where the initial state, prepared close enough to elliptic points, will periodically return to the initial configuration after q time steps, with q ≥ 2. Physically, this means that the response of the system will have a periodicity of qT , instead of T , thus breaking the discrete time-translation symmetry of the Hamiltonian. The connection between structural instabilities and subharmonic response was also seen in [50], where the author studied the emergence of robust subharmonic response in a spin chain with short-range (nearest-neighbor) two-body interactions in a regime which roughly coincides with the choice of τ = τ * 2,2 in the parametrization used in our work. This result, together with the general connection of the structural instability regions with discrete time crystals, indicates that this source of Trotter errors is not an isolated phenomenon happening only in long-range interacting spin models. A comprehensive study of all the Floquet time crystal phases present in the kicked p-spin system is part of an ongoing work and will be presented in an incoming publication [51].

ACKNOWLEDGMENTS
We would like to thank Poul Jessen and Kevin Kuper for insightful discussions, and Changhao Yi for his insights in Trotter formulas and quantum simulation. This work was supported by the U.S. National Science Foun-dation under grant numbers PHY-1820679 and PHY-2011582. This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator.

Appendix A: Calculation of Lyapunov exponents
Throughout the main text, we use the term Lyapunov exponent to describe the rate of exponential growth of a classical variable X. That is, up to some threshold time t th , this classical variable evolves according to X(t) = X(0)e Λt , with Λ being the Lyapunov exponent.
We have presented results in this work for two different scenarios where some components of the mean-spin in the thermodynamic limit evolve following the dynamics described in the previous paragraph. First, we have analyzed the chaotic instability of the delta-kicked p-spin as illustrated in Fig. 1c,d. Then, we analyze the short time evolution of the square commutator in Sec. V, whose physical origin can be traced to the exponential growth of the unstable separatrix branches of saddle points of the phase-space flow. Notice that we provide a more in-depth analysis of the instability around the saddle point in Appendix C.
Although, both scenarios represent exponential instabilities, they correspond to physically different situations. The chaotic instability is global, so any pair of points inside the chaotic region of the phase space will display exponential divergence of their distance at a rate given by the exponent. The ubiquity of this instability all over the chaotic region leads to the folding of trajectories at long times, and the decay of correlations, the later known as mixing. In fact, at long times and length scales, the motion inside the chaotic region resembles a diffusion process. This implies that the chaotic dynamics has a positive value of metric entropy, or Kolmogorov-Sinai entropy, which can be computed via Pesin's theorem [52]. For the instability of the saddle point, all of the above mentioned physical processes are absent. Furthermore, the exponential divergence of trajectories only occurs in a highly localized region of phase space including the immediate vicinity of the saddle point and the separatrix. For instance, see Ref. [42] for an in depth discussion on this matter.
The motivation behind Fig. 1 is two fold. On the one hand, we want to introduce the disimilarity, and recognize that this quantity identifies both types of instabilities. On the other hand, we also want to present a direct comparison with the Lyapunov exponent of the chaotic instability, allowing us to identify the region of parameter space (τ, s) in the similarity heat map whose origin is chaos.
In the following, we explain the details of calculating the Lyapunov exponent associated with the chaotic instability, as shown in Fig. 1. In the thermodynamic limit, the mean-spin evolves stroboscopically according to Eq. 12. For any point in phase space, the local dynamics of small increments in its vicinity is governed by the Jacobi matrix For a point inside the chaotic region, the exponential instability implies that the neighborhood of the point is getting exponentially strecthed in some of the principal directions of the Jacobi matrix and exponentially shrinked in the other principal directions. This takes place as one evolves the Jacobi matrix along the chaotic trajectory corresponding to the selected point. As such, the largest Lyapunov exponent can be computed as the largest eigenvalue of the Jacobi matrix evaluated along the trajectory, which follows from the celebrated Ergodic theorem of Oseledets [53,54].
For a map as the one in Eq. (12), Oseledet's ergodic theorem allows us to compute the largest Lyapunov exponent given by where n is the number of time steps, λ + is the largest eigenvalue of the matrix N m=1 M T (X m )M(X m ) and M(X m ) is the tangent map introduced before.
Naturally, the exponential growth of one of the eigenvalues leads to issues with the numerical computation of Eq. (A2). This can be avoided by looking at Eq. (A2) in a different basis other than the eigenbasis of the Jacobi matrix. This can be achieved via a QR decomposition, see for instance Ref. [55], which permits the numerical approximation of the asymptotic time limit in Eq. (A2).
A single point in the heat maps of Fig.(1)c,d is obtained by approximating Eq. (A2) via the QR method up to n = 10 6 time steps. This value is then averaged over 50 different initial points inside the chaotic region. This procedure is then repeated for a grid of points in the (τ, s) plane.

Appendix B: Long-time average of Jz
The long-time average of an operator A, assuming the time-evolution operator corresponding to one time-step has nondegenerate eigenphases, is given by where ρ (0) is the initial state and |φ r is the r th eigenstate of the system. The error in this observable due to Trotterization is given by where A ∞,id is the long-time average of A under the ideal Hamiltonian evolution, and A ∞,τ is the long-time average under Trotterized evolution. Assuming that the eigenstates of the system change under a perturbation |φ r → |φ , the expression for the long-time average to the first order is given by The error is then given by The above expression can be further simplified by expanding the summation in n and noticing the matrix elements of (J p x ) m,n are nonzero for n = m ± p, m ± p − 2, ..., m ± 0(1). Focusing on two particular terms with n = m ± p − q we obtain Manipulating the second term in above expression by first shifting the index of the second term in the above equation, m → m − (p − q), and then setting Re[z] = Re[z * ] in the second term results in the following expression for the error Focusing on the error in J z , the above expression further simplifies to the following Relabelling the index q → p − q results in the final expression, The equations of motion of the corresponding classical flow are given by this classical flow has two fixed points at the poles, X = Y = 0 and Z = ±1. Other fixed points satisfy Y = 0 and The range of values of ∆τ for which these new fixed points are real gives us the extent of the region of structural instability, from the above expression for X is easy to see that determines the width of the structural instability region.
Furthermore, it is not hard to see that the unstable point emerges as a consequence of the change in instability of one of the fixed points on the poles, depending on the sign of ∆τ . Hence, to compute the exponent we evaluate the Jacobi matrix on the poles and diagonalize it, finding that the two nonzero eigenvalues are given by where τ = ∆τ τ * 2,2 +∆τ . From the largest eigenvalue we obtain the expression for the value of the exponent λ saddle λ (2,2) saddle (∆τ ) = ∆τ (1 − s) sign(∆τ ) where we have included the appropriate prefactor accounting for the definition of the effective Hamiltonian.

4,2
In Sec. IV we mentioned that all the even kicked p-spin models have a structural instability centered at the same value as the p = 2 model, and the central value corresponds to the values of τ at which the period doubling bifurcation takes place. We now consider this structural instability region for the system with p = 4. The effective Hamiltonian is given by With the equations of motion for the associated classical flow given by where τ = ∆τ τ * 4,2 +∆τ . This flow has two fixed points on the poles, X = Y = 0 and Z = ±1, which are always stable. New fixed points can be found as the solution to dX dt = 0. Of particular interest to us are the ones whose coordinates satisfy Y = 0 and the x, z-coordinates are related via which leads to the cubic equation with the factor A = 1−s s τ . From the above relation between the x-and z-coordinates we see that X sd = ± 1 − Z 2 sd , then we can construct the tangent map and evaluate it in the saddle point. After doing so, we find that the exponent governing the exponential growth is given by λ (4,2) saddle (∆τ ) = s(τ * 4,2 + ∆τ ) 2 where again we have include the appropriate prefactor which accounts for the definition of the effective Hamiltonian.
c. The case of arbitrary even p Now that we have presented in detail the calculation of the exponents for the two models where explicit expressions can be obtained, let us briefly show some further insights which can be extracted from this approach, when we consider the cases of an arbitrary even p at the τ * p,2 structural instability region. Let us start with the equations of motion for the classical flow where now we have τ = δτ τ * p,2 +∆τ . This set of equations has two fixed points at the poles, X = Y = 0 and Z = ±1. The other fixed points satisfy Y = 0 and the x-and z-coordinates are related via Z = (1 − s)τ sX p−2 , and Z 2 = 1 − X 2 .
These two equations give rise to a algebraic equations of degree 2p − 2 for both coordinates. For instance, the one for the x-coordinate si given by Although a general solution to this algebraic equations is not available, we can use it to estimate the size of the structural instability region. Let us define the function S(X) = X 2p−2 − X 2p−4 + 1−s s 2 τ 2 , this function has extreme points at X = 0 and X ± = ± 2p−4 2p−2 , additionally we notice that S(1) = S(0) = 1−s s 2 τ 2 > 0 are always positive.
Then, the function S(X) will have at least one root in the interval X ∈ [0, 1] if there is a minimum in this interval and the function evaluated at this minimum is negative.
It is not difficult to see that X + is in fact a minimum of S(X). Thus, we want conditions on the parameters such that the inequality S(X + ) < 0 is satisfied. Solving for those conditions we find ∆τ ≤ sτ * p,2 F(p) 1 − s − sF(p) , with F(p) = (p − 1)(p − 2) p−2 − (p − 2) p−1 (p − 1) p−1 .
Thus we see that the size of the structural instability region around τ * p,2 decreases with increasing p, with the model with p = 2 having the most prominent one.
to investigate the width of the structural instability region, we look for the range of parameters such that the above equations has nontrivial real solutions. To do this we consider the function G(X) = 2X 2p−2 − X 2p−4 + 4 1−s s 2 τ 2 , the extreme points of this function are at X = 0 and X ± = p−2 2(p−1) . since G(0) and G(1) are positive, then if there is a minimum in the range X ∈ [0, 1] such that G(X) at this minimum is negative, then the algebraic equation has at least one nontrivial solution. In fact X + is a minimum, thus we want conditions on the function parameters such that the inequality G(X + ) < 0 is true. After solving for those conditions we find thus this region of structural instability shrinks with increasing p and is exponentially narrower that the region of structural instability around τ * p,2 .