Fragility to quantum fluctuations of classical Hamiltonian period doubling

We add quantum fluctuations to a classical period-doubling Hamiltonian time crystal, replacing the $N$ classical interacting angular momenta with quantum spins of size $l$. The full permutation symmetry of the Hamiltonian allows a mapping to a bosonic model and the application of exact diagonalization for quite large system size. In the thermodynamic limit $N\to\infty$ the model is described by a system of Gross-Pitaevskii equations whose classical-chaos properties closely mirror the finite-$N$ quantum chaos. For $N\to\infty$, and $l$ finite, Rabi oscillations mark the absence of persistent period doubling, which is recovered for $l\to\infty$ with Rabi-oscillation frequency tending exponentially to 0. For the chosen initial conditions, we can represent this model in terms of Pauli matrices and apply the discrete truncated Wigner approximation. For finite $l$ this approximation reproduces no Rabi oscillations but correctly predicts the absence of period doubling. Our results show the instability of time-translation symmetry breaking in this classical system even to the smallest quantum fluctuations, because of tunneling effects.


I. INTRODUCTION
The experimental discovery [1,2] of Floquet timecrystals few years after their theoretical prediction has been a real breakthrough. In analogy to ordinary crystals, time crystals appear as a consequence of breaking time-translation symmetry in the system [3][4][5]. Time crystals were first introduced in 2012 by Frank Wilczek [6]. Following earlier attempts to identify systems able to display time-translation symmetry breaking, in 2015, a no go theorem by Watanabe and Oshikawa showed that this is not possible in the ground state or in thermal equilibrium [7]. Among many possible non-equilibrium candidates, periodically periodically-driven (Floquet) systems have proven to be the most promising realization. Stimulated by the initial proposals [8,9], a large body of theoretical work has been performed [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. A common ingredient to all case is the presence in the dynamics of a sufficient number of constraints that introduce ergodicitybreaking, thus impeding the system to reach an effective infinite temperature.
Nearly all the attention, so far, has been devoted to quantum systems. Only few notable exceptions [25][26][27][28] consider classical dynamics. Especially interesting is the case of driven classical many-body Hamiltonian systems, where a long-lasting prethermal regime has been found [29,30], and period doubling (or period n-tupling with n > 2) can appear in the prethermal regime [31,32]. All these systems eventually thermalize after a transient, and this fact relies on their chaotic dynamics. Chaos is the generic situation for a finite number of coupled classical Hamiltonian systems [33], but the situation can drastically change in the thermodynamic limit for longrange interacting systems [33][34][35]. The phenomenon of sub-harmonic generation (period-doubling) in a classical Hamiltonian driven many-body system was recently considered in Ref. [36] and it was termed Hamiltonian synchronization (or classical Floquet time-crystals). One question is if this synchronization phenomenon is stable to fluctuations. In [36] this stability was discussed against thermal fluctuations, here we explore the stability against quantum fluctuations.
Besides addressing the problem of stability to fluctuations, the present work aims to make a first step towards a model that has a time-crystalline phase both in the classical and in the quantum regime, so to understand their difference. Most simply, we substitute the classical angular momenta with quantum spins of magnitude l and find that, whenever l is finite, the quantum fluctuations destroy the synchronized period-doubling motion. It is recovered only in the limit of infinite spin magnitude l → ∞, when the dynamics becomes classical again.
In all the paper we focus on the case where the interactions are all-to-all and the correlations are therefore very strong. Moreover, the kicking exactly flips the spins and we take the initial state as fully polarized up. If the quantum fluctuations destroy the period doubling in this most favorable situation, they will destroy it also in case of imperfect flipping and faster decaying fluctuations. What we find here is that adding even the smallest quantum fluctuations (l 1 finite), one spoils the time-translation symmetry breaking in this model. Due to quantum tunneling, some Rabi oscillations incommensurate with the driving period add on the period-doubling oscillations. The response is no more synchronous with the driving and there is no more a persistent period-doubling response, so there is no more time crystal.
The paper is organized as follows. In Sec. II we introduce a "period-doubling order parameter", a quantity which first vanishes at a time increasing with the system size if the system shows persistent period doubling in the thermodynamic limit. We add quantum fluctuations to the classical backbone of [36]. We do this in two ways, and we get two quantum models, both reducing to the classical one when l, a parameter we are going to ii describe, tends to infinity (we discuss this limit in some detail in Appendix A).
In the model-1 we simply substitute classical angular momenta with quantum spins of finite size l and discuss it in Sec. IV. By using a mapping to a bosonic model [18,37] (Appendix B) and exact diagonalization for finite system size, we see that the period-doubling order parameter first vanishes at a time not scaling with the system size, which marks the destruction of period doubling. By studying the average level spacing ratio in Sec. IV A we find that the dynamics leading to this result is related to quantum chaos.
In Sec. IV B we perform the thermodynamic limit and show that the system is here described by a system of Gross-Pitaevskii equations. In that limit we see the period-doubling order parameter performing Rabi oscillations, so there is no period doubling. We see that the period of these oscillations diverges with the spin magnitude l and in the limit of infinite spin the period doubling is recovered. This is in agreement with the fact that the quantum fluctuations disappear in this limit.
In agreement with the finite-size quantum dynamics, the classical infinite-size Gross-Pitaevskii dynamics is chaotic, as the largest Lyapunov exponent shows (Sec. IV B 1), but it is not fully ergodic and the Rabi oscillations can persist. Studying the amplitude and the frequency of the Rabi oscillations versus the parameter K for different values of l, we see that the curves show a crossing point for K ∼ 1, which corresponds to a transition from synchronized to trivial behaviour in the classical l → ∞ limit.
Rabi oscillations are related to the ones obtained in [17] for a single spin system. Coupling many of these systems with a small coupling K the oscillations are still there but with a renormalized period; a large coupling on the opposite leads to the destruction of the Rabi oscillations and to small chaotic oscillations of the period-doubling order parameter. The correlations induced by the coupling are never strong enough to stabilize the period-doubling order parameter to a persistent finite value, against the quantum fluctuations.
In Sec. V we study the model-2, where each classical angular momentum is substituted by an average of 2l Pauli matrices. We study this case by means of the discrete truncated Wigner approximation (DTWA), which we summarize in Sec. V A and is known to give good results for long-range interactions [38][39][40]. Also here we find find the disappearance of the period doubling (Sec. V B): the period-doubling order parameter decays as an exponential in time and the decay time scale does not scale with N . We see that the decay time increases for increasing value of l as a power law. So, for l → ∞, where the system behaves classically, the period doubling persists for an infinite time, as expected. In Appendix C we discuss a different way to estimate this decay time which gives consistent results and discuss some technical aspect related to DTWA.
We remark that, for the chosen initial state, the model-2 is equivalent to the first one but the DTWA gives results in quantitative agreement only in the limit l → ∞. For l finite it is only correct in predicting the absence of period doubling in the limit of large N but provides no Rabi oscillations.

II. THE MODELS
We introduce quantum fluctuations in the model studied in [36]. It is a chain of N coupled classical angular momenta undergoing a periodic pulsed driving. Here we will focus on the case with all-to-all interactions. These ones give rise to the strongest long-range correlations needed in order to stabilize a possible period-doubling phase. Indeed, in the classical case this model shows a phase with persistent period doubling in the thermodynamic limit, also in the all-to-all case. Adding the quantum fluctuations, we will show that the period doubling in the all-to-all interacting case disappears. This result implies the absence of period doubling also for faster decaying interactions (and smaller long-range correlations). The Hamiltonian is where δ τ ≡ n δ(t − nτ ) [41] and we put a factor N in the denominator in order to ensure extensivity. The m α j , α = x, y, z are the components of classical angular momenta which obey the angular-momentum Poisson brackets m µ i , m ν j = µ ν ρ δ i j m ρ j where µ ν ρ is the Ricci fully antisymmetric tensor. For K and h small enough and φ in a neighborhood of π/2 this classical Hamiltonian model shows a persistent period-doubling behaviour [36].
This classical model is such that when K = 0 it is equivalent to a single degree of freedom showing entrainment with the driving, that's to say it shows a response synchronized with the one of the driving, with a period doubled respect to the driving [17]. When K = 0 and N is finite, this response dies after a transient. For a region in the parameter space, the duration of this transient diverges with the system size going to infinity [36]. So, for N → ∞, the system shows persistent collective oscillations with a period double with respect to the driving, in which all the spins behave in a synchronous way. This is a form of period-doubling time crystal, as we discuss in Sec. III.
In order to add quantum fluctuations to this model we can quantize the angular-momentum variables replacing iii them with quantum spins. The resulting Hamiltonian iŝ x iŝ whereŝ α j , α = x, y, z are quantum spins of magnitude l (ŝ 2 j = l(l + 1)) obeying the commutation rules [ŝ µ ,ŝ ν ] = i µνρŝρ . We callĤ (1) (t) as the model-1. Another possibility, which should give results physically similar to the first one, is performing the following substitution So we replace the classical angular momenta in Eq. (1) with an average of 2l Pauli matrices, and then we multiply by 2l. We call the resultinĝ as the model-2. The parameter l has the same symbol here and in the model-1 on purpose. Indeed, alsom α j are spin variables and because we choose as initial state the one fully polarized up (see Eq. (16)), these are spins of size l, as well known from the rules of addition of angular momenta [42]. So, with our initialization, the variableŝ s α j of the model-1 and the variablesm α j of the model-2 are exactly equivalent. In some sense, the model-2 is a spin-1/2 representation of the first one, amenable to be described by means of DTWA.
For any finite l there are quantum fluctuations around the classical backbone Eq. (1). When l → ∞ the fluctuations become irrelevant and both the models tend to become classical. This can be seen, for instance, by using exactly the same methods discussed for the Lipkin-Meshkov-Glick model in [43,44]). For completeness, we give a sketch of this analysis in Appendix A.
In the rest of the paper we numerically study the two models. We study the model-1 in Sec. IV and the model-2 in Sec. V. In both cases we will consider the stroboscopic dynamics, that's to say we will focus on times which are an integer number of periods t = nτ . More precisely, we will chose the time nτ as the time immediately before the n-th kick. We will show that whenever there are quantum fluctuations -that's to say for any finite lthere is no period-doubling phase and in the limit l → ∞ one recovers the period doubling, consistently with the attaining of the classical limit.

