The effect of chaos on the simulation of quantum critical phenomena in analog quantum simulators

We study how chaos, introduced by a weak perturbation, affects the reliability of the output of analog quantum simulation. As a toy model, we consider the Lipkin-Meshkov-Glick (LMG) model. Inspired by the semiclassical behavior of the order parameter in the thermodynamic limit, we propose a protocol to measure the quantum phase transition in the ground state and the dynamical quantum phase transition associated with quench dynamics. We show that the presence of a small time-dependent perturbation can render the dynamics of the system chaotic. We then show that the estimates of the critical points of these quantum phase transitions, obtained from the quantum simulation of its dynamics, are robust to the presence of this chaotic perturbation, while other aspects of the system, such as the mean magnetization are fragile, and therefore cannot be reliably extracted from this simulator. This can be understood in terms of the simulated quantities that depend on the global structure of phase space vs. those that depend on local trajectories.


I. INTRODUCTION
Quantum simulators promise solutions to a wide variety of problems in many-body system physics that cannot be solved efficiently by analytic approximation or classical numerical simulation [1][2][3][4][5]. Quantum simulation is widely seen as a potential application for noisy intermediate-scale quantum (NISQ) devices, comprised of ∼ 100 qubits, but unable to perform fault-tolerant quantum computation [6]. The absence of error correction in these NISQ-era quantum simulators raises an important question about their reliability. That is, can one trust the output of these devices, which are subjected to noise and imperfections [7][8][9]? Errors arise from a variety of causes, which are typically characterized as control errors due to miscalibration or inhomogeneities, uncontrolled classical noise, and decoherence due to entanglement with a quantum reservoir. Recently there have been various studies about the reliability and quantum advantage of NISQ devices in the presence of such errors for applications including simulation [10][11][12][13][14], optimizations [16], and random sampling [17].
Another source of error is "dynamical complexity." NISQ devices, even in the circuit model, are fundamentally analog devices, operating with a continuous set of quantum maps [8]. As such, dynamical instabilities such as bifurcations and the onset of chaos can translate into a quantum map that leads to a proliferation of errors. Indeed, quantum chaos can be characterized as the hypersensitivity of quantum dynamics to external perturbations in the Hamiltonian [18,19] and thus we expect that quantum chaos may limit our ability to reliably extract the desired output of analog (not digitally error corrected) quantum simulators, especially in cases where such devices are expected to solve hard problems.
The effect of quantum chaos on quantum computing was considered in the early 2000s [20][21][22][23]. Georgeot et al. studied how static imperfections such as qubit level fluctuations and residual (always on) interactions between qubits can lead to quantum chaos that destroys the register states of the system in the absence of error correction [20]. Song et al. studied the effects of errors in quantum algorithms while simulating the chaotic regime of the kicked-rotor model and showed that the error in the diffusive constant, which is associated with chaotic dynamics, grows exponentially with the system size [21]. More recently, dynamical complexity in quantum simulation has come to the fore. In a series of papers, Heyl et al. [24] and Sieberer et al. [25] studied the simulation of a Ising-type Hamiltonian through Trotterization in a gate-model. In this approximation, the unitary map consists of a series of Floquet maps describing the dynamics of a delta-kicked system, which is quantum-chaotic in a particular parameter regime. They have shown that the resulting magnetization errors increase sharply in this quantum-chaotic regime.
In this work we study the effect of quantum chaos on different aspects of quantum phase transitions (QPTs) in a quantum simulation. For models describable in the thermodynamic limit by mean-field theories [57], QPTs are often associated with bifurcations in the phase space dynamics that governs the order parameters. These bifurcations lead to unstable fixed points and separatrix lines, and we expect chaos to develop in their vicinity in the presence of small perturbations. Thus, it is natural to consider how quantum chaos affects the reliability of the quantum simulation in such models. We study this in the simplest paradigm, the Lipkin-Meshkov-Glick (LMG) model, which describes the anharmonic motion of a collective spin.
The LMG model undergoes a ground state quantum phase transition (GSQPT) in the mean magnetization [27], and a dynamical quantum phase transition (DQPT) in the time-averaged magnetization of the quench dynamics [28]. We will characterize the emergence of chaotic dynamics in the LMG model when a weak time-dependent perturbation is applied and analyze the resulting robustness of the simulation protocols. Importantly, we will show that a key quantity such as the time-averaged magnetization, extracted from this quantum simulation, can be highly sensitive to chaos in certain regimes. Nonetheless, we find that other aspects of the quantum simulation, such as the estimation of critical points for both the GSQPT and the DQPT, are robust to the presence of the perturbation and chaos. This difference in fragility vs. robustness can be explained by the difference between "fine-grained" and "coarse-grained" information being extracted from the quantum simulator.
The remainder of this manuscript is organized as follows. In Sec. II we review the LMG model, and analyze its thermodynamic limit to explain the underlying mechanism associated with both the GSQPT and the DQPT present in this model. Following this, in Sec. III we present a protocol that allows us to extract the critical points associated with both the GSQPT and the DQPT from a unitary quantum simulator. Then, in Sec. IV.A we characterize the chaotic dynamics emerging from the simulation of the LMG model due to the presence of a background time-dependent perturbation. Finally, in Sec. IV.B we analyze the robustness of two particular quantities, the time-averaged magnetization and the critical point estimates of the GSQPT and the DQPT, to the presence of a time-dependent perturbation.

