Long-time behavior of periodically driven isolated interacting lattice systems

We study the dynamics of isolated interacting spin chains that are periodically driven by sudden quenches. Using full exact diagonalization of finite chains, we show that these systems exhibit three distinct regimes. For short driving periods, the Floquet Hamiltonian is well approximated by the time-averaged Hamiltonian, while for long periods the evolution operator exhibits properties of random matrices of a Circular Ensemble (CE). In-between, there is a crossover regime. Based on a finite-size scaling analysis and analytic arguments we argue that, for thermodynamically large systems and non-vanishing driving periods, the evolution operator always exhibits properties of CE random matrices. Consequently, the Floquet Hamiltonian is nonlocal and has multi-body interactions; and the driving leads to the equivalent of an infinite temperature state at long times. These results are connected to the breakdown of the Magnus expansion and are expected to hold beyond the specific lattice model considered.

Introduction. Periodically driven systems have a long history, one paradigmatic example is the Kicked-Rotor model of a particle moving on a ring subjected to timeperiodic "kicks" [1]. The behavior of such systems is very rich, e.g., they can display interesting integrabilityto-chaos transitions and dynamical Anderson localization [2,3] and counter intuitive effects such as dynamical stabilization [4,5] both in classical and quantum mechanics. Moreover, it has been recently shown that periodic perturbation can be used as a flexible experimental knob [6][7][8] to realize synthetic matter, i.e., matter with specific engineered properties. Two outstanding examples in this direction are the opening of a gap in graphene by using light or carefully selected lattice phonons [9][10][11] and the realization of the so called Floquet topological insulator in a material which, in absence of the driving, is topologically trivial [12,13].
With a few exceptions [14][15][16][17][18], recent studies mostly focus on driving simple (often single-particle) Hamiltonians. Furthermore, theoretical analyses often rely on approximations that allow one to recast the effect of the drive in, e.g., an effective hopping [19,20]. The validity of some of those approximations, such as the use of average Hamiltonians, has been a topic of interest for many years because of its relevance to NMR experiments [21,22]. Our goal in this Letter is to understand the long-time behavior of generic isolated quantum systems under a periodic drive. We pay special attention to the role of finite size effects, which may be relevant to experiments with ultracold gases. Combining theoretical arguments and numerical results we argue that finite systems exhibit three different regimes, while only one survives in the thermodynamic limit. In the latter, the system approaches infinite temperature at infinite time and the Floquet Hamiltonian becomes unphysical, i.e. it cannot be an extensive operator with few-body interactions.
Theoretical Analysis. The Floquet theorem states that the dynamics of periodically driven systems can be described by the time-independent Floquet Hamiltonian, H F , defined by: where T is the period of the driving,Û cycle is the exact evolution operator over a full cycle and |φ n and e −iθn are its eigenstates and eigenvalues. From Eq. (1) if follows that H F ≡ n |φ n ε n φ n | where the Floquet energies, ε n , can be shifted by integers of 2π /T without violating (1). We stress thatĤ F does not need to be a physical Hamiltonian. For example,Ĥ F may contain many-body interactions and may not be extensive. Moreover, in general, it is not possible to evaluateĤ F explicitly and one has to rely on approximations such as the Magnus expansion [23]. The latter is a perturbative expansion in T , which is guaranteed to converge for finite-dimensional Hamiltonians and sufficiently short periods [24].
In the limit T ∆ε n,m (τ )/ 1, one can expand the exponential in Eq. (5), and, to first order, the dependence on the period T disappears where we have used the fact that |φ n (τ ) is an orthonormal basis. The equation above is equivalent tô H F = 1 0 dτĤ(τ ) =Ĥ ave . In this limit, ε n are extensive, i.e. proportional to the volume of the system, V . Therefore, ∆ε n,m (τ ) between the lowest and highest eigenstates ofĤ F is ∝ V and, in order for the expansion of the exponential ∀ n, m, τ to be meaningful, one needs T const/V . In this limit, the phases of the evolution operator are proportional to the eigenenergies ofĤ ave . We note that the condition T const/V is only a sufficient condition forĤ F =Ĥ ave . For example, ifĤ(τ ) commutes at different times, [H(τ ), H(τ )] = 0, then it trivially follows thatĤ F =Ĥ ave independently of T .
Remarkably, evaluating Eq. (7) at τ = 1 allows one to connect the structure of the eigenstates ofĤ F to the statistics of the phases of the evolution operator. When the eigenstatesĤ F are chaotic, the RHS of (7) would be generically non-zero for any pair of states (this follows from recalling thatĤ(τ = 1) is a sum of local few-body operators) and the phases ofÛ cycle are expected to display repulsion [see LHS of Eq. (7)]. On the other hand, ifĤ F is noninteracting (or integrable) then the RHS of Eq. (7) can be zero for a large fraction of pairs of states and there needs not be repulsion in the phases. In this case, the spectrum ofÛ cycle could exhibit highly degenerate phases separate by finite gaps.
We should stress that, for T > T c = 2π/W , repulsion in the phases ofÛ cycle is incompatible withĤ F =Ĥ ave , where W = ave max − ave min ∝ V is the width of the spectrum of the average Hamiltonian. This because, for H F =Ĥ ave , the phases would be θ n = T ave n / and would not exhibit repulsion. The same reasoning leads to the conclusion that, for T > T c , repulsion in the phases is incompatible withĤ F being an extensive operator with few-body interactions. This type of operators do not have correlations between far-away energies levels and, therefore, can not produce repulsion for the phases. Actually, repulsion of the phases ofÛ cycle for T > T c hints [see Eq. (7)] thatÛ cycle should exhibit properties of random matrices belonging to the Circular Ensembles (CEs) [26][27][28][29]. CEs are the equivalent of Gaussian Ensembles (GEs) for unitary matrices. In CEs the phases display repulsion and the eigenstates are essentially random vectors [so that both sides of Eq. (7) are generically different from zero].
Numerical Results. We now turn to numerical simulations to understand what happens when one drives a realistic many-body system. Specifically, we use full exact diagonalization to study a spin-1/2 chain with periodic boundary conditions and Hamiltonian where σ α j , α = x, y, z, is the Pauli matrix for the j th spin, f (t) is a periodic function with zero average, f (t) = 1 for N < t/T < N + 1 2 and f (t) = −1 for N + 1 2 < t/T < N +1. The Hamiltonian (8) is relevant to nuclear magnetic resonance experiments in solids [30]. The next-nearest neighbor coupling J can be considered to be a phenomenological parameter that accounts for interactions that break the integrability of the spin-1/2 XXZ chain. We restrict our calculations to the sector with total momentum k = 0, magnetization m z = 1/3, and parity p = +1. For the largest system size considered (L = 24 spins) that sector contains D = 15581 states. The Hamiltonian parameters are J = 1, δJ = 0.2, and J = 0.8. They have been chosen to ensure that the average Hamiltonian is non-integrable and exhibits eigenstate thermalization [25,[31][32][33]. Given the time dependence chosen, we obtain the exact evolution operator over one cyclê We are first interested in understanding the properties of the phases θ n (obtained diagonalizingÛ cycle ) when compared to those obtained ifĤ F =Ĥ ave , i.e, θ ave n = T ave n / , as a function of the driving period T . In recent studies, the ratio of two consecutive energy gaps [34,35] has been proven to be a good quantum chaos indicator. We extend this idea to our setup [25] and measure the statistics of r = min(δn,δn+1) max(δn,δn+1) , where we have defined The values of r (symbols) and rave (dashed lines) are compared to the COE, GOE and POI predictions. The COE and GOE predictions are expected to coincide in the thermodynamic limit [25]. Vertical dotted lines depict the periods T1 = 2π /W , T2 = π /σ, and T3 = 2π /σ, where W is the width of the full spectrum and σ is the width of the Gaussian that best fits the density of states for L = 24 (W and σ are ∝ L [25]). (Inset) The period Tm at which r is minimum decreases with increasing system size.
the phase gap δ n ≡ θ n+1 −θ n . Values of r close to zero indicate the lack of phase repulsion. In Fig. 1, we show the average value of r for different values of T for the exact phases θ n (indicated by r ) and the phases θ ave n (indicated by r ave ) and compare them with the predictions of the Circular Orthogonal Ensemble (COE), the Gaussian Orthogonal Ensemble (GOE), and Poisson statistics (POI). The COE and GOE predictions in Fig. 1 have been computed for matrix of finite size and are expected to coincide in the thermodynamic limit [25].
In Fig. 1, we observe four regions in which r and r ave exhibit distinctive behavior. Note the vertical dotted lines that mark T 1 , T 2 , and T 3 (see Fig. 1 caption). i) For T < T 1 , r is independent of T , identical to r ave , and close to the GOE prediction. ii) For T 1 < T < T 2 , r and r ave are very close to each other and decrease with increasing T . iii) For T 2 < T < T 3 , r increases towards the COE prediction while r ave decreases towards the POI prediction. iv) For T > T 3 , r is again independent of T and close to the COE prediction, while r ave is very close and continues approaching the POI prediction. In the inset, we track the period T m at which r has its minimum. t can be seen to decrease with increasing system size.
Further insights on our driven system can be gained by studying the properties of the eigenvectors ofÛ cycle . Specifically, we would like to understand how they become different from the eigenvectors ofĤ ave as one increases T . We write the former in terms of the latter, |φ n = and D is the dimension of the Hilbert space, and compute the (Shannon) information entropy [33] S which measures the number of states |m ave that contribute to each |φ n . WhenĤ F =Ĥ ave , S n = 0. Aŝ H F deviates fromĤ ave , S n grows and it is expected to saturate at the COE prediction S n ≈ ln(0.48D). In Fig. 2, we show the average entropy S = D n=1 S n /D normalized by the COE prediction. We note that S grows monotonically with increasing the period and that the curves shift to the left as the system size increases. The latter is also seen in the inset, where we plot the period T inf at which each curve has an inflection point vs L. The variance (S n − S ) 2 1/2 is smaller than the size of the symbols used in Fig. 2. This makes apparent that, in average, eigenstates |φ n delocalize the same way in the basis {|m ave }. Contrary to r (in Fig. 1), S (in Fig. 2) is still changing for T ≥ T 3 . However, the COE prediction is reached for longer periods. This is not surprising since indicators based on the eigenvectors are known to suffer from stronger finite-size effects than level statistics indicators [33].
The fact thatÛ cycle shares properties of a COE has important consequences in the energy absorption of the system. To show it, we compute the expectation value ofĤ ave at time t = N T for N → ∞. This can be done using the so-called diagonal ensemble [32]  In the COE all eigenstates are essentially random vectors and the expectation value of few-body observables become n-independent [25]. This occurs in our system for φ n |Ĥ ave |φ n when T > T 3 , as shown in the inset in Fig. 3 for T = 6 and L = 24. Hence, for T > T 3 , the evaluation of Eq. (10) gives a result that corresponds to infinite temperature independently of |ψ 0 .
In Fig. 3, we show the energy absorbed at infinite time Q vs T . We have normalized Q so that −1 corresponds to no-energy absorption and 0 corresponds to the final energy being equal to that at infinite temperature, i.e., Q = Ĥ ave (t→∞)− Ĥ ave β=0 Ĥ ave β=0 − Ĥ ave (t=0) . The initial state is taken to be in thermal equilibrium (with respect toĤ ave ) with β = 0.5. Namely, |ψ 0 ψ 0 | = Z −1 m e −βE ave m |m ave m ave |, where Z is the partition function. Figure 3 shows that, with increasing T , Q approaches zero as expected. In addition, the value of T at which Q saturates close to zero decreases with increasing system size.
Discussion. Our numerical results support the initial conclusion thatĤ F =Ĥ ave for vanishingly small values of T (T T 2 in our case), whileÛ cycle shares properties of matrices of the COE for larger values of T (T T 3 in our case). In finite systems, there is a cross-over regime (for intermediate periods T ) in which some indicators (repulsion between consecutive levels) have already reached the COE predictions while others (the information entropy of eigenstates) have not (see Figs. 1 and 2). In this regime, a drive may lead to states in which the average energy in a cycle fluctuates between consecutive cycles or becomes stationary but not equal to the infinite temperature expectation value. This could be seen in current ex-periments with periodically driven cold-atomic systems.
The scaling of our results with system size (see Figs. 1-3) suggests that, in the thermodynamic limit, any nonvanishing driving period T leads to an evolution opera-torÛ cycle that shares properties with random matrices of the COE. As a result, the system reaches the infinite temperature limit for infinite driving times. When this occurs, the Floquet HamiltonianĤ F cannot be an extensive Hamiltonian with few-body interactions because the latter would not produce phase repulsion for nonvanishing values of T . Therefore finite interacting systems driven at intermediate periods offer a natural platform to investigate unique phenomena that may occur as a result of long-range and/or many-body interactions inĤ F .
While we expect that our conclusions apply to generic interacting many-body systems, there are specific classes of models that could lead toĤ F being a well-behaved local Hamiltonian over a finite range of driving periods. Examples are, (i) driven Hamiltonians that commute at different times, Ĥ (τ ),Ĥ(τ ) = 0, which lead toĤ F = H ave independently of T (in this case both the LHS and RHS of Eq. (7) are identically zero), ii) fictitious timedependence that can be removed by a proper choice of the reference frame (see, e.g., [11]), iii) special situations in which multiple commutators ofĤ(τ ) generate a finite algebra [14,37].
In conclusion we have shown that isolated finite interacting systems subjected to a periodic driving exhibit three possible physical regimes. For short driving period T ,Ĥ F converges to the time-averaged Hamiltonian H ave , while for large periods the evolution operatorÛ cycle shares properties of matrices of the COE of random matrix theory. For intermediate periods, we find that there is a cross-over regime in which some indicators have already approached the COE predictions while others have not. In the latter regime finite systems may exhibit unexpected behavior, e.g., they may not heat up to infinite temperature and may reach interesting stationary states that could be observed in current experiments with coldatomic systems. We have provided evidence that, in the thermodynamic limit, the only regime that occurs for nonvanishing driving periods is the one in whichÛ cycle exhibits properties of COE random matrices. This results in the system heating up to infinite temperature at infinite times (see Fig. 3). Furthermore,Ĥ F cannot be an extensive Hamiltonian with few-body interaction in that regime. An exciting question that we plan to explore next is the effect of a thermal bath on the dynamics of interacting periodically driven systems [22,[38][39][40].

in the SE equation (3) one obtains
Using the mathematical identity to compute the derivative of an exponential operator (note that this expression is required becauseĤ eff (τ ) at different times may not commute) results in Multiplying both sides in the equation above by exp −TĤ eff (τ ) i from the right one gets We now close this equation between the exact (time-dependent) eigenstates of H eff (τ ) [Ĥ eff (τ )|φ n (τ ) = ε n (τ )|φ n (τ ) ] and perform the integration over α to obtain Eq. (5). To obtain (6) it suffices to "bracket" Eq. (11) with the same eigenstate and use the Hellmann-Feynman theorem φ n (τ )|∂ τĤeff (τ )|φ n (τ ) = ∂ τ ε n (τ ).

DENSITY OF STATES OFĤave
In order to compute the density of states ofĤ ave , we take advantage of the symmetries of the model [Eq. (8) in the main text] and restrict ourselves to a single sector with total momentum k = 0, magnetization m z = 1/3, and parity p = +1. For the four system sizes considered in our study L = 15, 18, 21, and 24, the number of states in that sector are D = 111, 561, 2829, and 15581, respectively.
In Fig. 4, we show the density of states ofĤ ave for the three largest systems considered. The lines are fits that make apparent that the density of states are well described by Gaussians. In the inset, we plot the total width of the spectrum, W = max ave − min ave , and the width of the fitted Gaussian √ σ 2 vs the number of spins L. Their scaling is consistent with the expected linear dependence. Note that W and σ were used to define the periods T 1 = 2π /W , T 2 = π /σ, and T 3 = 2π /σ in the main text, which are all expected to collapse to zero in the thermodynamic limit as they scale ∝ L −1 .

PHASE REPULSION IN CIRCULAR ENSEMBLES
In what follows we are interested in understanding the functional form of the probability distribution P (r), where r ≡ min(δ n , δ n+1 ) max(δ n , δ n+1 ) ∈ (0, 1), in different ensembles. In Gaussian ensembles (GEs), δ n ≡ n+1 − n is the n-th energy gap while in Circular ensembles (CEs) δ n ≡ θ n+1 − θ n is the n-th phase gap. P (r) has been shown to be a sensitive and practical probe [34,35] (especially because no-unfolding procedure is necessary) of the integrability-to-chaos transition [36]. For the GEs, the distribution P (r) has been recently calculated [35] starting from the joint distribution of the eigenvalues of a 3 × 3 random matrix. We now repeat the calculation for a 3 × 3 random matrix in the CEs.
The joint distribution for the phases in the Circular Orthogonal Ensemble (COE), the Circular Unitary Ensemble (CUE), and Circular Symplectic Ensemble (CSE) are [26][27][28][29]: and β = 1, 2, 4 for COE, CUE and CSE, respectively. For a 3 × 3 matrix the expression simplifies to: where θ i ∈ (0, 2π) [this follows from the normalization 2π 0 dθ 1 dθ 2 dθ 3 ρ(θ 1 , θ 2 , θ 3 ) = 1]. We now choose one of the possible six ordering of the angles, θ 1 ≤ θ 2 ≤ θ 3 , define the variable of interest in one of the possible three ways δ r − θ2−θ1 θ3−θ2 , and perform the change of variables x = θ 2 − θ 1 , θ 2 , y = θ 3 − θ 2 to obtain: where the factor of 2π comes from the integration 2π 0 dθ 2 , the factor of 6 comes from the choice of the ordering, and the factor of 3 comes from the choice of the observable. We also used the identity sin 2π−x−y 2 = sin x+y 2 . The limits of integration in Eq. (13) come from the fact that 0 ≤ x + y ≤ 2π (see Fig. 5). In Eq. (13),r = δ 1 /δ 2 ∈ (0, ∞) and is different from r defined in Eq. (12). Luckily, the relation between P (r) and P (r) is a simple one, P (r) = 2P (r)Θ(1−r) [35]. The integration in (13) can be performed exactly, we have used Mathematica. For example, the distribution for the COE is: For the sake of brevity the expressions for the CUE and CSE are not reported. In Fig. (6), we plot the distributions P COE (r), P CUE (r), P CSE (r) and compare them with the corresponding distribution for the GEs [P GOE (r), P GUE (r), P GSE (r)] and the distribution for a Poisson sequence of random numbers, P POI (r). The distributions for the CE and the corresponding GE are very close for 3 × 3 matrices and are expected to coincide in the thermodynamic limit. In fact, the phases in the Circular Ensembles (CEs) are expected to have the same local fluctuation properties as the eigenvalues of the corresponding Gaussian Ensembles (GEs) [29]. In Table I, we report the average value of r for all the distributions considered. Finally, we report numerical results for the model studied in the main text [see Eq. (8)] when L = 24. In Fig. 7, we compare the full distribution P (r) for both the exact phases (obtained from the exact diagonalization ofÛ cycle ) and the phases one obtains ifĤ F =Ĥ ave , θ ave n = T ave n / , at the periods T 1 , T 2 , T 3 (see Fig. 1 in the main text) with the GOE, COE and POI predictions. These plots show that at T 1 and T 3 the numerical results for the exact phases  Fig. 1 in the main text. We compare the exact phases (obtained from the exact diagonalization ofÛ cycle ) (indicated by red circles) and the phases one obtains if HF =Ĥave, θn = T ave n / (indicated by squares), with the GOE, COE and POI predictions (lines).
-π -π/2 0 π/2 π θ n are virtually indistinguishable from the predictions for 3 × 3 matrices in the GOE and the COE. On the other hand the phases θ ave n are clearly described by Poisson statistics at T 3 .

EXPECTATION VALUES OFĤave IN THE FLOQUET EIGENSTATES
We now investigate the dependence of the expectation values φ n |Ĥ ave |φ n on the phases θ n for L = 24 and different driving periods T . |φ n and θ n are the exact eigenstates and eigenvalues of the evolution operator over a cycle, i.e. U cycle = n |φ n e −iθn φ n |. For short periods, the Floquet eigenstates, |φ n , are closely related to the eigenstates of H ave and the simple relation θ n = T φ n |Ĥ ave |φ n holds. We note that the phases θ n might or might not span the entire "Floquet Brillouin zone", i.e. (−π, π), depending on the value of the period T (see left panel in Fig. 8). For large periodsÛ cycle shares the properties of the matrix of the Circular Ensemble and the Floquet eigenstates are essentially random vectors. In this regime the expectation values φ n |Ĥ ave |φ n approach an n-independent value and they are not related to the values of the Floquet phases θ n (see right panel in Fig. 8). For even larger periods the eigenstate to eigenstate variation of φ n |Ĥ ave |φ n decreases (see inset in Fig. 3 in the main text). For intermediate periods, there is a crossover region in which the relation between θ n and φ n |Ĥ ave |φ n is progressively lost by increasing the period T (see central panel in Fig. 8).