III. PERIOD DOUBLING AND TIME-CRYSTAL BEHAVIOUR
In order to make our paper self contained, we briefly recap the main ideas about time-crystal behaviour, which appears as a period doubling in the classical limit of our model. Time-crystal behaviour is a synonym for timetranslation symmetry breaking: a driven system in the thermodynamic limit shows a response with a frequency multiple with the driving one. Thereby the discrete time translation symmetry of the driving is broken. In order to spot time-translation symmetry breaking -or its absence -it is very important to define precise criteria which are able to distinguish this complex collective phenomenon from analogous single particle effects. Summarizing the discussion of Refs. [8][9][10] -where the relevant criteria and conditions to have a Floquet time crystal were introduced -we can state that there must exist an observable O and a class of initial states |ψ such that, considering stroboscopic times t = nτ , the expectation value in the thermodynamic limit (N → ∞) satisfies all of the three conditions II) Rigidity: f (t) shows a fixed oscillation period τ B (for instance τ B = 2τ , the so-called "period doubling") without fine-tuned Hamiltonian parameters.

III)
Persistence: the non-trivial oscillation with fixed period τ B must persist for infinitely long time, when the thermodynamic limit N → ∞ in Eq. (5) has been performed.
We will focus here on period doubling, τ B = 2τ . In summary we seek for a quantity -called "order parameter" in analogy with standard symmetry breaking -such that it oscillates with frequency 2τ for an infinite time in the thermodynamic limit (when the size of the system N tends to infinity). In our model (model-1 and model-2 are essentially equivalent) there are some limits where such a quantity can be found. For instance, in the limit K → 0, our model reduces to the kicked Lipkin-Meshkov-Glick model of [17] and the order parameter is provided by Here the role of the system size is played by l which measures the number of interactingσ j, m spins which compose theŝ z j in the model-2 representation.
Another interesting limit is the l → ∞ limit (with K = 0). In this limit the model is classical (see Appendix A) and can show persisting period doubling in the thermodynamic limit (in this case N → ∞) [36]. In this case the order parameter is iv Taking l finite, it is quite natural that, if there were period doubling, it would appear in the finite-l version of s(t), namely In order to see if this quantity shows persisting oscillations with period 2τ (period doubling), we focus on its finite-N version and perform a finite-size scaling in N . We focus therefore on whereŜ z = N j=1ŝ z j . We put the multiplying factor (−1) t/τ , because a period-doubling is expected to imply a change of sign ψ(t)|Ŝ z |ψ(t) at every period [17,36]). Thanks to the multiplying factor, the period doubling would appear as a never-vanishing value of O(t), which is easier to study.
In order to probe if there is persistent period doubling in the thermodynamic limit, one should check the presence of the following finite-size scaling: If O(t) first vanishes after a time t * scaling with N towards infinity, then one has period doubling [36]. So, in the thermodynamic limit O(t) never vanishes and there is persistent period doubling. In the rest of the paper, we call for conciseness O(t) the "period-doubling order parameter", even if in the light of the discussion above this is a slight abuse of terminology.
For l finite, we will see that t * never scales with the system size, implying the absence of persistent period doubling and time-crystal behavior.