II. QUANTUM PHASE TRANSITIONS IN THE LMG MODEL
The LMG model describes an infinite-range transverse field Ising model on a completely connected graph. We consider the specific case of the Curie Model [29] governed by the Hamiltonian are the collective angular momentum operators. While originally studied in the context of nuclear physics [30], the LMG model has applications in a wide variety of contexts, including the dynamics of twomode Bose Einstein condensates for applications of quantum metrology [31][32][33], studies of quantum chaos [34,35], and quantum simulation [11,36,37]. For the remainder of this manuscript, we scale the energy and parametrize the LMG Hamiltonian with a single parameter s according to with 0 ≤ s ≤ 1. The total angular momentum J 2 is conserved, allowing us to focus on the dynamics within the symmetric subspace with J = N/2, which is spanned by the 2J + 1 = N + 1 Dicke states.
The LMG model has a second-order continuous GSQPT at s = 1 2 in the thermodynamic limit, N → ∞, resulting from the competition between the external magnetic field inducing paramagnetic order and the spin-spin interactions inducing ferromagnetic order [27]. A signature of this phase transition can be seen on the energy spectrum, shown in Fig. 1a, from the closing of the energy gap between the ground state and the first excited state at around s = s (gpt) c (N ), which is different from the thermodynamic-limit value s (gpt) c = 1 2 due to finitesize effects. In this phase transition, the ground state changes its character continuously from paramagnetic to ferromagnetic in nature. The order parameter associated with this phase transition is the magnetization along the x-axis, ψ gs |J x |ψ gs . This order parameter is zero in the paramagnetic phase and nonzero in the ferromagnetic phase. The presence of this second order GSQPT can also be seen by noting the nonanalyticity in the secondorder derivative of the ground-state energy at the critical point as shown in the inset of Fig. 1a. In addition, the LMG model has a continuum of excited-state quantum phase transitions (ESQPTs) along the critical line, E ESQPT /J = −(1 − s), in thermodynamic limit [26,27]. A signature of these ESQPTs can be seen on the spectrum in the local divergence of the density of states along this line.
The LMG model also exhibits a DQPT in a order parameter associated with nonequilibrium dynamics [64]. In particular, the order parameter associated with this phase transition is the time-averaged magnetization along the x-axis, In the paradigm of DQPT, an initial state |ψ 0 , which is the ground state of the Hamiltonian H(s i ) for an initial value of the control parameter is time-evolved under a quenched Hamiltonian H(s f ), and the long-time dynamics of the state is analyzed. In this work, we fix the initial value of the control parameter to s i = 1, and choose the initial state to be one of the fully polarized states along the x-axis. At s = s  of the state changes from precessesing around the external magnetic field to evolving closely around its initial position under the action of the quenched Hamiltonian, resulting in a DQPT [11,28]. The order parameter, Eq.
(3), is zero for s f < s (dpt) c and has a nonzero value for

