Chaos and subdiffusion in the infinite-range coupled quantum kicked rotors
CChaos and subdiffusion in the infinite-range coupled quantum kicked rotors
Angelo Russomanno, Michele Fava, and Rosario Fazio
3, 4 Max-Planck-Institut f¨ur Physik Komplexer Systeme,N¨othnitzer Straße 38, D-01187, Dresden, Germany Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory,University of Oxford, Oxford OX1 3PU, United Kingdom Abdus Salam ICTP, Strada Costiera 11, I-34151 Trieste, Italy Dipartimento di Fisica, Universit`a di Napoli “Federico II”, Monte S. Angelo, I-80126 Napoli, Italy
We map the infinite-range coupled quantum kicked rotors over a infinite-range coupled interactingbosonic model. In this way we can apply exact diagonalization up to quite large system sizes andconfirm that the system tends to ergodicity in the large-size limit. In the thermodynamic limitthe system is described by a set of coupled Gross-Pitaevskij equations equivalent to an effectivenonlinear single-rotor Hamiltonian. These equations give rise to a power-law increase in time of theenergy with exponent γ ∼ / I. INTRODUCTION
The kicked rotor is a paradigmatic model in classi-cal Hamiltonian dynamics. This simple model has beenwidely used to numerically study the development ofchaos when integrability is broken, in accordance withthe KAM theorem [1–3]. In this single-degree-of-freedommodel there can be weak chaos when the integrabilitybreaking is small and there are still many manifolds ofintegrable-like dynamics (conserved tori) in phase space.In this case of weak chaos, the energy does not increasein time, while in the case of strong chaos the dynam-ics is ergodic and diffusive in phase space and the en-ergy linearly increases in time. Considering the case ofmany coupled kicked rotors, on the opposite, there arenot enough conserved tori and the dynamics eventuallyshows a diffusive behaviour with energy linearly increas-ing in time [1, 2]. Nevertheless, the time scale after whichthis behaviour appears is very long and exponentiallylarge with the inverse of the perturbation from integra-bility [4, 5]. This is a particular case of a general theoremdue to Nekhoroshev [6]. So, a classical Hamiltonian sys-tem slightly perturbed from integrability shows a thermalergodic behaviour after a long non-thermal transient, instrict analogy with the well-known quantum prethermal-ization [7].The single kicked rotor is very interesting also fromthe quantum point of view. Indeed this model showsa lack of correspondence between classical and quantumbehaviour. Differently from other cases [8–14], the kickedrotor shows a diffusive dynamics with linearly increas-ing energy in the strong classical regime which has nocounterpart in the quantum domain. When the model isquantized due to quantum interference the classical en-ergy increase stops after a transient. This behaviour isknown as quantum dynamical localization [15–17], and has been interpreted as a dynamical disorder-free ana-log of Anderson localization [18–20] and has also beenexperimentally observed [21, 22].It is interesting to understand if dynamical localiza-tion persists in the case of many interacting rotors. Ithas been argued for different types of interactions, thatfor any set of parameters there is a threshold in systemsize beyond which the system becomes ergodic and theenergy increases in time without a bound [23, 24]. So, incontrast with other models [25–27], there is no dynam-ical localization in the many-body limit. Nevertheless,the system does not become equivalent to the classicalmodel.This is clear in the case of infinite-range coupling wherethe thermodynamic-limit dynamics can be exactly com-puted [24] and one finds that the energy increases asa power law with exponent smaller than 1 (subdiffu-sively): Quantum effects make the classical linear energyincrease slower. A power-law increase of the energy hasalso been observed in the quantum few-rotor interactingcase [28] and in many non-linear generalizations of thesingle quantum kicked rotor [29–31] and related disor-dered models [32–34]. On the opposite, in the classicalchaotic many-rotor case, the energy increases linearly intime [4, 5]. One can observe a linear increase of the en-ergy also in a single quantum kicked rotor if the kickingamplitude is modulated by a noise [35].In this work we focus again on the quantum infinite-range coupled quantum kicked rotors. We restrict to thesubspace even under all the site-permutation transforma-tions, we apply a method used before in another infinite-range coupled context [36], and map the model over abosonic infinite-range interacting model over a lattice.From one side this fact allows us to apply exact diago-nalization for larger system sizes and larger truncationswith respect to what previously done. In this way we can a r X i v : . [ qu a n t - ph ] F e b probe ergodicity by means of the average level spacingratio for larger system sizes N and find further confirma-tion for the generalized tendency towards ergodicity forincreasing N previously demonstrated [24].This bosonic mapping is convenient also because allowsto show that, in the limit N → ∞ , the model is describedby a system of Gross-Pitaevskij equations. In the limit ofvanishing interactions, these equations are equivalent tothe Schr¨odinger equation of the single kicked rotor rep-resented in the momentum basis. They are exact in thelimit N → ∞ , and we can show that in general they areequivalent to the Schr¨odinger dynamics of the non-linearsingle kicked rotor effective Hamiltonian found in [24].The energy increases in time as a power law with expo-nent γ < γ ∼ / /
3. Moreover, we predict that the coarse-grainedsquared non-linear modulation of the kick depends ontime as t − / and we numerically verify this fact. Ourapproach can be applied also to the case of a kicking mod-ulated by a noise with properties invariant under timetranslations. This gives rise to a diffusive behaviour ofthe energy, in agreement with [35].Our master-equation approach is possible thanks tothe chaotic behaviour of the Gross-Pitaevskij equations.In order to probe the chaoticity properties of these equa-tions, we evaluate the largest Lyapunov exponent, whichgives a measure of the rate of exponential divergence ofnearby trajectories [39]. We see that, when we considerparts of the phase space with larger and larger momen-tum, the largest Lyapunov exponent decreases as a powerlaw towards 0. So, for increasing momentum, the systemis still chaotic but becomes asymptotically regular in thelimit of infinite momentum. This is in agreement withthe results of the average level spacing ratio at finite sizesuggesting full ergodicity in the large-size limit.The largest Lyapunov exponent results are relevantalso for finite N (cid:29)
1. Here, due to Heisenberg principle,the relevant dynamical variables have an initial uncer-tainty of order 1 / √ N . For finite N we can use exponen-tial divergence of nearby trajectories to show that theGross-Pitaevskij equations are valid for a time increas-ing as log N , with a coefficient equal to the inverse ofthe relevant Lyapunov exponent. So, in the limit of infi-nite momentum, the Gross-Pitaevskij equations tend tobe valid for an infinite time, as one naively would expect, being in that limit the integrability-breaking interactionterm vanishingly small compared with the undriven term.Moreover, the validity of the Gross-Pitaevskij equationsfor a time logarithmic in N in the case of a chaotic dy-namics is a fact generally valid in bosonic mean field dy-namics. This is relevant for instance in the context oftime crystals [40].Chaos in this model obeys the Nekhoroshev theo-rem [6]. In the noninteracting case the dynamics isthe Schr¨odinger dynamics of the single rotor which isintegrable and the resulting constants of motion of thestroboscopic dynamics are related to the Floquet states.Turning on the interaction, the Schr¨odinger dynamics be-comes the nonlinear Gross-Pitaevskij one, integrabilityis broken and the Floquet constants of motion are nomore conserved. They deviate in time from their initialvalue over a time scale exponential in the inverse of theintegrability-breaking interaction term.The paper is organized as follows. In Section II we in-troduce the model Hamiltonian and we perform the map-ping to the infinite-range bosonic model. We describe allthe details of the construction of this bosonic representa-tion in Appendix A. In Sec. II A we perform exact diago-nalization on this Hamiltonian, with an appropriate trun-cation, and by means of the average level spacing ratio weshow the existence of a generalized tendency to ergodicityfor increasing system size. In Sec. III we show that in thelimit N → ∞ the dynamics of the bosonic Hamiltonianis described by a system of Gross-Pitaevskij equationscompletely equivalent to a non-linear single-rotor effec-tive Hamiltonian. In Sec. IV we use these equations tonumerically study the evolution of the energy. We findthat it increases in time with a power law. We analyt-ically explain the power-law exponent γ = 2 / II. HAMILTONIAN AND MAPPING TO THEBOSONIC MODEL
So, this is the quantum infinite-range coupled kickedrotor ˆ H ( t ) = 12 N (cid:88) l =1 ˆ p l + V (ˆ θ ) δ ( t ) ; (1)where we define [41] δ ( t ) ≡ (cid:80) + ∞ n = −∞ δ ( t − n ) and V (ˆ θ ) ≡ ¯ kK N (cid:88) l =1 cos ˆ θ l − ¯ kK(cid:15) N − (cid:88) l, l (cid:48) cos(ˆ θ l − ˆ θ l (cid:48) ) (2)(here we have replaced the physical coupling K with ¯ kK ,in order to somewhat simplify the subsequent formulae).The commutation relations[ˆ θ l , ˆ p l (cid:48) ] = i ¯ kδ l, l (cid:48) (3)are valid, with ¯ k related to the physical parameters ofthe Hamiltonian (one arrives at Eq. (1) after an appro-priate rescaling [24]). In all the paper we will focuson the stroboscopic dynamics, looking at the system inthe instant n + , that’s to say immediately after the n th kick (in the text we will omit the superscript + ). Mostimportantly this model is invariant under all the site-permutation transformations. The subspace even underall these transformations is therefore an eigenspace of theHamiltonian, a point which will be crucial in the nextsection.We restrict to the subspace even under all the per-mutation transformations. Using the methods explainedin [36] we get the effective Hamiltonian (see Appendix A)ˆ H b ( t ) = ¯ k ∞ (cid:88) m = −∞ m ˆ n m + ¯ kK (cid:104) ∞ (cid:88) m = −∞ (cid:16) ˆ b † m ˆ b m +1 + H . c . (cid:17) − (cid:15) N − ∞ (cid:88) m,m (cid:48) = −∞ (ˆ b † m +1 ˆ b m ˆ b † m (cid:48) ˆ b m (cid:48) +1 + H . c . ) (cid:105) δ ( t ) , (4)where ˆ b m are bosonic operators obeying the commutationrelations [ˆ b m , ˆ b † m (cid:48) ] = δ m, m (cid:48) , [ˆ b m , ˆ b m (cid:48) ] = 0 and we defineˆ n m ≡ ˆ b † m ˆ b m . The bosons obey the constraint (cid:80) m ˆ n m = N . We notice that for (cid:15) = 0 (no interactions betweenthe rotors) the model reduces to a quadratic integrablebosonic model. This is in agreement with the knownfact [16] that the quantum kicked rotor behaves as anintegrable model with a Poisson level statistics [24]. Inthe interacting many-rotor case there can be an ergodicbehaviour and the level statistics behaves differently aswe are going to elucidate in the next section. A. Average level spacing ratio
This bosonic representation is quite convenient froma technical point of view because it is possible to ap-ply exact diagonalization to the Hamiltonian Eq. (4) forsystem sizes and truncations of the Hilbert space signif-icantly larger than those considered in [24]. A very im-portant object in a periodically-driven dynamics is thetime-evolution operator over one period which for theHamiltonian in Eq. (1) isˆ U ≡ e − iV (ˆ θ ) / ¯ k e − i (cid:80) Nl =1 ˆ p l / (2¯ k ) . We report the expression of V (ˆ θ ) and (cid:80) Nl =1 ˆ p l in thebosonic representation in Eq. (A9). We further restrict tothe subspace even under the m → − m symmetry [42]. Inthis subspace, we compute ˆ U , diagonalize it, and get the many-body Floquet quasienergies µ α as the phases of theeigenvalues of ˆ U [43, 44]. Of course, this is possible onlyimposing a truncation to the Hilbert space, restrictingto the states for which | m | ≤ M . In order to probeergodicity of the system, we can evaluate the averagelevel spacing ratio r [45] defined as r = 1 N M − N M − (cid:88) α =1 min( µ α +2 − µ α +1 , µ α +1 − µ α )max( µ α +2 − µ α +1 , µ α +1 − µ α ) (5)where the quasienergies are restricted to the first Flo-quet Brillouin zone (see for instance [46]) and taken inincreasing order. N M is the dimension of the truncatedHilbert space. It is important that on the subspace evenunder all the permutations we have imposed the furtherconstraint of being even under the m → − m symmetry.In this way we are restricting to an irreducible invari-ant subspace of the Hamiltonian, a condition requiredin order that the level spacing distribution (and the re-lated ratio r ) is a meaningful ergodicity indicator [47].When the driven system is ergodic, i.e. locally thermal-izing with T = ∞ [48], the Floquet operator ˆ U belongsto the circular orthogonal ensemble (COE) of symmet-ric unitary matrices (because of the time-reversal invari-ance) [28, 49–51]. In this case, the level spacing distri-bution is of the COE type and the average level spacingratio acquires the value r COE (cid:39) . r P (cid:39) . r versus K for a fixed (cid:15) in Fig. 1. For each N we choose M so that r has attained convergence (fixing N , we see a quite fast convergence of r with the trunca-tion M of the Hilbert space). We see that from a certain K on, r attains the COE value: The system becomesergodic. For increasing N the values of r genericallyincreases and r attains the COE value for smaller andsmaller values of K . This suggests a tendency towardsergodicity for increasing size. This is in agreement withthe analytical predictions of [24] of a completely ergodicsystem in the limit N → ∞ . III. GROSS-PITAEVSKIJ EQUATIONS IN THE N → ∞ LIMITA. Derivation
We start by considering the limit N → ∞ where thedynamics is described by effective classical equations.Everything is based on the definition ˆ β m ≡ ˆ b m / √ N andon the observation that[ ˆ β m , ˆ β † m (cid:48) ] = δ m, m (cid:48) /N , (6) r K N = 3, M = 18N = 4, M = 12N = 5, M = 9N = 6, M = 7r P r COE
FIG. 1. r versus K for different values of N . Numericalparameters: K = 6, (cid:15) = 0 .
52, ¯ k = 1 . so in the limit N → ∞ , the ˆ β m behave as classical vari-ables. We write the Heisenberg equations for the ˆ β m i ¯ k dd t ˆ β m H ( t ) = [ ˆ β m H ( t ) , ˆ H b ( t )]and exploit the fact that in the limit N → ∞ the ˆ β m areuncorrelated, due to the vanishing commutators in thislimit [Eq. (6)]. By defining β m ( t ) = (cid:104) ψ ( t ) | ˆ β m | ψ ( t ) (cid:105) , inthe N → ∞ limit we get from the Heisenberg equations the following system of Gross-Pitaevskij equations i dd t β m ( t ) = ¯ k m β m ( t ) + K δ ( t ) (cid:20) β m +1 ( t ) + β m − ( t ) − (cid:15) (cid:16) (cid:88) m (cid:48) β ∗ m (cid:48) ( t ) β m (cid:48) +1 ( t ) (cid:17) β m − ( t ) − (cid:15) (cid:16) (cid:88) m (cid:48) β m (cid:48) ( t ) β ∗ m (cid:48) +1 ( t ) (cid:17) β m +1 ( t ) (cid:21) . (7)The energy per site of the unkicked part of the model canbe written as e ( t ) = ¯ k (cid:88) m m | β m ( t ) | . (8)We see that for (cid:15) = 0 this system of equations is equiva-lent to the Schr¨odinger equation of a single rotor [50].The Gross-Pitaevskij equations can be obtained fromthe effective classical Hamiltonian H ( { β m , β ∗ m } , t ) = ∞ (cid:88) m = −∞ ¯ k m | β m ( t ) | + K (cid:104) ∞ (cid:88) m = −∞ ( β ∗ m β m +1 + H . c . ) − (cid:15) ∞ (cid:88) m,m (cid:48) = −∞ ( β ∗ m +1 β m β ∗ m (cid:48) β m (cid:48) +1 + H . c . ) (cid:105) δ ( t ) . (9) In order to get the Gross-Pitaevskij equations, one writes dd t β m = −{ β m , H} where {· · · } are the Poisson brack-ets [53], and the elementary Poisson brackets (needed toevaluate all the other ones) are { β m , β ∗ n } = iδ m n and { β m , β n } = 0 ∀ n, m . (10)These Poisson brackets will have a relevant role in theanalysis related to the Nekhoroshev theorem in Sec. VI.In the limit M → ∞ we can actually map the semiclas-sical model Eq. (7) into the self-consistent single-particlemodel studied in [24]. In order to do that, we rewriteEq. (7) in a sort of continuum-limit approximation. Iden-tifying x ≡ m and ψ ( x ) ≡ β m , expanding to all orders in1 /m as in [54], we get i ¯ k dd t ψ ( x, t ) = ¯ k x ψ ( x, t ) + ¯ kK δ ( t ) (cid:20) (cid:18) cosh (cid:18) dd x (cid:19) − (cid:19) ψ ( x, t ) − (cid:15) (cid:18) (cid:90) ∞−∞ ψ ∗ ( x (cid:48) , t ) exp (cid:18) dd x (cid:19) ψ ( x (cid:48) , t )d x (cid:48) (cid:19) exp (cid:18) − dd x (cid:19) ψ ( x, t ) − (cid:15) (cid:18) (cid:90) ∞−∞ ψ ∗ ( x (cid:48) , t ) exp (cid:18) dd x (cid:19) ψ ( x (cid:48) , t )d x (cid:48) (cid:19) ∗ exp (cid:18) dd x (cid:19) ψ ( x, t ) (cid:21) , (11)where we have exploited that (cid:82) ∞−∞ | ψ ( x, t ) | = 1 and usedthe formal relation ψ ( x +1) = exp (cid:0) dd x (cid:1) ψ ( x ), involving allthe perturbative orders [54]. Introducing the appropriate coordinate and momentum operatorsˆ θ = i dd x ˆ p = ¯ kx , (12)such that [ˆ θ, ˆ p ] = i ¯ k , we can rewrite this formula as i ¯ k dd t | ψ t (cid:105) = ˆ H eff ( t ) | ψ t (cid:105) . (13)with the effective self-consistent Hamiltonianˆ H eff ( t ) ≡ ˆ p kKδ ( t ) (cid:20) cos(ˆ θ ) − (cid:15) (cid:104) ψ t | cos ˆ θ | ψ t (cid:105) cos ˆ θ − (cid:15) (cid:104) ψ t | sin ˆ θ | ψ t (cid:105) sin ˆ θ (cid:21) . (14)This is exactly the effective mean-field Hamiltonian fol-lowing Eq. (1), found in [24] through a much more in-voluted analysis. We can easily see that, if the ini-tial state is symmetric under ˆ θ (cid:55)→ − ˆ θ , this symme-try will be preserved during the time-evolution, so that (cid:104) ψ t | sin ˆ θ | ψ t (cid:105) = 0. Consequently F ( t ) = (cid:88) m β ∗ m ( t ) β m +1 ( t ) = (cid:104) ψ t | cos ˆ θ | ψ t (cid:105)− i (cid:104) ψ t | sin ˆ θ | ψ t (cid:105) (15)is purely real. Thus, the effective self-consistent Hamil-tonian Eq. (16) acquires the formˆ H eff ( t ) ≡ ˆ p kKδ ( t ) [1 − (cid:15)F ( t )] cos ˆ θ . (16)This effective Hamiltonian will be very relevant for us inthe next section, where we will use it as the starting pointto construct a master-equation model for explaining theenergy subdiffusion with exponent 2 / IV. ENERGY SUBDIFFUSION
We numerically solve equations (7) by means of afourth order Runge-Kutta method with adaptive timestep scheme [55, 56]. We initialize them in the sym-metrized m = 0 state for the rotors, therefore we choose β m (0) = (cid:26) m = 00 otherwise . and we can compare the re-sults with those of the mean field analysis performed ina different way in [24]. In order to implement them wehave to impose a truncation, fixing some M > m such that m ≤ M . We showsome plots of e ( t ) versus t with stroboscopic time ( t isinteger) for an interacting case with (cid:15) (cid:54) = 0 in Fig. 2. Thehorizontal line corresponds to the T = ∞ value of the un-kicked energy in the truncated subspace. This quantitycan be readily evaluated as e T = ∞ ( M ) = ¯ k M + 1 N (cid:88) m =1 m = ¯ k M + 1 (cid:18) M M M (cid:19) . (17)In Fig. 2 we can see convergence towards e T = ∞ , andtherefore thermalization.
10 100 1000 10000 100000 1 10 100 1000 10000 e ( t ) t M = 8M = 13M = 28 M = 53M = 68M = 103 M = 253
FIG. 2. Stroboscopic evolution of e ( t ) versus t with trunca-tion (bilogarithmic plot). We plot also the T = ∞ value in thetruncated subspace (horizontal lines). We first see a power-law increase and then a convergence to the T = ∞ value.Without a truncation the power-law increase would last for-ever. Numerical parameters: K = 6, (cid:15) = 0 .
52, ¯ k = 1 . Most importantly, we see a power-law increase in timeof the energy, still in agreement with the findings of [24].We can see this power-law increase until saturation setsin, but with M → ∞ it would last forever, as one can eas-ily convince by looking at Fig. 2. We show some examplesof power-law increase for different choices of parametersin Fig. 3. The curves with (cid:15) ∈ [0 . ,
1] are consistent atlong times with a power-law increase of the form e ( t ) ∼ t γ with γ ∼ / . This numerical result was already found in [24] (seeFigs. 9 and 10 of that paper). So, in a wide range of (cid:15) , the exponent of the power law appears to be indepen-dent of the choice of the parameters. In the following weprovide an analytical argument for explaining this find-ing.Before deriving the equation giving rise to momentumsubdiffusion, we explain the rationale behind our anal-ysis. The key observation is that energy absorption ispurely driven by the fluctuation of the self-consistent field F ( t ), in fact, if F ( t ) were constant, the system would bein a dynamical-localized phase with an asymptoticallyconstant energy. We can then expect that the evolutionwith the effective Hamiltonian Eq. (14) gives a ddt (cid:104) ˆ p (cid:105) t proportional to the variance of F ( t ) over time — as wewill show in Sec. IV B. Furthermore, as the momentum-range m t = (cid:112) (cid:104) ˆ p (cid:105) t / ¯ k grows in time, the variance of F ( t ) decreases towards 0 as1 /m t (see Sec. IV A). As a consequence the diffusion inmomentum is slowed down as the wave-function spreads,giving rise to the subdiffusive behaviour observed above.
10 100 1000 10000 100000 1x10
1 10 100 1000 10000 100000 e ( t ) tK = 4 ε = 0.01 ε = 0.03 ε = 0.17 ε = 0.52 ε = 1.0t
10 100 1000 10000 100000 1x10
1 10 100 1000 10000 100000 e ( t ) tK = 6 ε = 0.01 ε = 0.03 ε = 0.17 ε = 0.52 ε = 1.0t
100 1000 10000 100000 1x10
1 10 100 1000 10000 100000 e ( t ) tK = 8 ε = 0.01 ε = 0.03 ε = 0.2 ε = 0.52 ε = 1.0t FIG. 3. Stroboscopic evolution of e ( t ) versus t with truncation M = 20153 for different values of the parameters (bilogarith-mic plot). Notice the power-law increase which at long timestends towards a behaviour consistent with the law e ( t ) ∼ t / when (cid:15) ∈ [0 . , k = 1 . A. Statistical properties of F ( t ) While the non-linear equation Eq. (7) is hard to charac-terize analytically, a simplifying assumption is to consider F ( t ) as an effective noise. This is reasonable because, as a consequence of the chaotic diynamics in Eq. (7), F ( t )shows random oscillations symmetric around 0, as wecan see in Fig. 4. Before starting our analysis we discusswhich minimal set of self-consistency properties shouldthe effective noise have in order to adequately describe F ( t ).Firstly, looking at Eq. (15), we can split the sums over m into pieces, each of one including O ( ξ loc ) adjacent m s,where ξ loc is the single-rotor localization length (in mo-mentum space). It is now reasonable to assume that thecontribution of each of these pieces will be largely un-correlated with the other ones. Then, when the state isspread over a momentum range m t (cid:29) ξ loc , by the cen-tral limit theorem, the sum F ( t ) will be distributed likea Gaussian with mean zero.To estimate the variance of the Gaussian noise, we usethe normalization condition (cid:80) m | β m | = 1 to infer thatthe magnitude of | β m | scales like 1 / √ m t for m (cid:46) m t , andis approximately zero otherwise. Then summing over aregion of length 2 ξ loc around ¯ m , such that | ¯ m | (cid:46) m t , weapproximately haveVar (cid:88) | m − ¯ m | <ξ loc β ∗ m ( t ) β m +1 ( t ) ∼ (cid:18) ξ loc m t (cid:19) (18)assume perfect correlation inside that region. Consider-ing, instead, pieces from different regions as uncorrelated,we have thatVar [ F ( t )] ∼ m t ξ loc (cid:18) ξ loc m t (cid:19) ∼ ξ loc m t (19)Estimating the localization length in terms of the single-rotor parameter (see e.g. Ref. [20]) we then haveVar [ F ( t )] ∼ K m t . (20)Finally, since the evolution of β is chaotic (see Sec. V),we will assume that the autocorrelation of F ( t ) decays intime over a finite scale τ .These assumptions are consistent with the numericalsolution of the equation of motion. We show an examplein Fig. 4, where we can see that the random oscillationsare symmetric around 0 and F ( t ) is short-range corre-lated in time. B. Effective master equation
We start by considering the self-consistent Hamilto-nian Eq. (16) and we define ˆ H ( t ) ≡ ˆ p + ¯ kKδ ( t ) cos ˆ θ ,ˆ V ≡ cos ˆ θ . First of all we move to the interaction rep-resentation defining ˆ U ( t ) ≡ ←− exp (cid:16) − ( i/ ¯ k ) (cid:82) t dt (cid:48) ˆ H ( t (cid:48) ) (cid:17) ,and ˆ O I ( t ) ≡ ˆ U † ( t ) ˆ O ˆ U ( t ), | ψ ( t ) (cid:105) I ≡ ˆ U † ( t ) | ψ ( t ) (cid:105) for ageneric operator ˆ O . In this representation Eq. (13) takes -0.3-0.2-0.1 0 0.1 0.2 0.3 0 20000 40000 60000 80000 100000 F ( t ) t FIG. 4. Example of evolution of F ( t ). Notice the chaoticbehaviour symmetric around 0. Numerical parameters: K =4, (cid:15) = 0 .
52, ¯ k = 1 . M = 20153. the form ddt | ψ ( t ) (cid:105) I = − i(cid:15)Kδ ( t ) F ( t ) ˆ V I ( t ) | ψ ( t ) (cid:105) I . (21) Treating F ( t ) as a Gaussian noise, we can derive fromEq. (21) a master equation similar to the ones resultingfrom the coupling to an environment [37]. We write firstthe evolution equation for the density matrix ˆ ρ ψI ( t ) ≡| ψ ( t ) (cid:105) I (cid:104) ψ ( t ) | i dd t ˆ ρ ψI ( t ) = − (cid:15)Kδ ( t ) F ( t )[ ˆ V I ( t ) , ˆ ρ ψI ( t )] . (22)Integrating Eq. (22) once we getˆ ρ ψI ( t ) = ˆ ρ ψI (0) + i(cid:15)K (cid:90) t δ ( t (cid:48) ) F ( t (cid:48) )[ ˆ V I ( t (cid:48) ) , ˆ ρ ψI ( t (cid:48) )]d t (cid:48) . (23)Substituting this formula into Eq. (22) and averagingover the Gaussian noise (the assumption of Gaussianityallows to neglect the correlations between ˆ V I and ˆ ρ ψI ) wegetdd t ˆ ρ I ( t ) = − (cid:15) K (cid:90) t (cid:104) F ( t (cid:48) ) F ( t ) (cid:105) δ ( t ) δ ( t (cid:48) ) (cid:104) ˆ V I ( t ) , [ ˆ V I ( t (cid:48) ) , ˆ ρ I ( t (cid:48) )] (cid:105) d t (cid:48) , (24)where (cid:104)· · ·(cid:105) marks the average over the Gaussian ensem-ble and we have defined ˆ ρ I ( t ) = (cid:104) ˆ ρ ψI ( t ) (cid:105) . At this point weapply a coarse-graining average over a mesoscopic time-scale ∆ t (cid:29) τ,
1. Furthermore, we assume that averagesof functions of F ( t ) over time can be substituted withthe corresponding averages over the Gaussian ensembleproducing the noise (see Sec. IV A). We expect this stepto be valid as long as the system is effectively ergodic onthe mesoscopic time-scale ∆ t , which is reasonable in thiscase, given that the system is chaotic (see Sec. V). Thus,we have (cid:104) F ( t (cid:48) ) F ( t ) (cid:105) (cid:39) F ( t ) g ( t − t (cid:48) ) (25)where ( · · · )( t ) marks the average over the coarse-grainingtime interval [ t − ∆ t, t + ∆ t ] and g ( t − t (cid:48) ) is a functionwhich decays to 0 when | t − t (cid:48) | > τ . Integrating over d t (cid:48) and applying the coarse-graining average, the δ func-tions disappear and give rise to an average over timesmultiple of the kicking period (stroboscopic times). Thekicking period 1 is much smaller than ∆ t , the resolutionin time after coarse graining, then we can approximatethe average over stroboscopic times with an average overcontinuous times. Moreover, assuming τ (cid:46) ∆ t , we cansubstitute t (cid:48) with t because the function g ( t − t (cid:48) ) acts asa delta function at the level of time resolution given by ∆ t . Defining σ ( t ) ≡ g(cid:15) K F ( t ) , (26)where the coefficient g results from the integration of g ( t − t (cid:48) ) and putting everything together, we getdd t ˆ ρ I ( t ) = − σ ( t ) (cid:104) ˆ V I ( t ) , [ ˆ V I ( t ) , ˆ ρ I ( t )] (cid:105) , (27)where, due to the presence of the δ ( t ), the coarse-graining average results being over discrete values of time( t is integer) and we have approximated the coarse-grained derivative as ˆ ρ I ( t +∆ t ) − ˆ ρ I ( t − ∆ t )2∆ t (cid:39) dd t ˆ ρ I ( t ) becausewe are interested in phenomena occurring over timescaleslarger than ∆ t . Moreover, σ ( t ) is already coarse grainedand the coarse-graining average does not affect it.Finally, since we are interested in the asymptotic prop-erties at large times, we can assume that m t (cid:29)
1. Conse-quently, by Eq. (20), we can assume that σ ( t ) (cid:28)
1. Thisassumption is crucial for our analysis since it allows usto make a separation-of-timescales approximation [37].With this assumption, we can see from Eq. (27) that thenoise induces significant changes in ˆ ρ I ( t ) over a time scale t r much larger than the typical order-1 timescale t of thedynamics induced by ˆ H . So, if we take t (cid:28) ∆ t (cid:28) t r ,ˆ ρ I ( t ) in Eq. (27) is not affected by the coarse-grainingaverage. In essence, ˆ ρ I ( t ) can be approximated by itscoarse-grained average, ˆ ρ I ( t ) (cid:39) ˆ ρ I ( t ). This fact will haveimportant consequences in the following analysis of sub-diffusion. C. Subdiffusion
In order to derive subdiffusion, we need to considerthe dynamics of the operator ˆ O = | m (cid:105) (cid:104) m | where | m (cid:105) isan eigenstate of the operator ˆ p with eigenvalue ¯ km and m ∈ Z . Its expectation is the occupation probabilityof the state | m (cid:105) and we will eventually find a diffusionequation for these occupations.We can expand this operator in the basis of the non-interacting Floquet modes [43, 44]. They are definedin terms of the noninteracting Floquet states whichare solutions of the noninteracting Schr¨odinger equation i ¯ k∂ t | ψ ( t ) (cid:105) = ˆ H ( t ) | ψ ( t ) (cid:105) which are periodic up to aphase | ψ j ( t ) (cid:105) = e − iµ j t | φ j ( t ) (cid:105) . The periodic part | φ j ( t ) (cid:105) are the Floquet modes. We have defined ˆ U ( t ) as thetime-evolution operator of ˆ H from time 0 to time t , andwe find ˆ U ( t ) | φ j (0) (cid:105) = e − iµ j t | φ j ( t ) (cid:105) . (28)So, we get the expressionˆ O ≡ (cid:88) i j O i j ( t ) | φ i ( t ) (cid:105) (cid:104) φ j ( t ) | , (29)where O i j ( t ) ≡ (cid:104) φ j ( t ) | ˆ O| φ i ( t ) (cid:105) is a periodic quantity.Now let us write its expectation at time t in the interac-tion representation as O ( t ) = Tr[ ˆ U † ( t ) ˆ O ˆ U ˆ ρ I ( t )]= (cid:88) i j O i j ( t )e i ( µ i − µ j ) t (cid:104) φ j (0) | ˆ ρ I ( t ) | φ i (0) (cid:105) (30)where we have used Eqs. (28) and (29). We apply thecoarse-graining average. Thanks to the separation oftimescales, it does not act on ˆ ρ I ( t ). It acts over manyperiods and, assuming that the µ j are incommensuratewith the driving frequency 2 π [38], we get O ( t ) = (cid:88) i O i i (cid:104) φ i (0) | ˆ ρ I ( t ) | φ i (0) (cid:105) . (31)Being O i j ( t ) periodic, O i i does not depend on the coarse-grained time. Using Eq. (27) and applying the cyclicproperty of trace, we getdd t O ( t ) = − σ ( t ) (cid:88) i O i i Tr (cid:110)(cid:2) [ | φ i (0) (cid:105) (cid:104) φ i (0) | , ˆ V I ( t )] , ˆ V I ( t ) (cid:3) ˆ ρ I ( t ) (cid:111) . (32) Expanding in the Floquet-mode basis ˆ V = (cid:80) k l V k l ( t ) | φ k ( t ) (cid:105) (cid:104) φ l ( t ) | [as in Eq. (29)] we obtaindd t O ( t ) = − σ ( t ) (cid:88) i k (cid:54) = i O i i | V k i | (cid:104) φ k (0) | ˆ ρ I ( t ) | φ k (0) (cid:105) (33)(we report the detailed derivation from Eq. (32) in Ap-pendix B). Dynamical localization implies Anderson lo-calization of the Floquet modes | φ k (0) (cid:105) in the momen-tum basis {| m (cid:105)} [18, 19, 24]. Moreover, ˆ V is localin momentum space because it only connects | m (cid:105) with | m + 1 (cid:105) . These two facts imply that | V k i | is local (orshort range) and is nonvanishing only if | j − k | is smallerthan some localization length λ . At this point we cansubstitute ˆ O = | m (cid:105) (cid:104) m | and define p m ( t ) ≡ | m (cid:105) (cid:104) m | ( t ) = (cid:104) m | ˆ ρ I ( t ) | m (cid:105) . Concerning this quantity, Anderson local-ization of Floquet states has another important conse-quence. If we make a coarse-graining average in momen-tum space, then (cid:104) φ k (0) | ˆ ρ I ( t ) | φ k (0) (cid:105) can be approximatedby (cid:104) m (cid:48) | ˆ ρ I ( t ) | m (cid:48) (cid:105) for some m (cid:48) , and we get ddt p m ( t ) (cid:39) − σ ( t ) (cid:88) m (cid:48) C m, m (cid:48) p m (cid:48) ( t ) , (34)where C m m (cid:48) is a coupling local in m and m (cid:48) . In orderto better specify the form of this local coupling, we con-sider that, thanks to the coarse graining in momentumspace, m can be assumed continuous. A constraint comesfrom the normalization condition (cid:80) m p m ( t ) = 1, whichis automatically enforced by writing the right-hand sideof Eq. (34) as a total derivative in mddt p m ( t ) = − σ ( t ) ∂ m J m ( { p m (cid:48) ( t ) } ) , (35)where J m ( { p m (cid:48) ( t ) } ) is linear in the p m (cid:48) ( t ). Due to the lo-cality of C m m (cid:48) , we can apply to J m ( { p m (cid:48) ( t ) } ) a gradientexpansion J m ( { p m (cid:48) } ) = C ( m ) p m − D ( m ) ∂ m p m +higher derivatives . (36)Using the assumption of average uniformity of the m -lattice and the symmetry m → − m of the model [seeEq. (15)] we can conclude that C ( m ) = 0 and that D ( m )does not depend on m .We then conclude that ddt p m ( t ) = σ ( t ) D ∂ m p m ( t ) (37)plus higher derivatives term, whose contribution goes to0 as the time goes on.This is a diffusion equation. If we initialize witha p m (0) symmetric for m → − m , the dynamics willpreserve this symmetry. This is our case because weinitialize with p m (0) = δ m . Therefore, the quantity m t ≡ (cid:80) m m p m ( t ) is the variance of the p m ( t ) distri-bution. Notice that e ( t ) is proportional to m t being e ( t ) = ¯ k m t . It is well known that the variance of adistribution obeying a diffusion equation like Eq. (37)obeys the equation dd t m t = 12 σ ( t ) D . (38)Notice that if σ ( t ) = const . , we find a diffusive behaviourof the energy, with e ( t ) ∝ m t linearly increasing in time,in agreement with the results known in literature for adriving undergoing a noise with properties invariant un-der time translations [35].In the case of our dynamics, Eq. (26) is valid and thenwe need then to estimate F ( t ) = Var [ F ( t )]. ApplyingEq. (20), we then havedd t m t ∝ Dm t . (39)Solving this simple differential equation we find e ( t ) = ¯ k m t ∝ D / t / . (40)This time dependence has been observed for long timesand (cid:15) ∈ [0 . ,
1] (see Fig. 3 and Ref. [24]).This theory gives rise to another prediction which wecan numerically check. In fact, combining the scaling ofthe energy with Eq. (20), we find σ ( t ) ∝ F ( t ) ∝ t − / . (41)Whenever there is power-law increase of the energy withexponent 2 /
3, this is exactly the time dependence of σ ( t )that we numerically observe (see Fig. 5). So, σ ( t ) de-creases with time and then the coarse-graining approxi-mation improves when larger and larger times are con-sidered. This fact explains why the curves in Fig. 3 tendtowards a power law with exponent 2 / σ ( t ) (cid:28)
1. Plugging Eq. (20) into thedefinition of σ ( t ), we have that σ ( t ) ∼ (cid:15) K m t . (42)Thus, for σ ( t ) (cid:28) m t (cid:29) (cid:15) K .This implies that for large K , there could be a very longtransient regime where the energy absorption could fol-low a different scaling in time. This could explain thelinear energy increase found in [24] for large K and (cid:15) ,which might be just a transient behaviour. For (cid:15) (cid:28)
1, incontrast, the chaotic behaviour of the system becomesapparent after a transient exponentially large in 1 /(cid:15) ,as a consequence of the Nekhoroshev theorem [6] (seeSec. VI). In this transient the system looks like integrableand shows dynamical localization. This fact could ex-plain the dynamical localization without energy increasethat we observe in Figs. 3 and 5 for (cid:15) = 0 .
01. In orderto see if at some point an energy increase with exponent2 / /(cid:15) ) = exp(10 ) (cid:39) . × , which is beyond ournumerical capabilities. σ ( t ) / ( g ε ) t K = 4 ε = 0.01 ε = 0.03 ε = 0.17 ε = 0.52 ε = 1.0t -1/3 σ ( t ) / ( g ε ) t K = 6 ε = 0.01 ε = 0.03 ε = 0.17 ε = 0.52 ε = 1.0t -1/3 σ ( t ) / ( g ε ) t K = 8 ε = 0.01 ε = 0.03 ε = 0.17 ε = 0.52 ε = 1.0t -1/3 FIG. 5. Bilogarithmic plot of σ ( t ) / ( g(cid:15) ) versus t ( t is coarsegrained with ∆ t = 50). Notice the power-law decrease con-sistent with the law σ ( t ) ∝ t − / when e ( t ) ∼ t / in Fig. 3.Other numerical parameters: ¯ k = 1 . M = 20153. V. LARGEST LYAPUNOV EXPONENT
Chaos was relevant in the last section in order to treat F ( t ) as an effective noise. We concentrate here in charac-terizing the chaotic properties of this model and we evalu-ate the largest Lyapunov exponent [39] λ ( T ) as a strobo-scopic average over T periods, and study its convergencewith T , using the method explained in [57]. We com-pare two trajectories with nearby initializations, one with β (0) = 1, β m (cid:54) =0 (0) = 0, another with β (cid:48) (0) = √ − δ , β (0) = δ , β m (cid:54) =0 , = 0. We choose δ = 10 − and we0study the rate of exponential divergence in each periodin the following way. At the end of the first period theinitial distance √ δ has become d . Before evolving overanother period, we leave ( . . . β m ( n ) . . . ) unchanged andmodify ( . . . β (cid:48) m ( n ) . . . ): We move it along the line joiningit with ( . . . β m ( n ) . . . ) so that the distance between thetwo trajectories become again δ . We iterate this proce-dure over T periods obtaining the sequence of distances d , d , . . . , d T . The largest Lyapunov exponent is evalu-ated as the average rate of exponential divergence λ ( T ) = 1 T T (cid:88) j =1 log (cid:18) d j √ δ (cid:19) . (43)This quantity reaches a limit λ over a finite T wheneverthe phase space is a bounded set.The phase space is bounded for any finite M . Foreach choice of M we take T so that we have reachedconvergence, and so we find that λ decays as a power lawwith 2 M +1 (Fig. 6). We also fit log( λ ) versus log(2 M +1)with a straight line and the slope of the fitting line is − . ± .
05. Therefore, in the limit M → ∞ the systemshould tend towards a regular behaviour with vanishing λ . This looks like a paradox for a thermalizing system(see Fig. 2), but we should not worry so much. We plot λ ( T ) versus T for different values of the truncation M inFig. 7. We see that λ ( T ) decreases as a power law with T , until it saturates to a plateau decreasing with 2 M + 1as a power law (see Fig. 6). For M → ∞ the power lawdecay would last forever: the phase space is not bounded,the dynamics wanders away towards m → ∞ and thesystem explores parts of the phase space with larger andlarger m and smaller and smaller exponential divergenceof the trajectories. This looks reasonable: the ratio ofthe perturbation inducing chaos and the unkicked partof the Hamiltonian is (cid:15)/m , so for large m there is lesschaos.For any finite T , λ ( T ) is nonvanishing because thedynamics is restricted to finite values of m and gives ameasure of the average Lyapunov exponent in that rangeof m . We can see that when you probe parts of thephase space with larger m , you get a smaller value of theLyapunov exponent. We notice that λ ( T ) relaxes to theplateau at the same time when the dynamics saturatesto the maximum possible attainable value of m and theenergy in Fig. 2 attains the T = ∞ value at finite M .An important information is that different timeregimes in Fig. 7 correspond to different ranges of m :The larger T , the larger m , the smaller the largest Lya-punov exponent λ ( T ). We can see this fact more ex-plicitly if we plot λ ( T ) versus m ( T ) ≡ (cid:112) e ( T ) / ¯ k , aswe do in Fig. 8. We find further confirmation of thisfact if we evaluate the largest Lyapunov exponent con-sidering different initializations. Beyond the one con-sidered up to now, we take another one where for onetrajectory β m (0) = 1 / √ M + 1 ∀ m and for the other β (cid:48) m (cid:54) =0 1 (0) = 1 / √ M + 1, β (cid:48) m (cid:54) =0 1 (0) = √ − δ/ √ M + 1, -8-7.5-7-6.5-6-5.5-5-4.5-4-3.5 3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 5.8 6 l og ( λ ) log(M) ε = 0.52Linear fit ε = 0.17Linear fit FIG. 6. Largest Lyapunov exponent λ versus 2 M + 1. Nu-merical parameters: K = 6, (cid:15) = 0 .
52, ¯ k = 1 . T ≥ .FIG. 7. Largest Lyapunov exponent λ ( T ) versus T for differ-ent values of M . Same numerical parameters as in Fig. 6 β (cid:48) m (cid:54) =0 1 (0) = √ δ/ √ M + 1. In the second initializa-tion scheme the energy expectation is equal to the T = ∞ value from the beginning and then the dynamics exploreslarge values of m from the beginning. We show an exam-ple of the situation in Fig. 9. In the upper panel we showthe evolution of λ ( T ). We see that with the second ini-tialization λ ( T ) is small from the beginning, consistentlywith the fact that the relevant values of m are larger andour assumption of a λ m decreasing with m . For both ini-tializations λ ( T ) reaches the same limit, and this occursat the time where the energy has reached the T = ∞ value for both the initializations (lower panel of Fig. 9). N (cid:29) A positive Lyapunov exponent has an important im-plication for the dynamics at N (cid:29) β m have fluctuations of order 1 / √ N has we can see in1 FIG. 8. Largest Lyapunov exponent λ ( T ) versus m ( T ) ≡ (cid:112) e ( T ) / ¯ k for different values of M . Same numerical param-eters as in Fig. 6 λ ( T ) β (0) = 1, β m ≠ = 0 β m (0) = 1 / (2M+1) e ( T ) T FIG. 9. Largest Lyapunov exponent λ ( T ) versus T for dif-ferent initializations (upper panel) and corresponding en-ergy evolution (lower panel). Numerical parameters: M =153 , K = 4 , (cid:15) = 0 . , ¯ k = 1 . Eq. (6). Because of chaotic dynamics, this initial uncer-tainty increases exponentially fast in time with a rategiven by the Lyapunov exponent. The larger m , thesmaller the Lyapunov exponent, with a dependence λ m resembling the one in Fig. 8. So, the Gross-Pitaevskijequations Eq. (7) are valid until the initial uncertaintybecomes of order 1. This occurs for a time t such thate λ m t / √ N ∼
1, that’s to say t ∼ λ m log N . (44)For m → ∞ our results suggest that λ m →
0, and so inthat limit the Gross-Pitaevskij equations are true for atime tending to infinity, even for finite N . VI. CONSERVED QUANTITIES AND THENEKHOROSHEV THEOREMA. Nekhoroshev theorem in a nutshell
Nekhoroshev theorem [6] applies when an integrablesystem is perturbed with a term breaking the integrabil-ity. A classical integrable system is such if it has as manyintegrals of motion I j as degrees of freedom and theseintegrals of motion are in involution, this gives to thedynamics peculiar regularity properties [53]. Perturb-ing the integrable system with an integrability breakingterm, these quantities are no more conserved and they de-pend on time: I j ( t ) deviate from their initial value. Letus call (cid:15) the strength of the perturbation. Nekhoroshevtheorem says that there are two positive real numbers a , b such that | I j ( t ) − I j (0) | ≤ (cid:15) b whenever t ≤ exp( C(cid:15) − a ) (45)for some positive constant C . This is a very important in-formation because it tells us that for (cid:15) (cid:28) I j ( t ) areapproximately conserved for a time exponentially largein the inverse perturbation, that’s why the Nekhoroshevtheorem is also called “classical prethermalization”. Weare going to show that it is valid also for our Gross-Pitaevskij equations. B. Conserved quantities at (cid:15) = 0
Let us start focusing on (cid:15) = 0. In this case, Eqs. (7) co-incide with the Schr¨odinger equation of the single kickedrotor in the momentum basis (see for instance [50]). Nowwe are going to show that this noninteracting model, ifprobed stroboscopically in time, shows infinitely manyconserved quantities local in m . They can be con-structed by means of the Floquet diagonalization ofEq. (7) at (cid:15) = 0. This procedure is equivalent to eval-uating the single-rotor Floquet modes [43, 44] in mo-mentum basis. To find them, one gets from Eq. (28)the eigenvalue equation ˆ U (1) | φ j (0) (cid:105) = e − iµ j t | φ j (0) (cid:105) and solves it, expanding ˆ U (1) and | φ j (0) (cid:105) in the mo-mentum eigenbasis. More specifically, we are focus-ing here on the time immediately after the kick, sowe consider the form ˆ U (1 + ) = e − iK cos ˆ θ e − i ˆ p / (2¯ k ) .The eigenvalues are of the form e − iµ j ; the eigenvec-tors are the Floquet modes and they have the form V j ≡ (cid:0) . . . U m − j U m j U m +1 j . . . (cid:1) T (in the notationof Sec. IV C we have | φ j (0) (cid:105) = (cid:80) m U m j | m (cid:105) ).From one side we see that if we prepare the system inthe condition β m (0) = U m j ∀ m , the evolution reducesto β m ( n ) = e − inµ j U m j being V j a Floquet mode of thesingle-rotor dynamics. From the other side, we see thatthe quantity O j ( n ) = ∞ (cid:88) m = −∞ U ∗ m j β m ( n ) (46)2with n integer evolves as O j ( n ) = e iµ j n O j (0), so that | O j | is conserved, whatever are the initial conditions inthe β m . We notice that the | O j | are local in m due to theAnderson localization of the Floquet states [18, 19, 24]Moreover, we can show that the | O j | are in involu-tion, that’s to say their Poisson bracket vanishes. Thiscan be easily seen by evaluating the Poisson bracket {| O j | , | O j (cid:48) | } and showing that it vanishes by using theelementary Poisson brackets Eq. (10) and the orthonor-mality condition (cid:80) m U ∗ m j U m j (cid:48) = δ j j (cid:48) .This is an important remark because it implies thatthe (cid:15) = 0 system is an integrable system in the classicalHamiltonian sense [53]. From one side the integrals ofmotion are as many as the degrees of freedom (imposinga truncation one sees that the j are as many as the m ).From the other side the Poisson brackets of the integralsof motion vanish so that they are in involution. So, ifwe apply a perturbation to it, the Nekhoroshev theoremshould be valid, as we are going to explicitly show now. C. Nekhoroshev estimate for (cid:15) > It is very interesting to study the time dependence of | O j ( t ) | (cid:15)> when (cid:15) > β m (0) = 1 / √ M + 1. Weconsider δ j, (cid:15) ( t ) = || O j ( t ) | (cid:15)> − | O j | (cid:15) =0 | . (47)We study the evolution in time of this quantity. In or-der to get rid of any influence of the initial state, weaverage over N diso initial conditions taken as β m (0) =e − iφ m / √ M + 1 with φ m ∈ [0 , π ] uniformly distributedrandom variables and call the average (cid:104) δ j, (cid:15) (cid:105) ( t ). We showsome examples of evolution of η j, (cid:15) ( t ) ≡ (cid:104) δ j, (cid:15) (cid:105) ( t ) (cid:104)| O j | (cid:15) =0 (cid:105) (48)in Fig. 10. We see a initial power-law increase followedby a saturation to a plateau. We can fit the power lawby means of a linear fit of the bilogarithmic plotlog η j, (cid:15) ( t ) = A j, (cid:15) log t + B j, (cid:15) . (49)We average both A j, (cid:15) and B j, (cid:15) over j and get A (cid:15) = 12 M + 1 M +1 (cid:88) j =1 A j, (cid:15) B (cid:15) = 12 M + 1 M +1 (cid:88) j =1 B j, (cid:15) . (50)We plot A (cid:15) and B (cid:15) versus (cid:15) in Fig. 11. We see twoimportant properties which will be important in the fol-lowing. The first one is that A (cid:15) is almost constant in (cid:15) and equals ∼ .
53. The second one, emphasised by the η ε , j t ε = 0.03, j = 1 ε = 0.03, j = 49 ε = 0.17, j = 1 ε = 0.17, j = 49 ε = 0.52, j = 1 ε = 0.52, j = 49 FIG. 10. η j, (cid:15) ( t ) versus t for different values of j and (cid:15) . Nu-merical parameters: K = 6, ¯ k = 1 . M = 153, N diso = 624. A ε - B ε ε C/ ε a FIG. 11. A (cid:15) and − B (cid:15) versus (cid:15) . The errorbars are evaluatedas the averages over j of the errorbars given by the fittingprocedure: they are much larger than the standard deviationsover j of respectively A j, (cid:15) and B j, (cid:15) . In the lower panel weplot also the result of the linear fit of the bilogarithmic plot.Numerical parameters: K = 6, ¯ k = 1 . M = 153, N diso =624. bilogarithmic plot, is that − B (cid:15) decays in (cid:15) as a powerlaw, − B (cid:15) ∼ C/(cid:15) a with a = − . ± . . . . . (51)where a has been numerically obtained through the lin-ear fit of the bilogarithmic plot. Being the errorbars quitesmall in relation to the values of A (cid:15) and B (cid:15) , we can focuson an average value of η defined aslog η (cid:15) ( t ) ≡ A (cid:15) log t + B (cid:15) . (52)3We give an estimate of the time when the deviation fromthe initial value gets significant enough. Consistentlywith the Nekhoroshev analysis, we estimate this time t ∗ as the time when the condition η (cid:15) ( t ∗ ) = (cid:15) b is verified. If b > t ∗ is still inthe regime of power-law increase, using Eq. (49) we getthe result t ∗ = (cid:15) b/A (cid:15) e − B (cid:15) /A (cid:15) . (53)We have noticed that A (cid:15) (cid:39) .
53. So, let us choose any b > .
6. Moreover, taking a C (cid:48) slightly larger than C inEq. (51), we have − B (cid:15) ≤ C (cid:48) /(cid:15) a , as we can easily see inthe lower panel of Fig. 11. In this way we have that η (cid:15) ( t ),the average deviation of the O j from the (cid:15) = 0 conservedvalues, is smaller or equal than (cid:15) b when t ∗ ≤ (cid:15) b/A (cid:15) e C (cid:48)(cid:48) /(cid:15) a ≤ e C (cid:48)(cid:48) /(cid:15) a (54)where we have used Eq. (51) and defined C (cid:48)(cid:48) = C (cid:48) / . VII. CONCLUSIONS
In conclusion we have studied the dynamics of theinfinite-range coupled quantum kicked rotors. By map-ping it over a model of interacting bosons we haveperformed exact diagonalization for quite large systemsizes and truncations. We have analyzed the averagelevel spacing ratio and we have found a generalized ten-dency towards ergodicity for increasing system sizes N ,in agreement with previous analytical demonstrations.Then we have moved to the thermodynamic limitwhere the model is described by a system of Gross-Pitaevskij equations which reduces to the Schr¨odingerequation of the non-interacting model when the magni-tude (cid:15) of the interaction term vanishes. For (cid:15) (cid:54) = 0, theseequations are equivalent to the dynamics of a single-rotornonlinear effective Hamiltonian. This system gives riseto a power-law increase of the energy in time with expo-nent γ ∼ / /
3. Moreover, we have predictedthat the nonlinear modulation of the kicking, squaredand coarse-grained in time, showed a power-law time de-pendence of the form ∼ t − / which we have numericallyverified. Remarkably, considering a the kick modulatedby a noise with properties invariant under time trans-lation, we could get the diffusion of the energy foundin [35]. We have shown that the Gross-Pitaveskij equations wehave found are equivalent to a previously existing mean-field analysis of this model and we have studied chaos inthem by means of the largest Lyapunov exponent. Thisquantity is a measure of chaos being an estimate of therate of exponential divergence of nearby trajectories. Wefind that it decreases towards zero as a power law whenlooking to portions of the phase space with larger andlarger momentum. From this fact follows that, for finite N , the time of validity of the Gross-Pitaevskij equationsdiverges with the momentum. Indeed, we have shownthat this time is proportional to log N and inversely pro-portional to the rate of exponential divergence of the tra-jectories.Later, we have considered the limit of (cid:15) = 0 and usedFloquet theorem to construct integrals of motion of thestroboscopic dynamics. These integrals of motion are asmany as the degrees of freedom and are in involution witheach other, so the system is integrable in a classical sense.Taking 0 < (cid:15) (cid:28) (cid:15) and consistent with the predictionsof the Nekhoroshev theorem.Future perspectives include the study of the modelwhere the infinite-range interactions are replaced by long-range power-law interactions. In this case the full per-mutation symmetry is broken and a strong qualitativechange of the Hilbert space structure occurs [58] whichmight turn the subdiffusion into diffusion. On the otherhand, one might expect that the 2 / ACKNOWLEDGMENTS
A. R. thanks D. Rossini and A. Tomadin for the ac-cess to the late GOLDRAKE cluster, where part of thenumerical analysis for this project was performed.4
Appendix A: Bosonic representation
We can now discuss the bosonic mapping of the notes.This mapping was invented by [36] for a similar infinite-range coupled model. Let us consider a system of N sites, and we explicitly work out the representation ofthe operator ˆ O = N (cid:88) l =1 e − i ˆ θ l (A1) (for the other operators the procedure is very similar).We consider the basis {| m (cid:105) l } of eigenstates of the oper-ator ˆ p l with eigenvalues ¯ km with m ∈ Z . In this basis,the operator ˆ O readsˆ O = N (cid:88) l =1 ∞ (cid:88) m = −∞ | m + 1 (cid:105) l l (cid:104) m | . (A2)Because the system is fully symmetric under permuta-tions, we can restrict to the states even under all thepossible N ! permutations. If we call ˆ P the sum of allthe permutation operators, we can take as basis of ourHilbert space the states | . . . n m − n m n m +1 . . . (cid:105) ≡ (cid:113) N ! (cid:81) ∞ m = −∞ n m ! ˆ P | . . . ( m − . . . m − (cid:124) (cid:123)(cid:122) (cid:125) n m − ( m . . . m ) (cid:124) (cid:123)(cid:122) (cid:125) n m ( m + 1 . . . m + 1) (cid:124) (cid:123)(cid:122) (cid:125) n m +1 . . . (cid:105) . (A3)There are N sites, so (cid:80) ∞ m = −∞ n m = N . The factor infront is for normalization, you divide by √ N ! becausethere are N ! possible permutations. And you divide alsoby all the √ n m !. Consider for instance m = 1. Fixing everything else, there are n ! ways of rearrange the siteswith m = 1 and this increases the norm by a factor n !.You divide by √ n ! and the norm is again 1.Consider for instance the application of the operatorˆ O [Eq. (A2)]. You findˆ O | . . . n m − n m n m +1 . . . (cid:105) = 1 (cid:113) N ! (cid:81) ∞ m = −∞ n m ! ˆ P N (cid:88) l =1 ∞ (cid:88) m (cid:48) = −∞ | m (cid:48) + 1 (cid:105) l l (cid:104) m (cid:48) | | . . . ( m − . . . m − (cid:124) (cid:123)(cid:122) (cid:125) n m − ( m . . . m ) (cid:124) (cid:123)(cid:122) (cid:125) n m ( m + 1 . . . m + 1) (cid:124) (cid:123)(cid:122) (cid:125) n m +1 . . . (cid:105) . (A4)Using permutation symmetry, it is quite immediate to see thatˆ O | . . . n m − n m n m +1 . . . (cid:105) = 1 (cid:113) N ! (cid:81) ∞ m = −∞ n m ! ˆ P N (cid:88) l =1 ∞ (cid:88) m (cid:48) = −∞ n m (cid:48) | . . . ( m (cid:48) . . . m (cid:48) ) (cid:124) (cid:123)(cid:122) (cid:125) n m (cid:48) − ( m (cid:48) + 1 . . . m (cid:48) + 1) (cid:124) (cid:123)(cid:122) (cid:125) n m (cid:48) +1 ( m (cid:48) + 2 . . . m (cid:48) + 2) (cid:124) (cid:123)(cid:122) (cid:125) n m (cid:48) +2 . . . (cid:105) . (A5)Slightly rearranging we findˆ O | . . . n m − n m n m +1 . . . (cid:105) = ˆ P N (cid:88) l =1 ∞ (cid:88) m (cid:48) = −∞ (cid:112) n m (cid:48) ( n m (cid:48) +1 + 1) | . . . ( m (cid:48) . . . m (cid:48) ) (cid:124) (cid:123)(cid:122) (cid:125) n m (cid:48) − ( m (cid:48) + 1 . . . m (cid:48) + 1) (cid:124) (cid:123)(cid:122) (cid:125) n m (cid:48) +1 ( m (cid:48) + 2 . . . m (cid:48) + 2) (cid:124) (cid:123)(cid:122) (cid:125) n m (cid:48) +2 . . . (cid:105) (cid:113) N ! (cid:81) m (cid:48) − m = −∞ n m !( n (cid:48) m − n m (cid:48) +1 + 1)! (cid:81) ∞ m = m (cid:48) +2 n m ! . (A6)From a formal point of view this is entirely equivalent to writeˆ O | . . . n m − n m n m +1 . . . (cid:105) = ∞ (cid:88) m (cid:48) = −∞ ˆ b † m (cid:48) +1 ˆ b m (cid:48) | . . . n m − n m n m +1 . . . (cid:105) (A7)where ˆ b m (cid:48) are bosonic operators and | . . . n m − n m n m +1 . . . (cid:105) is an eigenstate of the corresponding number operatorsˆ n m (cid:48) with eigenvalues n m (cid:48) . In a fully analogous way one can show that N (cid:88) l =1 e i ˆ θ l | . . . n m − n m n m +1 . . . (cid:105) = ∞ (cid:88) m (cid:48) = −∞ ˆ b † m (cid:48) ˆ b m (cid:48) +1 | . . . n m − n m n m +1 . . . (cid:105) N (cid:88) l =1 ˆ p l | . . . n m − n m n m +1 . . . (cid:105) = ¯ k ∞ (cid:88) m (cid:48) = −∞ m (cid:48) ˆ n m (cid:48) | . . . n m − n m n m +1 . . . (cid:105) . (A8)5In the subspace even under all the permutation symmetries, we can indeed write the two parts of the HamiltonianEq. (1) as 12 N (cid:88) l =1 ˆ p l = ¯ k ∞ (cid:88) m (cid:48) = −∞ m (cid:48) ˆ n m (cid:48) V (ˆ θ ) = ¯ kK (cid:104) ∞ (cid:88) m = −∞ (cid:16) ˆ b † m ˆ b m +1 + H . c . (cid:17) − (cid:15) N − ∞ (cid:88) m,m (cid:48) = −∞ (ˆ b † m +1 ˆ b m ˆ b † m (cid:48) ˆ b m (cid:48) +1 + H . c . ) (cid:105) , (A9)where V (ˆ θ ) is defined in Eq. (2). Appendix B: From Eq. (32) to Eq. (33)
We start considering the nested commutator inEq. (32). Expanding ˆ V in the Floquet-mode basis asˆ V = (cid:80) k l V k l ( t ) | φ k ( t ) (cid:105) (cid:104) φ l ( t ) | with V k l ( t ) time-periodicwith period 1 we get (cid:2) [ | φ i (0) (cid:105) (cid:104) φ i (0) | , ˆ V I ( t ) (cid:3) , ˆ V I ( t )] = (cid:88) k l k (cid:48) l (cid:48) V k l ( t ) V k (cid:48) l (cid:48) ( t ) (cid:2) [ | φ i (0) (cid:105) (cid:104) φ i (0) | , ˆ U † ( t ) | φ k ( t ) (cid:105) (cid:104) φ l ( t ) | ˆ U ( t )] , ˆ U † ( t ) | φ k (cid:48) ( t ) (cid:105) (cid:104) φ l (cid:48) ( t ) | ˆ U ( t ) (cid:3) = (cid:88) k l k (cid:48) l (cid:48) V k l ( t ) V k (cid:48) l (cid:48) ( t ) ˆ U † ( t ) (cid:2) [ | φ i ( t ) (cid:105) (cid:104) φ i ( t ) | , | φ k ( t ) (cid:105) (cid:104) φ l ( t ) | ] , | φ k (cid:48) ( t ) (cid:105) (cid:104) φ l (cid:48) ( t ) | (cid:3) ˆ U ( t ) , (B1)where we have applied Eq. (28). Exploiting the orthonormality of the Floquet-mode basis we can evaluate in a lengthybut straightforward way the double commutator and get (cid:2) [ | φ i (0) (cid:105) (cid:104) φ i (0) | , ˆ V I ( t )] , ˆ V I ( t ) (cid:3) = (cid:88) k l ˆ U † ( t ) (cid:110) V i k ( t ) V k l ( t ) | φ i ( t ) (cid:105) (cid:104) φ l ( t ) | − V k i ( t ) V i l ( t ) | φ k ( t ) (cid:105) (cid:104) φ l ( t ) | + V l i ( t ) V k l ( t ) | φ k ( t ) (cid:105) (cid:104) φ i ( t ) | − V k i ( t ) V i l ( t ) | φ k ( t ) (cid:105) (cid:104) φ l ( t ) | (cid:111) ˆ U ( t )= (cid:88) k l (cid:110) V i k ( t ) V k l ( t )e i ( µ i − µ k ) t | φ i (0) (cid:105) (cid:104) φ l (0) | − V k i ( t ) V i l ( t )e i ( µ k − µ l ) t | φ k (0) (cid:105) (cid:104) φ l (0) | + V l i ( t ) V k l ( t )e i ( µ k − µ i ) t | φ k (0) (cid:105) (cid:104) φ i (0) | − V k i ( t ) V i l ( t )e i ( µ k − µ l ) t | φ k (0) (cid:105) (cid:104) φ l (0) | (cid:111) . (B2)At this point we perform the coarse-graining average. Weaverage over ∆ t which lasts many periods and we as-sume that the µ k are incommensurate with the drivingfrequency 2 π [38]. In this way we obtain (cid:2) [ | φ i (0) (cid:105) (cid:104) φ i (0) | , ˆ V I ( t )] , ˆ V I ( t ) (cid:3) = 2 (cid:88) k (cid:54) = i | V k i | | φ k (0) (cid:105) (cid:104) φ k (0) | . (B3) Here | V k i | does not depend on the coarse-grained timebecause | V k i ( t ) | is periodic of period 1 and the averagingtime ∆ t spans many periods. Substituting this expres-sion into Eq. (32) we get Eq. (33). [1] B. V. Chirikov, Physics Reports , 263 (1979).[2] A. Lichtenberg and M. Lieberman, Regular and ChaoticMotion (Springer, 1992).[3] M. V. Berry, in
Topics in Nonlinear Mechanics , Vol. 46,edited by S. Jorna (Am.Inst.Ph., 1978) pp. 16–120. [4] K. Kaneko and T. Konishi, Phys. Rev. A , 6130 (1989).[5] T. Konishi and K. Kaneko, Journal of Physics A: Math-ematical and General , L715 (1990).[6] N. N. Nekhoroshev, Functional Analysis and Its Appli-cations , 338 (1971). [7] M. Ueda, Nat Rev Phys , 669–681 (2020).[8] M. V. Berry, Journal of Physics A: Mathematical andGeneral , 2083 (1977).[9] P. Pechukas, Phys. Rev. Lett. , 943 (1983).[10] M. Srednicki, Phys. Rev. E , 888 (1994).[11] M. Feingold and A. Peres, Phys. Rev. A , 591 (1986).[12] T. Prosen, Annals of Physics , 115.[13] O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev.Lett. , 1 (1984).[14] B. Eckhardt and J. Main, Phys. Rev. Lett. , 2300(1995).[15] B. Chirikov, F. Izrailev, and D. Shepelyansky, PhysicaD , 77 (1988).[16] B. V. Chirikov, in Chaos and quantum mechanics , LesHouches Lecture Series, Vol. 52, edited by M.-J.Giannoni,A. Voros, and J. Zinn-Justin (Elsevier Sci. Publ., Ams-terdam, 1991) p. 443–545.[17] G. Casati, B. V. Chirikov, J. Ford, and I. F. M., in
Lecture Notes in Physics , Vol. 93 (Springer, 1979) p. 334.[18] S. Fishman, D. R. Grempel, and R. E. Prange, Phys.Rev. Lett. , 509 (1982).[19] D. R. Grempel, R. E. Prange, and S. Fishman, Phys.Rev. A , 1639 (1984).[20] D. Delande, “Kicked rotor and anderson lo-calization (lectures 1 and 2),” (2013),http://boulderschool.yale.edu/sites/default/files/files/Delande-kicked rotor lectures 1 and 2.pdf.[21] F. L. Moore, J. C. Robinson, C. Bharucha, P. E.Williams, and M. G. Raizen, Phys. Rev. Lett. , 2974(1994).[22] F. L. Moore, J. C. Robinson, C. F. Bharucha, B. Sun-daram, and M. G. Raizen, Phys. Rev. Lett. , 4598(1995).[23] P. Qin, A. Andreanov, H. C. Park, and S. Flach, Scien-tific Reports , 41139 (2017).[24] S. Notarnicola, F. Iemini, D. Rossini, R. Fazio, A. Silva,and A. Russomanno, Phys. Rev. E , 022202 (2018).[25] E. B. Rozenbaum and V. Galitski, Phys. Rev. B ,064303 (2017).[26] C. Rylands, E. Rozenbaum, V. Galitski, and R. Konik,Phys. Rev. Lett. , 155302 (2020).[27] M. Fava, R. Fazio, and A. Russomanno, Phys. Rev. B , 064302 (2020).[28] S. Notarnicola, A. Silva, R. Fazio, and A. Russomanno,Jour. Stat. Mech. , 024008 (2020).[29] G. Gligori´c, J. D. Bodyfelt, and S. Flach, EPL , 30004(2011).[30] L. Ermann and D. L. Shepelyansky, Journal of PhysicsA , 335101 (2014).[31] A. S. Pikovsky and D. L. Shepelyansky, Phys. Rev. Lett. , 094101 (2008).[32] D. L. Shepelyansky, Phys. Rev. Lett. , 1787 (1993).[33] B. P. Nguyen and K. Kim, Eur. Phys. Jour. B , 79(2011).[34] N. Cherroret, B. Vermersch, J. C. Garreau, and D. De-lande, Phys. Rev. Lett. , 170603 (2014).[35] B. G. Klappauf, W. H. Oskay, D. A. Steck, and M. G.Raizen, Physical review letters , 1203 (1998).[36] F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva,R. Fazio, and F. Iemini, Phys. Rev. B , 104303 (2019). [37] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Atom-photon interactions: basic processes and applica-tions , Wiley-Interscience publication (J. Wiley, 1992).[38] A. Russomanno, S. Pugnetti, V. Brosco, and R. Fazio,Phys. Rev. B , 214508 (2011).[39] E. Ott, Chaos in dynamical systems ( nd Ed.) (Cam-bridge University Press, 2002).[40] P. Matus and K. Sacha, Phys. Rev. A (2019),10.1103/physreva.99.033626.[41] B. Chirikov and V. Vecheslavov, Journal of StatisticalPhysics , 243 (1993).[42] This symmetry comes from the symmetry ˆ θ j → − ˆ θ j ∀ j of Eq. (1).[43] J. H. Shirley, Phys. Rev. , B979 (1965).[44] H. Sambe, Phys. Rev. A , 2203 (1973).[45] A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010).[46] A. Russomanno, A. Silva, and G. E. Santoro, Phys. Rev.Lett , 257201 (2012), 1204.5084.[47] M. V. Berry, in Les Houches, Session XXXVI, 1981 —Chaotic Behaviour of Deterministic Systems , edited byR. S. G. Ioos, R. H. G. Helleman (North-Holland Pub-lishing Company, 1983) pp. 174–271.[48] A. Lazarides, A. Das, and R. Moessner, Phys. Rev. E , 012110 (2014).[49] L. D’Alessio and M. Rigol, Phys. Rev. X , 041048(2014).[50] F. Haake, Quantum Signatures of Chaos ( nd ed) (Springer, 2001).[51] B. Eynard, T. Kimura, and S. Ribault, (2018),arXiv:1510.04430 [math-ph].[52] M. V. Berry and M. Tabor, Proc. Roy. Soc. A349 , 101(1976).[53] V. I. Arnol’d,
Mathematical Methods of Classical Me-chanics (Springer, 1989).[54] B. Sciolla and G. Biroli, J. Stat. Mech.: Theor. and Ex-per. , P11003 (2011).[55] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery, Numerical Recipes in Fortran 90. The Art ofParallel Scientific Computing , 2nd ed., Vol. 2 (CambridgeUniversity Press, Cambridge, 1996).[56] The evolution between two kicks simply amounts to themultiplication of the ˆ β m by a phase, β m (( n + 1) − ) = β m ( n + )e − i ¯ k m / . The evolution across the n -th kickamounts to numerically solve the system of nonlinear dif-ferential equations i dd t β m ( t ) = K (cid:20) β m +1 ( t ) + β m − ( t ) − (cid:15) (cid:16) (cid:88) m (cid:48) β ∗ m (cid:48) ( t ) β m (cid:48) +1 ( t ) (cid:17) β m − ( t ) − (cid:15) (cid:16) (cid:88) m (cid:48) β m (cid:48) ( t ) β ∗ m (cid:48) +1 ( t ) (cid:17) β m +1 ( t ) (cid:21) (B4)with t ∈ [0 ,
1] and { β m ( n − ) } as initial condition.[57] G. Benettin, L. Galgani, and J.-M. Strelcyn, Phys. Rev.A , 2338 (1976).[58] A. Russomanno, M. Fava, and M. Heyl, (2020),arXiv:2012.06505 [cond-mat.stat-mech].[59] B. ˇZunkoviˇc, M. Heyl, M. Knap, and A. Silva, Phys.Rev. Lett.120