IV. ANALYSIS OF MODEL-1
The Hamiltonian is given in Eq. (2). In order to probe the existence of a possible persistent period doubling, we initialize the system in the state where all the spins are in an eigenstate of the correspondingŝ z j with eigenvalue l. This is the most favorable condition for the appearance of a persisting period doubling.
We perform the explicit derivation of the mapping in Appendix B and we find the effective bosonic Hamiltonian to bê with the constraint l m=−ln m = N and In the bosonic representation the initial state has the form It is very important to remark that here the bosons jump on a linear chain of length 2l + 1, while in the clock model they used to jump over a ring. This difference in topology makes impossible the realization of the period n-tupling of [18] using spin variables. We choose parameters where the classical model Eq. (1) shows period doubling and we study its fate for finite l in Fig. 1(a-c). Here we plot some examples of stroboscopic evolution of O(t) versus t/τ with t = nτ .
We see that, fixing l, O(t) oscillates. Especially interesting is the stroboscopic time t * when O(t) crosses 0 for the first time. If this time increases with the system size N , the period-doubling oscillations persist in the thermodynamic limit and there is a period doubling. If this time saturates with N , the period-doubling oscillations are a transient phenomenon and there is no period doubling. We plot t * versus N for the values of l we have considered in Fig. 1(d). For l = 1 and l = 3/2, t * saturates quite clearly with N . For l = 2 there is a sudden drop and also here there is no period doubling. We see from Fig. 1(d) that t * increases with l. This is entirely consistent with the fact that for l → ∞ the model tends to the classical limit of [36] where there is a period doubling and O(t) persists indefinitely for N → ∞.