A. Semiclassical description of phase transitions
The LMG model is an example of a quantum meanfield model in which the thermodynamic limit, N → ∞, is equivalent to the dynamics of the mean field [57]. We obtain these dynamics from the expectation values of the Heisenberg equations of motion, neglecting all fluctuations by approximating AB ≈ A B . The resulting mean-field dynamics are thus equivalent to the classical dynamics of a massless "top." Defining the unit vector (X, Y, Z) ≡ lim J→∞ 1 J ( J x , J y , J z ), the classical equations of motion are [38] The above equations of motion can be solved analytically. As expected, d dt (X 2 + Y 2 + Z 2 ) = 0, indicating that the classical LMG Hamiltonian describes a system with a single degree of freedom, whose phase space is the unit sphere. Since the LMG Hamiltonian is timeindependent, the energy of the system is conserved, and thus the system is integrable and motion is regular.
The essential features of the classical (meanfield/thermodynamic) limit of the LMG Hamiltonian can be captured using the phase-space diagrams, such as the ones shown for different values of s in Fig. 1b. Note that for values of s smaller than 1 2 , all the phase-space trajectories precess around the z axis. At s = 1 2 a bifurcation occurs. The topology of the phase-space trajectories changes as the fixed point located at θ = 0, which is stable for s < 1 2 , becomes unstable for s > 1 2 and two new stable fixed points are formed, which are located at (θ, φ) = ± cos −1 ( 1−s s ), 0 (Here θ and φ are the usual angular coordinates on the sphere, θ = cos −1 (Z) and φ = tan −1 (Y /X)). This major reconfiguration of phasespace structure signifies the onset of the GSQPT. As a result of bifurcation of the fixed point, there are two kinds of trajectories present on the phase space for s > 1 2 , the ones that precess around the stable fixed points, and the ones that revolve around the whole sphere. These two kinds of trajectories are separated by a separatrix layer, which includes the unstable fixed point located at θ = 0. As the value of s is increased from 1 2 to 1, these stable fixed points move farther away from the unstable fixed point and the separatrix layer bounds more trajectories on the phase space.
The DQPT also admits a clear description in terms of the classical trajectories. The system is first initialized in the fully polarized state along the x-axis, which is a ground state at s = 1, and then allowed to evolve under the action of the quenched Hamiltonian, H(s f ). In the thermodynamic limit, as the value of s f is varied, the position of the initial state changes with respect to the separatrix layer, resulting in different dynamical phases Here, the initial conditions are spin coherent states centered roughly at the points of the classical protocols (see text for details). [28]. For s f < 2 3 , the state precesses around the z-axis resulting in a zero value for the order parameter (see Fig.  2a and Fig. 2b). At the dynamical critical point, s f = 2 3 , the state is exactly located on the separatrix layer, while for s f > 2 3 , the state precesses around one of the stable fixed points, as shown in Fig. 2c, resulting in a nonzero value of the order parameter for these s values.
We gain further insight into the dynamics in the semiclassical limit by mapping the LMG to a problem of a fictitious particle moving in a one-dimensional well in the limit N 1 [27,31,[39][40][41][42][43]. Following [31], one can derive a time-dependent Schrodinger equation in the Fock basis, obtained by expressing the spin operators in the Schwinger representation. Ignoring terms of order O(N −3 ) and smaller, the time-dependent Schrodinger equation along the x-axis can be written as This can be interpreted as the Schrodinger equation for a particle with a position-dependent mass, moving in a plus Ndependent corrections. Importantly, as shown in Fig. 1b, this potential well V s (x) changes its shape from a single well to a double well at s = 1 2 as s is increased from 0 to 1 in this large-N limit.
Classical motion of a particle in this potential also provides intuition to explain the presence of the GSQPT and the DQPT in the LMG. At s = 1/2 the potential changes its shape from a single well to a double well, corresponding to the critical point of the GSQPT in the thermodynamic limit, which is accompanied by the spontaneous symmetry breaking in accordance with Landau-Ginzberg theory [44]. Note that for a finite value of N , the critical point, s (gpt) c (N ), moves to a higher value of s compared to the thermodynamic limit value, s (gpt) c , due to the finite-size zero-point energy of the ground state, which is higher than the barrier height of the shallow well for finite-N values. Therefore, one can identify the value of s at which the zero-point energy equals the barrier height of the double well as the finite-size critical point. For more details refer to Appendix A, where we derive the finite-size scaling of the ground-state critical point.
The DQPT can be understood as a transition across the separatrix line. For the specially chosen initial condition, at s = 2/3 the system in the classical description is exactly on the separatrix in phase space; its energy is equal to the height of the barrier of the double well. For s < 2/3 the pseudoparticle is above the barrier and executes librations in the double well. This corresponds to the top precessing around the external magnetic field in the LMG model, resulting in the order parameter X = 0 in the long-time limit. In the second dynamical phase, for s > 2 3 , the energy is lower than the barrier height of the double-well potential. As a result, the pseudoparticle becomes trapped on one side of the double well, corresponding to precession in phase space around one of the stable fixed points, and X = 0. This explains why the critical value of s is greater for the DQPT than the GSQPT. The double-well forms first at the bifurcation point, yielding the GSQPT; for larger values of s the barrier height becomes equal to the energy of the specially chosen initial condition yielding the DQPT.
Moreover, the double-well potential in the classical picture also explains the presence of the ESQPTs [27]. A classical particle with the energy equal to the separatrix energy spends a very long time in the vicinity of the unstable fixed point, which implies that the probability of finding the associated quantum state is very high in such neighborhood. Hence, the state undergoing an ESQPT is localized around the unstable fixed point. Since the barrier height of the double well increases with the value of the control parameter, the energy of the state undergoing ESQPT increases as well, and the energy of the phase space separatrix is E ESQPT /J = −(1 − s), in the classical limit.

III. ACCESSING QPTS WITH DYNAMICAL QUANTUM SIMULATORS
We are interested in studying QPTs in a quantum simulator whose ideal operation involves implementing a desired unitary transformation on a chosen initial state, followed by measurement and data processing to extract a desired output [11,36]. On this type of device, some states, such as spin-coherent states, can be prepared with high fidelity, and the expectation values of certain observables are accessible from measurements after the system evolves under the Hamiltonian of interest for a chosen duration. With this setting in mind, we seek to study the critical phenomena associated with both the GSQPT and the DQPT in the LMG model, and in particular we want to extract the critical points associated with these phase transitions. Recently, protocols have been proposed to identify the critical points in some of the Ising-like models by measuring out-of-time-ordered (OTO) correlators [45], or studying the spectrum of multiple quantum coherences [46], or through adiabatic quenches [47]. In this section, we propose a protocol for identifying the critical point of the GSQPT, which only requires us to have the ability to prepare a particular spin-coherent state and the ability to measure the collective magnetization in the longitudinal direction, ψ(t)|J x |ψ(t) , as a function of time. We first motivate our protocol based on the classical phase-space trajectories and then provide justification on why this protocol would still work in the case of a finite-dimensional system. Moreover, we will see that changing the initial state in the above mentioned protocol also allows us to measure the order parameter associated with the DQPT and thereby allowing us to extract the DQPT critical point.

A. Classical bifurcation
As discussed in Sec. II A, in the thermodynamic limit a stable fixed point at θ = 0 bifurcates into two other stable fixed points resulting in the change of topology of the phase-space trajectories at the critical point of the GSQPT, in accordance with the Landau-Ginzberg theory. Therefore, as shown in Fig. 1b, the single-well potential changes into a double-well potential at s = 1 2 as s is increased from 0 to 1. Based on this semiclassical picture, we propose a protocol to identify the critical point of the GSQPT by computing the time-averaged magnetization, X, of an initial condition located at an angle slightly off the z-axis, (θ 0 = T L , φ 0 = 0), where T L is a small angle. We show in the Appendix B that the critical point estimate of the GSQPT is robust to perturbations in this angle. In this work, we choose T L = π 60 . This time-averaged magnetization, when plotted as a function of s, exhibits a pitchfork bifurcation at the GSQPT critical point [48], as shown in Fig. 2b. In this figure, the zero values of the time-averaged magnetization correspond to the values of s where the initial condition is oscillating in the single-well potential centered at θ = 0 as shown in Fig. 1b. The nonzero values correspond to the values of s where the initial condition's dynamics is constrained to one side of the double well centered at (θ = cos −1 ( 1−s s ), φ = 0 . Likewise, the time-averaged magnetization for the initial condition located at (θ 0 , φ 0 ) = ( π 2 , 0), which is the order parameter associated with the DQPT, also shows a pitchfork bifurcation at s = 2 3 , as shown in Fig. 2c, because the initial state can access the whole well, singlewell or double-well potential, for s < 2 3 while its dynamics is constrained to one side of the double well for s > 2 3 [28]. Note that this bifurcation is expected because X is the order parameter for the DQPT, which is zero in one phase and nonzero in the other phase by definition. As a result, this bifurcation protocol can be used to determine the critical points of both the GSQPT and the DQPT in the thermodynamic limit provided appropriate initial conditions are used, (θ 0 = T L , φ 0 = 0) for GSQPT and (θ 0 = π 2 , φ 0 = 0) for DQPT. Note that, as a consequence of the parity symmetry, one could also construct the negative branch associated with both the bifurcation diagrams in Figs. 2b and c by starting the protocol in initial states that are rotated by π about the z-axis from the initial states corresponding to the positive branch.
Since LMG model is a system with one degree of freedom, energy conservation allows us to gain complete knowledge of the system. In particular, we can derive an analytic expression for the time-averaged magnetization for all the initial conditions of the form (θ 0 , φ 0 = 0), which is given by where Λ(θ 0 , s) = − 4(1−s) s sin 2 (θ0) cos(θ 0 ) − 1−s s , and K(x) denotes the complete elliptic integral of the first kind. We present the full derivation and further details in Appendix B. The time-averaged magnetization in the LMG model has been explored in various previous works [11,49,50] and a similar calculation is presented in Ref. [51]. Note that the above result also provides us with the bifurcation points along with their respective bifurcation curves for all initial conditions of the form (θ 0 , φ 0 = 0). The two cases of interest for us are the θ 0 = T L (GSQPT) and θ 0 = π 2 (DQPT), whose bifurcation points are given by (1 + cos 2 ( T L 2 )) −1 0.5 and (1 + cos 2 ( π 4 )) −1 = 2 3 respectively. Furthermore, the above expression for the bifurcation point also explains its robustness with respect to the choice of initial θ for the GSQPT case, which can be inferred from the absence of the first order term in the Taylor expansion of (1 + cos 2 ( θ 2 )) −1 around θ = 0.

B. Quantum bifurcation
In this subsection, we show that even for a finitedimensional system, the protocol described above allows us to estimate the finite-size critical points of both the GSQPT and the DQPT. For the GSQPT, a spincoherent state centered at an angle slightly off the classical unstable fixed point, (θ = (N ), φ 0 = 0), labelled by |θ = (N ), φ 0 = 0 is time-evolved under the LMG Hamiltonian for different values of s. Here, we set (N ) = π 60 + 1 √ N , which accommodates the variance of the spin-coherent state and approaches T L in the thermodynamic limit. The resulting time-average from the timeevolution, X = J x /J, is then analyzed as a function of s, as shown in Fig. 2d. The finite-size critical point can be estimated as the value of s in which X changes from zero to a nonzero value. An important point here is the choice of the averaging time, denoted by T , which in this case should be larger compared to the time-scales associated with the periods of oscillation localized in a single-well, but also smaller compared to the time scales associated with the tunnelling time between the two sides of the double well. This can be inferred from the following argument.
Consider the expectation value of J x averaged over time, J x , which can be written as where the initial state |ψ(0) in Eq. (7) has been expressed in the eigenbasis of H, |ψ(0) = n c n |u n resulting in Eq. (8). Note that the above expression becomes J x = n |c n | 2 u n |J x |u n for the averaging time T >> 2π En−Em because the time-averaging integral in Eq. (8) becomes a Kronecker delta function, δ nm , assuming that the energy states are non-degenerate. Since J x is odd under the action of the parity operator, e iπJz , the above expression finally becomes zero for T >> 2π En−Em . In the LMG Hamiltonian, the gap between opposite parity energy eigenstates decreases as s is increased from 0 to 1 and starts to close exponentially with the system size in the ferromagnetic phase [52], s > s (gpt) c (N ), signalling the onset of the GSQPT. Therefore, we find that choosing the average time T such that its much larger than the typical inverse energy gap for s < s c , and using the same value of averaging time, T , for all s values does a good job at reproducing the expected bifurcation.
For the case of the DQPT, the ground state at s = 1, which is a spin-coherent state centered at (θ = π 2 , φ 0 = 0), is time-evolved under the action of the LMG Hamiltonian for different values of s to obtain the the bifurcation diagram shown in Fig 2e. Since the true order parameter of the DQPT is being measured in this case, we expect this protocol to give us an estimate of the finite-size critical point for all cases.

A. Chaos in dynamical quantum simulation
The LMG Hamiltonian describes a system with one degree of freedom with no explicit time dependence implying that the classical system is integrable. This integrability can be broken upon addition of a perturbation that does not preserve the original symmetries of the system. For instance, in Ref. [53] the authors explore the impact of adding nearest-neighbor interactions in the LMG model. This perturbation, which breaks the permutational invariance of the system, leads to a new dynamical phase, where slight changes in the parameter of the Hamiltonian results in a very different dynamical behavior of the system. In this work we study the effect of a weak time-dependent perturbation that renders the system nonintegrable, with the goal of studying the impact of chaos on the quantum simulation of critical quantities associated with both the GSQPT and the DQPT. Here we choose a perturbation of a weakly-oscillating B-field along the y-axis, which results in the following effective Hamiltonian Mourik et al. showed that a Hamiltonian of this form shows chaotic behavior for appropriate values of ε 0 and s [35]. To evaluate the emergence of chaos in the perturbed system, we computed the fraction of phase space that becomes chaotic in the presence of this J y perturbation term for a range of frequencies and s values, which is shown in Fig. 3. Each data point on the heat map was obtained by identifying the fraction of initial conditions whose distance from their corresponding neighboring points separated exponentially for a given value of frequency and s using the appropriate classical equations of motion. This data is in excellent agreement with the associated Poincaré sections, similar to that seen in Ref. [35]. Note that this system becomes chaotic mostly only for s > 1/2 as shown in Fig. 3. This can be explained by the presence of the separatrix layer on the phase space only for s > 1/2. As it is well known, when an integrable system is made chaotic by adding a perturbation, chaos in the system first originates in the vicinity of this separatrix layer through the so-called homoclinic tangle [54,55]. As a result, there is a very thin layer of chaos around this separatrix even for very small perturbations, reminiscent of the conservative Duffing oscillator [54], which is a tilted-oscillating double well [65].
It should also be noted that the existence of phase transitions is closely related to the presence of unstable fixed points on the classical phase space [28,41]. A crucial question then stands out: given that these points are also the regions of the classical phase space where chaos first develops when a nonintegrable perturbation is added to the system, how does chaos affect our ability to simulate quantum phase transitions in a noisy quantum device? We address this question in the following subsection by analyzing the effects of a nonintegrable perturbation on the bifurcation diagrams.

B. Sensitivity and robustness to perturbations in the simulation of QPTs
The bifurcation diagrams associated with the GSQPT and the DQPT, whose bifurcation points provide us the critical-point estimates, are shown in Figs. 4a and b in the presence of the aforementioned chaotic perturbation for J = 100 (N = 200). The black curves in these figures represent the ideal bifurcation diagrams, whereas the red and the green curve represent the perturbed cases associated with the GSQPT and the DQPT correspondingly. Note that, for the case of the GSQPT, the perturbation has significantly modified the bifurcation curve for all values of s s (gpt) c (N ), whereas the perturbed-DQPT bifurcation diagram has been altered significantly only in a small range of s values. This difference in behavior of the perturbed bifurcation diagrams can be attributed to the way chaos emerges in the system when a nonintegrable perturbation is added. Note that the initial state associated with the GSQPT protocol is always in the immediate vicinity of the unstable fixed point and therefore very close to the separatrix layer for all s > s (gpt) c (N ). As a result, the initial state in this case is always in the chaotic sea for s > s . In contrast, the initial condition associated with the DQPT is closer to the separatrix layer only around the DQPT critical point, and therefore present in the chaotic sea only for these s values as shown by the green dot in the above-mentioned figures.
The different behavior of the GSQPT and the DQPT bifurcations upon addition of a perturbation can be further corroborated by analyzing the Lyapunov exponents for the associated initial conditions as a function of s as shown in Fig. 4c, where the red curve and the green curve label the Lyapunov exponents associated with the . The Lyapunov exponents here are computed using the method introduced in [56]. Also, note that the depth of the effective potential well increases as s is increased, which means that the perturbation does not affect the potential well as much for larger values of s compared to smaller s values. This can also be observed in the Poincaré sections shown in Fig. 4f, where the potential wells are still intact surrounded by a chaotic sea. This explains the robustness of the DQPT bifurcation diagram for larger values of s.
So far, we've focused our analysis on a specific perturbation amplitude ε 0 = 0.05, but the behavior observed here is fairly general. That is, for any perturbation amplitude, there is a significant error in the time-averaged magnetization for all s > s The results shown in Figs. 5a,b demonstrate that the onset of chaos makes the system more sensitive to pertur-bations around the critical regions. However, a crucial question is, what is the particular quantity one wants to extract from a quantum simulation? For some applications, one is not interested in the exact value of the time-averaged magnetization, but in extracting the value of the critical point, i.e, the value of the control parameter s which separates the two phases involved in the QPT. In our protocol, this can be done by identifying, for each bifurcation diagram, the value of s at the time-averaged magnetization becomes larger than a small threshold X th << 1. The resulting critical-point estimates associated with both the quantum phase transitions are then shown in Fig. 5c as a function of the perturbation strength, ε 0 . Importantly, we observe that these critical point estimates are not significantly affected by the perturbation, and remain close to the ideal values even when the perturbed system is highly chaotic and the associated time-averaged magnetization have been modified significantly with respect to the ideal case. This can be understood to be a consequence of the fact that the perturbation renders the dynamics of the system chaotic only for s s c (for both types of QPTs). As a result, the time-averaged magnetization undergoes a drastic change, with or without the presence of the perturbation, between being zero before the transition (with dynamics largely unaffected by perturbation) to being nonzero after the transition, either due to the change in phase for ε o = 0 or the system turning chaotic for ε o > 0. In this way, the robustness of the critical point estimates can be attributed to the fact that they signal a global change in the structure of the LMG phase space, which occurs with or without the perturbation.
Although the perturbation changes the value of the critical point estimates slightly for both the phase transitions, it does significantly affect the sensitivity of critical point estimates to the value of the threshold X th = sin( δ √ 2J ) used in identifying these points. The critical point estimates as a function of threshold value are plotted in Figs. 6a,b, where the ideal-critical-point estimates have been labelled by the black curve and the perturbedcritical-point estimates by the blue and the magenta colored curves for both the GSQPT and the DQPT. Notice that the perturbed-critical-point estimates vary over a wide range of s values, particularly for the GSQPT, whereas the ideal estimates are reasonably bounded. This sensitivity could play an important role in experiments where the lack of resolution in the time-averaged magnetization might make the critical point estimates more sensitive to the perturbation.

V. SUMMARY AND DISCUSSION
In this work, we have studied the effects of chaos arising from a perturbation on the quantum simulation of QPTs in the LMG Hamiltonian. One expects chaos to play an important role in the reliability of such simulations given the connection between QPTs, bifurcations, unstable fixed points, and instabilities that arise in the semiclassical behavior of an order parameter in the thermodynamic limit. We have proposed a protocol that allows us to extract the critical points of both the GSQPT and the DQPT by measuring the time-averaged expectation values of the collective spin operators. This is motivated by the dynamics in the semiclassical limit, which can be described by the motion of a fictitious particle in a one dimensional potential well. Both the GSQPT and the DQPT are intuitively understood in terms of the motion in the potential well. Given the integrability in the unperturbed case, we are able to derive an analytic expression for the time-averaged magnetization, which describes the output of the ideal simulation in the thermodynamic limit. This protocol also allows us to estimate the finite-size critical point in a finite-dimensional system provided the appropriate time-averaging is performed.
To study the effect of chaos, we added a weak external time-dependent perturbation during the simulation of the LMG Hamiltonian, which becomes chaotic in a manner reminiscent of the conservative Duffing oscillator. The time-averaged magnetization associated with both the GSQPT and the DQPT are affected significantly, particularly at the values of the control parameter s closer to the critical point. Furthermore, we have shown that the range of s values where time-averaged magnetization have been affected the most can be identified based on the location of the initial state of the protocol with respect to the separatrix layer on the classical phase space, as the homoclinc tangle forms in the transition to chaos. Finally, despite of the sensitivity of the time-averaged mag-netization to the presence of the chaotic perturbation, we showed that the critical point estimates, obtained from the analysis of the perturbed evolution, are robust for both the GSQPT and the DQPT. This robustness can be attributed to the fact that the critical points signal a change in the global structure of the LMG phase space that is also captured by the perturbed evolution.
We have studied here the effect of chaos on quantum simulation of GSQPT and DQPT for the simplest case, the LMG model, a completely connected Ising model with two-body interactions. In future, it would be interesting to explore the effect of chaos on the signatures of ESQPT. Given that the states undergoing ESQPTs are localized around the unstable fixed point, we would expect the associated signatures of these phase transitions to be fragile to the presence of the chaos, particularly in the vicinity of their critical points. In addition, we also plan to analyze more complex many-body systems. A first rich generalization is to analyze p-spin models, with all-to-all p-body interactions [57][58][59]. In contrast to the case p = 2 (the LMG model), the cases of p > 2 exhibit discontinuous first-order QPTs and are characterized by a gap that closes exponentially with N [58], making quantum simulation more challenging in some protocols. Moreover, it has been shown that the phasespace associated with the kicked-version of these models undergo rapid changes at various bifurcation points compared to the LMG model [59]. We expect these new classes of dynamical instabilities to impact the reliability of the quantum simulation of p-spin models. Beyond quantum mean-field models, describable by a collective spin with one degree of freedom, it is important to study the impact of many-body chaos on more general manybody models.
In the context of NISQ quantum simulators that do not have access to error correction, it is important to identify quantities that are robust with respect to the error but also hard to simulate classically [6][7][8]. In this context, what is the relationship between the hardness of quantum simulation and its robustness to perturbations [7,8]? In particular, are computationally hard analog quantum simulations limited by many-body chaos, while computationally simulable problems are robust to dynamical perturbations? The work presented here gives some initial indications in this direction. Identifying the quantum critical points is robust in this case because these correspond to changes in global structure of phase space, which was possible to identify even without keeping track of the exact quantum state. On the other hand, other quantities that depend on simulating the exact trajectory are fragile. In a general many-body system, the question is then, which structural changes in the effective phase space can be efficiently simulated, and which emergent critical phenomena cannot, and how robust are they in the presence of dynamical instabilities?

VI. ACKNOWLEDGEMENTS
We thank Poul Jessen and Manuel Muñoz for numerous insightful discussions. This work was supported by the U.S. National Science Foundation under grant numbers PHY-1820679 and PHY-1630114.
Appendix A: Finite-size scaling In this subsection, we derive an analytic expression for the scaling of the critical point with the system size for the GSQPT, s (gpt) c (N ). We obtain this expression using the effective Schrodinger equation introduced in section II A, which is a good approximation in the semi-classical limit. The main idea in this derivation is to identify the critical point of the GSQPT with the value of s at which the zero-point energy of the double well equals the barrier height of the double-well. The potential of the semi-classical well is given by The above expression for the potential allows us to compute the barrier height as shown below In addition, the zero-point energy can also be computed using the potential energy expression in Eq. (A1). The main idea is to notice that the potential energy is of the form V (z) = V 0 + 1 2 mω 2 z 2 , and therefore the zero-point energy is given by 1 . As a result, the zero-point energy is given by Note that the above expression agrees with the predictions of the finite-size scaling ansatz [60], sc(N )−sc sc ∝ N − 1 ν * where s c (N ) is the effective critical point for finitesize system N and ν * is the exponent associated with the divergence of the coherence number at the critical point [63]. In the case of fully connected models, this exponent ν * is related to the mean-field exponent ν M F through the upper critical dimensionality of the corresponding finiteranged model, ν * = d c v M F . For the case of LMG, d c = 3 and v M F = 1 2 , therefore v * = 3 2 [63]. Note that this result v * = 3 2 is also consistent with the results in [40,[61][62][63].

Appendix B: Classical bifurcation
In this section, we first derive the critical point of the bifurcation diagrams and then derive an analytic expression for the time-averaged magnetization. Given an initial condition of the form (θ = θ 0 , φ 0 = 0), energy conservation leads to the following relationship between θ(t) and φ(t): which is the value of s when the second solutions starts to exist. It is interesting to note that the bifurcation point corresponding to GSQPT is more robust to deviations in θ ini than the one for DQPT as shown in Fig. 7.
For a given initial condition, the analytic expression for the time-averaged magnetization, X, can be obtained by averaging X(t) over one time period of the classical trajectory as shown below Using the above equations of motion, the integrand can be expressed as dt = dθ s 2 sin(θ) sin(2φ) . Also conservation of energy allows all the φ(t) terms in the integral to be expressed as a function of θ(t). Evaluating the integral for an initial condition of the form (θ = θ 0 , φ = 0) we obtain, The time period, T , can then be evaluated as follows because the limits of the integral were chosen assuming this condition holds true. Note that ∓ has been added to the above expression to account for the fact that the direction of trajectory switches as the fixed point passes through the initial condition, which happens at s = 1 1+cos θ0 . Using various identities, the above expression can be reexpressed as where K denotes the complete elliptic integral of the first kind and Λ(θ 0 , s) = − 4(1−s) s sin 2 (θ0) cos(θ 0 ) − 1−s s