A. Quantum chaos
We can study if this dynamics is regular or quantum chaotic. "Regular" means similar to an integrable model where the (classical or quantum) dynamics is constrained by as many local and commuting integrals of motion as degrees of freedom [45][46][47]. "Quantum chaotic" means that the Hamiltonian is equivalent to random matrix and this leads in general to thermalization of local observables [48][49][50]. In order to probe the regular or quantumchaotic behavior, we use the average level spacing ratio, defined as [51] r where µ α are the Floquet levels [52] and H is the relevant Hilbert subspace (more details below). The µ α are obtained from the eigenstates e −iµατ of the timeevolution operator over one periodÛ (τ, 0) of the Hamiltonian Eq. (2) and they are taken in increasing order [53]. If r 0.5269 the level-spacing distribution is of the COE type and the dynamics is ergodic (the Floquet states are like eigenstates of a random matrix) while if r 0.386 the level-spacing distribution is of the Poisson type and the model is integrable (see for instance [54]). We can evaluate r for the Hamiltonian in Eq. (7) provided we restrict to H, the subspace even under the mirror symmetry m → −m, which is an irreducible eigenspace of U (τ, 0) [55]. We can see that r reaches the quantumchaotic COE value for l = 1 and K ≥ 2, while for l = 2 the system shows always quantum chaos (see Fig. 2). This closely mirrors the classical-chaotic behaviour of the corresponding N → ∞ Gross-Pitaevskii equations observed through the Lyapunov exponent (see Sec. IV B 1).
In the next subsection we consider the limit N → ∞ and show that the model is described there by a system of Gross-Pitaevskii equations. In this case we will see persisting oscillations for O(t), for any l, and we will argue that they are Rabi oscillations between the states with angular momenta N l and −N l.
We can writeb m = √ Nβ m . We see that So, in the limit N → ∞, these are classical variables and have vanishing correlations. Using this fact, evaluating the expectation over the initial state of Eq. (10) (we define β(t) ≡ ψ(0)|β m, H (t)|ψ(0) ), and performing the limit N → ∞, we get the Gross-Pitaevskii equations with β m (t) ≡ 0 for m < −l or m > l. These equations are pretty simple to simulate numerically even for quite large values of l [56] and we do it using 4th order Runge-Kutta [57]. The initialization is β m (0) = δ m l . The expectation of the operatorŜ z /N (see Eq. (8)) is easily written as We show some examples of stroboscopic Gross-Pitaevskii evolution compared with the finite N cases in Fig. 3 (ac). For N → ∞, we see very clear Rabi oscillations of O(t) with no decay. These oscillations come from the resonance between the state with z angular momentum l (β m = δ l, m ) and the one with z angular momentum −l (β m = δ −l, m ). At finite N these states are |ψ ↑ = 1 and correspond to z angular mo-vii mentum N l and −N l, respectively. When K, h 1, we expect that these states are connected in perturbation theory at order ∼ 2l + 1, so the frequency ω Rabi of the Rabi oscillations of O(t) should be of order [17] . (14) From our numerics we find exactly this exponential scaling [see Fig. 3(d)]. We evaluate ω Rabi frequency by performing the Fourier transform of the signal of s z (t), finding the frequency ω peak corresponding to the maximum of the power spectrum and then evaluating ω Rabi = π − ω peak . The vanishing of ω Rabi for l → ∞ implies the existence of persisting period-doubling oscillations in this limit, which is equivalent to the classical case (see Sec. II). In agreement with that, for the parameters of Fig. 3(d), the classical case Eq. (1) shows period doubling, as one knows from Ref. [36]. We further remark that the Rabi oscillations for uncoupled spins (K = 0) in this same model have been already observed in Ref. [17].
We consider also the amplitude of the Rabi oscillations ∆O. We define them as square deviation of s z (t) [Eq. (13)] over time. We call it ∆O because it is also the mean square deviation of O(t), as it is easy to show. In order to make a comparison between different values of l possible, we consider ∆O/l. We plot this quantity versus l in Fig. 4(a). For every l, we see a crossing point between the curve for l and the one for l + 0.5. We see that the crossing moves towards the right for increasing l and for l = 2 the crossing is at K * ∼ 0.7. For K < K * the value of ∆O/l increases with l, for K > K * it decreases. This suggests that there is a phase transition in the limit l → ∞, as actually occurs [36]. Moreover, also the curves for ω Rabi versus K show a crossing [ Fig. 4(b)]. This crossing occurs for K = 1 and there is no contradiction with the result for the amplitude because in that case K * increases with increasing l and K * < 1. The crossing in the Rabi frequency is a strong evidence of a transition in the limit l → ∞ between a period-doubling and a trivial phase, and corresponds to what is observed in the dynamics of Eq. (1).

Largest Lyapunov exponent
We evaluate here the largest Lyapunov exponent, which is a probe of exponential divergence of nearby trajectories and therefore a probe of chaotic dynamics [58]. The largest Lyapunov exponent is approximated as λ(T ), a stroboscopic average over T periods tending to λ for T → ∞. We compute λ(T ) evaluating the rate of exponential increase in each period and averaging over periods. In practice, we consider two points in the phase space with distance d 0 , we evolve over a period and consider the value of the distance d 1 . Then we move the phase-space point of one of the trajectories along the segment joining the two so that we get again a distance d 0 1, and evolve again for one period getting a distance d 2 . Repeating T times, we get a sequence {d n } of distances [59] and we evaluate Taking T = 2 · 10 5 we already see convergence of λ(T ) and show the result in Fig. 5. What is remarkable is that this exponent is always positive, although it can get very small values (< 10 −2 ) for K < 1, marking thereby the existence of chaos. This classical chaos is fully mirrored by the quantum chaos occurring for finite N and appearing for any value of K if l is large enough (see Fig. 2). Only for l = 1 and K < 2 there is a lack of correspondence between the quantum behaviour (not quantum chaotic) and the classical nonvanishing Lyapunov exponent. Nevertheless, right at K = 2 the Lyapunov exponent shows a discontinuity mirroring thereby the crossover in the quantum finite-N behaviour. For small K, the system is chaotic but not ergodic. Indeed, it can support a regular behaviour as the one in Fig. 3. And we have checked that this behaviour is not due to an isolated regular trajectory: we see the same oscillations even if we take a slightly different initial state (β m (0) = δ m 0 + √ 1 − 2 δ m 1 ), see Fig. 6. Nevertheless, this is just a finite-time analysis and a chaotic behaviour might manifest at a time exponentially large in 1/K [60].
The largest Lyapunov exponent plotted in Fig. 5 allows to estimate the time scale over which the Gross-Pitaevskii description is valid for finite N . We see from Eq. (11) that for a finite-N system the width of the quantum fluctuations of β m (t) is at best ∼ 1/ √ N . Due to chaotic dynamics, this initial uncertainty increases exponentially in time with rate λ. The time the uncertainty reaches order 1 is After this time, the dynamics is quantum.

V. ANALYSIS OF MODEL-2
We get this model by applying the substitution Eq. (3) into Eq. (1) and then multiplying the resulting Hamiltonian by 2l. The resulting Hamiltonian is given in Eq. (23). Similarly to what we have done above, we de-fineŜ and, in order to understand if there is period doubling, we study the evolution of the period-doubling order parameter Eq. (6). Our initial state is given by We notice that for K = 0 this model reduces to the kicked Lipkin model of Ref. [17]. This model showed period doubling for h and φ − π small enough. In particular, O(t) showed Rabi oscillations with a frequency ∼ (h/J) 2l . In the limit l → ∞ (which in that context was the thermodynamic limit) the frequency of these oscillations tended to 0 and the period-doubling order parameter O(t)/l persisted to keep a finite value up to t → ∞. Now we couple many of these models with each other by means of the coupling K. As we have seen in the discussion for the model-1, which is equivalent to this one, this coupling is not strong enough to stabilize the order parameter to a value different from 0 for any finite l. At most, if K is small enough, the order parameter still shows Rabi oscillations with a renormalized frequency.
We study here the model-2 by means of the DTWA, an approximation which has proved to work fine in a long-range context [38][39][40]. We see that the DTWA is unable to reproduce the Rabi oscillations, but correctly gets the fact that there is no period doubling for finite l in the limit of large N . We get period doubling, in agreement with the exact dynamics, only in the classical l → ∞ limit. We briefly outline the DTWA approach in the next subsection.

A. Discrete truncated Wigner approximation in a nutshell
This is an approximation method especially convenient for long-range interacting spin models. All the details can be found in [38,39,61]. Here we just outline the application to our case. We start by expanding the expectation of a generic operatorB in a basis of operators in the form where w β ≡ 1 2 Tr Â βρ is the Wigner function, B w β (t) = Tr Â βB (t) are the Weyl symbols andB(t) ≡ U † (t, 0)ÔÛ (t, 0) withÛ (t, 0) the time-evolution operator form 0 to t of the Hamiltonian Eq. (23). We can take a basis of operators factorized over the siteŝ  where we can take over each site [62] A β = 1 + s β ·σ 2 (19) where s β can take the values 1 1 1 , −1 1 −1 , 1 −1 −1 and −1 −1 1 andσ = σ xσyσz . The approximation amounts to take the evolution ofÂ β as factorizedÛ The s α j, m, βj, m (t) have as initial values the ones given in Eq. (19), for the corresponding β j, m , and obey the evolution equationṡ where µνρ is the usual Ricci tensor, the elementary Poisson brackets are s µ j, m, βj, m , s ν i, m , β i, m = µ ν ρ δ i j δ m m s ρ j, m, βj, m and we have defined x In our case we can implement a Monte Carlo sampling procedure to approximate the sum of 4 N terms in Eq (17) in a numerically feasible way. Being the initial state Eq. (16) given by the density matrix in Eq. (17) we have that w β = 1/2 N for all the products of operators in Eq. (19) containing only β j, m = (1 − 1 − 1) and β j, m = (1 1 1). So, one can approximate the expectation of any operator with a Monte Carlo sampling of the uniform distribution w β , with the desired accuracy.
More specifically, we focus on the expectation and we evaluate it as the average over n r random initializations where each s j, βj is initialized with probability 1/2 in the condition 1 1 1 and probability 1/2 in the condition −1 −1 1 . Remarkably, the error bars do not scale with the system size, so this method is feasible also in the case of large systems [61]. The errorbars are evaluated as 1/ √ n r times the mean square deviation over randomness. In our analysis we have found that already for n r = 800 we have a satisfying convergence (see Appendix C). We are going to apply the DTWA method in the next subsection to study the period-doubling dynamics of the model-2.

B. Results
Consistently with the results found in the case of the model-1 (Sec. IV) we find here no period doubling. Indeed the period doubling order parameter O(t) decays to 0 in a finite time, independent of the system size N [see an example for l = 3/2 in Fig. 7(a)], for a set of parameters where the classical model Eq. (1) shows period doubling. As in the model-1, the limit l → ∞ corresponds to the classical case Eq. (1). We show this fact in Fig. 7(b) where we fix N and show the stroboscopic evolution of O(t) versus t/τ for different values of l. We qualitatively see that O(t) decay over a longer time as l increases. We plot for comparison also the stroboscopic evolution of (−1) t/τ 1 N j m z j (t) in the classical case Eq. (1). This quantity persists for an infinite time and O(t) tends to this curve when l → ∞. We notice that already for l = 3 the quantum dynamics is very near to the classical one, at least until t/τ = 4 · 10 3 .
Let us move to study the decay of the period-doubling order parameter in a more quantitative way. First of all, we plot O(t)/l versus t/τ with a logarithmic scale along the vertical axis [see Fig. 8(a)] and we see that O(t)/l decays exponentially in time. We can find the rate of this decay by fitting the curve of log O(t) versus t with a  Fig. 8(b). We see that δ decays with l as a power law, δ ∼ 1/l γ . Fitting the bilogarithmic plot with a straight line we find the decay exponent to be γ 2.51. So, extrapolating, we find that δ → 0 when l → ∞ and so in this limit the classical model and the period doubling are recovered. In Appendix C we discuss another method to estimate the decay time of O(t) which gives similar results.
We remark that the model-2 is a different representation of the model-1 for the chosen initialization, as we have discussed in Sec. II. DTWA therefore gives results which are not quantitatively correct (it does not catch any Rabi oscillation) and become so only in the limit l → ∞. Nevertheless, when l is finite, this approximation correctly gets the absence of persistent period doubling in the limit of large N .

VI. CONCLUSION
In conclusion we add quantum fluctuations to a classical and Hamiltonian model of interacting classical angular momenta showing synchronized period doubling in the thermodynamic limit. We consider the case of allto-all interactions where the long range correlations are stronger and the synchronized period doubling is most robust in the classical case. We study the robustness of synchronized period doubling adding quantum fluctuations in two different ways, realizing two different quantum models. In both the quantum models we find that the synchronized period doubling is fragile to quantum fluctuations and disappears. We perform our analysis by means of the so-called period doubling order parameter. When the system shows period doubling in the thermodynamic limit, the first zero of this order parameter occurs at at a time t * scaling to infinity for increasing system size. In both the quantum models t * does not increase with the system size, and so there is no period doubling, whenever the quantum fluctuations are significant.
We construct the quantum model-1 by replacing the classical angular momenta with quantum spins of size l. For any finite l there are quantum fluctuations and we show that the model becomes classical in the limit l → ∞. We restrict to the subspace even under all the permutation symmetries of the Hamiltonian, performing a mapping over a bosonic model. Due to the moderate Hilbert subspace dimension, we perform exact diagonalization for quite large system sizes and do the finite-size scaling of t * . For all the accessible values of l, we find no scaling, and so there is no period doubling.
This result is confirmed in the limit of infinite system size (N → ∞), where the bosonic model is described by a system of Gross-Pitaevskii equations. In this limit, the period-doubling order parameter performs Rabi oscillations related to the existence of resonant states in the spin model. For increasing l, these states are connected at higher orders in perturbation theory and consistently the frequency of the Rabi oscillations exponentially decreases in l. In particular, for l → ∞ the Rabi frequency goes to 0, so t * tends to infinity, and the classical period doubling is restored, consistently with the model becoming classical in this limit. Studying the dependence of the amplitude and the frequency of the Rabi oscillations on the parameter K, we find that the curves for different l cross at a point around K ∼ 1. This point corresponds to the transition from synchronized to unsynchronized behaviour in the classical l → ∞ limit.
For any finite N and the accessible values of l, we observe quantum chaos in this model, as shown by the average level spacing ratio being Wigner-Dyson. This is true for K 2 for l = 1 and for any value of K for l = 2. Analogously, in the N → ∞ limit, the Gross-Pitaevskii equations show a positive largest Lyapunov exponent λ (with a discontinuity at K = 2 for l = 1) and then the dynamics are chaotic. The Lyapunov exponent spans many orders of magnitude as K increases. In particular, for the values K < 1 corresponding to a synchronized period doubling in the limit l → ∞, we see λ ≤ 10 −2 for all the considered values of l. In this case we have a chaotic but not ergodic dynamics, as we verify by checking the existence of the Rabi oscillations of the period doubling order parameter also for a different initial condition. For any finite size N , we show that the dynamics is correctly described by the Gross-Pitaevskii equations up to a time scaling with log N .
Then we move to introducing the model-2. We substitute the classical angular momenta with sums of 2l Pauli matrices. We argue that also this model reduces to the classical one in the limit l → ∞. We find convenient to study this model by an approximation method called DTWA, which is known to give good results for longrange interacting models, and describe it in some detail. We focus on a set of parameters where the classical model shows period doubling and we use DTWA to study the evolution of the period-doubling order parameter. We find that it decays to zero as an exponential and the time scale of this decay does not scale with the system size, marking the absence of period doubling. Nevertheless, the time scale of the exponential decay scales as a power law with l. So, in the limit of l → ∞ there is period doubling, consistently with the model being classical in this limit. This model is equivalent to model-1 for the chosen initial conditions, so we see that DTWA provides quantitatively correct results only for l → ∞. For l finite xii it is not correct (it does not provide Rabi oscillations) but correctly predicts the absence of persistent period doubling in the limit of large N .
Therefore, we find that the period doubling in this model is fragile to the smallest quantum fluctuations. Prospects of future research include the exploration of Hamiltonian synchronization in different models, with stronger long-range correlations, and the study of its stability under quantum and thermal fluctuations.

ACKNOWLEDGMENTS
Part of the computational resources were provided by CINECA on the Marconi-A3 partition through the HPC agreement between CINECA and the ICTP. We acknowledge T. Mendes who provided us the access to these resources.
Appendix A: Limit l → ∞ as a classical limit Due to the equivalence of the two models for the chosen initial conditions, we discuss only the model-1. Conclusions apply also to the second one, not in general but for the chosen initial conditions. The analysis strictly resembles the one leading to the Gross-Pitaevskii for the bosonic model in the N → ∞, as we have discussed in Sec. IV B. We rescale the spin variables asŜ α j = 1 lŝ α j . Their commutator is [Ŝ α j ,Ŝ β i ] = i l αβγŜ γ j δ i j , so these variables are classical in the limit l → ∞. One can write the Heisenberg equations for theŜ α j variables d dtŜ then evaluate the expectation over the initial state and, performing the limit l → ∞, neglect any quantum correlation due to the vanishing commutator. Performing this calculation one gets l → ∞ evolution equations for the expectations ofŜ α j which exactly coincide with the classical evolution equations for m α j (t) obtained with Eq. (1) using the classical Poisson brackets.
There is another method, nearer to the analysis of [44], which we are going to sketch. Take for instance the operatorŝ x j and apply it to the state |l, m j , the eigenstate with eigenvalue l(l +1) ofŝ 2 j and eigenvalue m ofŝ z j . One gets [42] Defining q j = m/l and using the translation operator of shift a along q j , exp(a d dqj ), one finds for l 1 Definingp j = − i l d dqj , we see that we have written this object in terms of two canonical variables,q j andp j whose commutator is [q j ,p j ] = i/l. In the limit l → ∞, therefore, they are classical canonical variables obeying the canonical Poisson bracket {q j , p j } = −1. The classical m α j can then be obtained as m x j = b 1 − q 2 j cos p j , with b > 0 arbitrary giving the size of the classical spin. In the same way one gets [44] m z j = bq 2 j and m y j = b 1 − q 2 j sin p j . Appropriately fixing b = 1/2 one gets the classical Hamiltonian Eq. (1). The canonical Poisson brackets of q j and p j give rise to the angular momentum Poisson brackets for the m α j which are stated immediately below Eq. (1). In this way one gets back the classical dynamics for l → ∞.

Appendix B: Mapping onto the bosonic model
We can now discuss the bosonic mapping of the notes. This mapping was introduced in [18,37] for similar infinite-range models. Let us consider a system of N sites, and let us take the local spins with value l = 1 for clarity (the generic case is exactly identical). With l = 1, we can have m = −1, 0, 1. Because the system is fully symmetric under permutations, we can restrict to the states even under all the possible N ! permutations. If we callP the sum of all the permutation operators, we can take as basis of our Hilbert space the states There are N sites, so n −1 + n 0 + n 1 = N . The factor in front is for normalization, and the √ N ! is there because there are N ! possible permutations. The factors √ n m ! at the denominator are there for the following reason. Consider for instance m = 1. Fixing everything else, there are n 1 ! ways of rearrange the sites with m = 1 and this increases the norm by a factor n 1 !. One divides by √ n 1 ! and the norm is again 1. Consider for instance the application of the operator S + = jŝ + j , whereŝ + j =ŝ x j + iŝ y j . One findŝ S + |n −1 n 0 n 1 = 1 N ! (n −1 !n 0 !n 1 !)P jŝ If for instanceŝ + j acts over a site with m = −1, one gets a factor 1(1 + 1) = √ 2 [42]. Moreover, there are n −1 of these sites and they are equivalent due to the permutation operator. This gives rise to a factor n −1 in xiii front. Moreover, in this way one decreases n −1 by 1 and increases n 0 by 1. One has a similar situation for m = 0 and m = 1, sô S + |n −1 n 0 n 1 = 1 N ! (n −1 !n 0 !n 1 !) n −1 √ 2P |(−1 . . . The exponential decay found in Fig. 8(a) does not last forever and at some point the period-doubling order parameter starts oscillating around 0, as we have seen in Fig. 7(a). Let us call t * the first value of the stroboscopic time where the order parameter vanishes. We define . (C1) As we can see in Fig. 9(a), t d increases with l as a power law, t d ∼ l γ . From a linear fit of the bilogarithmic plot we find γ 2.20, in perfect agreement with the finding of Sec. V B. The errorbars for t d come from the errorbars for t * . We evaluate the latter from the errorbar in O(nτ ) (evaluated as described in Sec. V A) which gives rise to an error in the time of first vanishing t * . In Fig. 9(b) we plot t d versus n r . In the cases where we can numerically afford n r > 800, we see that for n r = 800 we have already attained convergence. For larger values of l we cannot go beyond that value, but the clear scaling with l suggests that a satisfying convergence has been already attained for this value of n r .