Aperiodically driven integrable systems and their emergent steady states

Does a closed quantum many-body system that is continually driven with a time-dependent Hamiltonian finally reach a steady state? This question has only recently been answered for driving protocols that are periodic in time, where the long time behavior of the local properties synchronize with the drive and can be described by an appropriate periodic ensemble. Here, we explore the consequences of breaking the time-periodic structure of the drive with additional aperiodic noise in a class of integrable systems. We show that the resulting unitary dynamics leads to new emergent steady states in at least two cases. While any typical realization of random noise causes eventual heating to an infinite temperature ensemble for all local properties in spite of the system being integrable, noise which is self-similar in time leads to an entirely different steady state, which we dub as"geometric generalized Gibbs ensemble", that emerges only after an astronomically large time scale. To understand the approach to steady state, we study the temporal behavior of certain coarse-grained quantities in momentum space that fully determine the reduced density matrix for a subsystem with size much smaller than the total system. Such quantities provide a concise description for any drive protocol in integrable systems that are reducible to a free fermion representation.


I. INTRODUCTION AND MOTIVATION
Statistical mechanics can be derived from Jaynes' principle of maximum entropy [1,2], given the constants of motion for any generic many-body system. Remarkably, this principle may also be extended to the level of a single eigenstate in a thermodynamically large, isolated quantum system, leading to the eigenstate thermalization hypothesis [3][4][5]. The eigenstate thermalization hypothesis implies that the expectation values of all operators that are local in space are equal to those obtained from a maximum entropy statistical description when appropriate conserved quantities are taken into account for constructing the ensemble. The Hamiltonian is the only such conserved quantity in a generic chaotic system. Exceptions to the eigenstate thermalization hypothesis are provided by disordered systems that show localization, both for noninteracting models [6] and with interactions [7] (i.e., many-body localized systems [8]). Integrable models provide examples of systems that are not localized and yet do not satisfy the eigenstate thermalization hypothesis. However, Jaynes' principle may again be applied here by considering the extensive number of conservation laws due to integrability that then leads to a generalized Gibbs ensemble (GGE) [9][10][11] instead of the standard Gibbs ensemble. These statistical descriptions are expected to hold only for local properties (i.e., properties that are determined by the degrees of freedom in a subsystem that is much smaller than the rest of the system, as shown schematically in Fig. 1) and not for the wave function of the entire system, which is, after all, in a pure state [12].
Tremendous experimental progress in producing and manipulating well-isolated quantum systems such as ultracold quantum gases [13][14][15][16] has also led to great interest in driven closed quantum systems. Removing the time translational invariance of closed quantum systems via a timedependent Hamiltonian leads to the much richer possibility of stabilizing nonequilibrium steady states with purely unitary dynamics. Such nonequilibrium steady states may even have properties otherwise forbidden in thermal equilibrium, such as spontaneous time-translational symmetry breaking in many-body localized systems [17,18] and dynamical topological ordering [19,20]. Many-body Floquet systems, where the time-dependent Hamiltonian satisfies with T being the time period of the drive and n being integer valued, are well-studied cases where such nonequilibrium steady states have been shown to exist. More precisely, the local properties of Floquet systems eventually synchronize with the period of the drive, which then allows the possibility of a "periodic" ensemble [21][22][23][24].
Propagating the wave function jψi of the system in a stroboscopic manner, where UðTÞ is the time-evolution operator for one time period, generates a discrete quantum map indexed by n, and the steady state emerges as n → ∞. This steady state can again be described using Jaynes' principle [21], where only quantities that are stroboscopically conserved need to be taken into account. In the nonintegrable case, it is expected that no such conserved quantity exists [8] since the Hamiltonian is no longer conserved; hence, the system should locally mimic an infinite-temperature ensemble in the steady state by entropy maximization. This is indeed consistent with the results of several recent studies [22][23][24][25][26][27][28], but a complete understanding is still lacking. The situation is markedly different for integrable systems. For systems that are reducible to free fermions, e.g., the prototypical one-dimensional (1D) spin-1=2 transverse-field Ising model, it has been shown [21] that an extensive number of quantities that are stroboscopically conserved continue to exist even in the driven case; these prevent the system from approaching the infinite-temperature ensemble locally and instead lead to a periodic GGE (p-GGE) description of the synchronized local properties at late times. The fate of other types of integrable systems like Bethe-integrable ones (for instance, the 1D spin-1=2 XXZ model) under periodic drives is still an open issue [29].

A. Synopsis
In this work, we consider integrable models that are reducible to free fermions, and our objective is to understand whether well-defined nonequilibrium steady states exist for continually driven protocols that are not periodic functions in time. We restrict ourselves to drives that can be represented by discrete quantum maps similar to Eq. (2) for ease of analysis. For concreteness, we work here with the 1D transverse-field Ising model, which is a standard manybody integrable quantum system [30][31][32]. The timedependent Hamiltonian is given by where s x;y;z are the usual Pauli operators, L denotes the number of spins in the system under consideration, and we assume periodic boundary conditions in space. Here, gðtÞ represents the time-dependent transverse magnetic field through which the system is continually driven. We start with a reference function g ref ðtÞ defined for 0 ≤ t ≤ T. Repeatedly using the corresponding timeevolution operator UðTÞ for all values of n > 0 in Eq. (2) results in periodic driving with period T. Consider now a different drive protocol where gðtÞ is constructed by patching together rescaled versions of g ref ðtÞ in time as a function of n such that the period of g ref ðtÞ is stretched or reduced by dT depending on n. The corresponding time-evolution operator is then denoted by UðT þ dTÞ or UðT − dTÞ. The new drive protocol gðtÞ can now be represented by a sequence in n that takes a value of either þ1 [representing UðT þ dTÞ] or −1 [representing UðT − dTÞ] at each n. There are an infinite number of ways to choose the sequence of AE1 as a function of n, each of which represents a different purely unitary dynamics, such that the resulting driving protocol gðtÞ cannot be represented by any periodic function in time. We note that the protocol described above involving two time periods T þ dT and T − dT is mathematically equivalent to a protocol in which the time period T is kept fixed but the Hamiltonian is scaled globally by factors of 1 þ dT=T and 1 − dT=T, respectively. The important point about the protocol is that it involves two possible time-evolution operators that do not commute with each other.
We work out the nature of the resulting nonequilibrium steady states as n → ∞ in two cases here. First, we consider a typical realization of a random sequence of AE1, where the sign is chosen with probability 1=2 at each n with dT ≪ T. This mimics a periodic drive protocol perturbed with random noise due to the lack of perfect control over the time variation of gðtÞ in an experimental setup. Second, we consider a different sequence of AE1 that is neither random nor periodic in n but has a fractal structure instead and thus can be thought of as an example of scaleinvariant noise.
For a periodic drive protocol perturbed with random noise where dT ≪ T, we show that the system, after initially approaching the p-GGE that corresponds to dT ¼ 0, eventually heats up to an infinite-temperature ensemble as n → ∞ in spite of the underlying integrability of the model. The time spent by the system close to the p-GGE, which can now be thought of as a quasisteady prethermal state, and the eventual manner of heating up to an infinite-temperature ensemble can both be understood from a simple random walk argument, provided that dT ≪ T. We further show that a coarse-graining procedure in momentum space explains how the reduced density matrix of any subsystem with l adjacent spins approaches an infinite-temperature ensemble as long as l ≪ L, even though the dynamics is completely unitary. It is well known l FIG. 1. The separation of a closed system into a subsystem with l sites (indicated in black) and the environment with L − l sites (indicated in red). Local properties only involve degrees of freedom on the sites of a subsystem that satisfies l ≪ L with L ≫ 1. that weakly interacting systems (i.e., systems close to an integrable point) after a sudden perturbation show a prethermalization regime that can be understood by a GGE based on perturbatively constructed constants of motion [33], before eventually thermalizing on a much longer time scale [34,35] with the energy being the only constant of motion. In the perturbed Floquet integrable systems with random noise that we consider, the integrability of the model is never broken, but the stroboscopically conserved quantities at dT ¼ 0 no longer remain conserved when dT ≠ 0.
However, it is important to note that the lack of conservation of these quantities does not fix the form of the final nonequilibrium steady state for an aperiodic drive protocol. For a periodic drive protocol perturbed with scaleinvariant noise, we show that if the local properties are observed not as a function of n but instead as a function of 2 n (thus, geometrically instead of stroboscopically), a new steady state emerges when n → ∞, which we call "geometric" GGE (g-GGE). Interestingly, this steady state can also be understood using an effective p-GGE construction even though the drive is not periodic in time because of a fixed-point property of products of SUð2Þ matrices. An extensive number of conserved quantities appear when the system is viewed after every 2 n drives but only when n is sufficiently large, which shows that these conservations are emergent.
The rest of the paper is arranged in the following manner. In Sec. II, we review some results relevant for our work and set the notations for the rest of the paper. Section II A details the pseudospin representation for the 1D transversefield Ising model that allows the wave function of the system to be denoted in terms of points on Bloch spheres in momentum space. Section II B details the aperiodic driving protocols that we study. Section II C reviews the way to obtain the reduced density matrix of a subsystem of l adjacent spins using appropriate correlation matrices. In Sec. III (including Secs. III A and III B), we introduce certain coarse-grained quantities in momentum space which fully determine the reduced density matrix of a subsystem when l ≪ L and thus all local quantities, and we discuss their behavior for time-periodic drive protocols. In Sec. IV, we present results for perturbed Floquet systems with random noise and use a random walk argument to show how the system heats up to an infinite-temperature ensemble in spite of its integrability. In Sec. V, we consider perturbed Floquet systems with noise that is scale invariant in time; then, we show that for a certain type of scaleinvariant noise, the system reaches a steady state that is not an infinite-temperature ensemble, but it does so only after an astronomically large number of drives. The final steady state is explicitly constructed using a fixed-point property of strings of SUð2Þ matrices. Finally, we summarize our results and conclude with some future directions in Sec. VI.

II. PRELIMINARIES
A. Pseudospin representation of dynamics for the 1D transverse-field Ising model To solve the 1D transverse-field Ising model [Eq. (3)], we introduce the standard Jordan-Wigner transformation of spin-1=2's to spinless fermions [36], Writing H in Eq. (3) in terms of these fermions, we obtain where P AE are projectors onto subspaces with an even (þ) and odd (−) number of fermions, and with antiperiodic (periodic) boundary conditions for the fermions for H þ (H − ). We now focus on the case of even L with antiperiodic boundary conditions for the fermions (which corresponds to periodic boundary conditions for the spins). To go to k space, we define e −ikn c n ; where k ¼ 2πm=L with m ¼ −ðL − 1Þ=2; …; −1=2, 1=2; …; ðL − 1Þ=2 for even L. Rewriting H þ in terms of c k and c † k , we get This Hamiltonian connects the vacuum j0i of the c fermions with the two-particle state jk; We now introduce a "pseudospin representation" σ k [which is different from the original s matrices in Eq. (3)], where j↑i k ¼ jk; −ki ¼ c † k c † −k j0i and j↓i k ¼ j0i [37]. Writing H k in Eq. (9) in this basis, we see that Then, the pseudospin state jψ k i for each k mode evolves independently as where jψðtÞi ¼⊗ k>0 jψ k ðtÞi; Thus, specifying (u k ðtÞ; v k ðtÞ) T [where the superscript T denotes the transpose of the row (u k ðtÞ; v k ðtÞ)] for each allowed k > 0 for a finite L specifies the complete wave function jψðtÞi through Eq. (12). The wave function of the system can be equivalently characterized by a two-dimensional state space at each allowed positive k, and (u k ðtÞ; v k ðtÞ) T can then be represented as a point that evolves with time according to Eq. (11) on the corresponding Bloch sphere.
In the rest of the paper, we take ðu k ; v k Þ T ¼ ð0; 1Þ T to be the starting state at each k from which the unitary dynamics proceeds. This corresponds to the system being initially prepared in the pure state where all the spins have σ x ¼ þ1 (which corresponds to the ground state for g → ∞).

B. Details of the driving protocol
While our basic conclusions are protocol independent and do not depend on the exact nature of g ref ðtÞ, we take it to have a square pulse variation in time for mathematical convenience, namely, The corresponding unitary time-evolution operator U k ðTÞ for the pseudospin mode with momentum k equals where H k has the form given in Eq. (10). For the perturbed Floquet dynamics that we consider in this work, we need to define U k ðT AE dTÞ, which follows in a straightforward manner from Eq. (14) by replacing T → T AE dT. Rewriting Eq. (2) for the perturbed Floquet dynamics for the mode with momentum k, we get where the sequence τ i ¼ τ 1 ; τ 2 ; τ 3 ; … (with each τ i being equal to either þ1 or −1) is the same for all the allowed k modes at a finite L. [In the last expression in Eq. (15), T denotes time ordering]. For Floquet systems perturbed with random noise, the sequence τ i can be taken to be any typical realization of a random process where each element is chosen to be AE1 randomly and independently with probability 1=2; we discuss the results for such a realization later in the paper. For Floquet systems perturbed with scale-invariant noise, we take the well-known Thue-Morse sequence [38][39][40], which is an infinite sequence of τ i ¼ AE1 that is obtained by starting with τ 1 ¼ −1 and successively appending the negative of the sequence obtained thus far. The first few steps of this recursive procedure yield At each recursion level m, we thus obtain the first 2 m elements of this infinite sequence. The self-similar structure of the Thue-Morse sequence can be verified by removing every second term from it, which then results in the same infinite sequence. This sequence is thus an example of a quasiperiodic sequence in time which is neither periodic nor random. Note that the average time period of this perturbed Floquet system equals T (the time period of the unperturbed Floquet system), just like in the random noise case.

C. Reduced density matrices and the distance measure
Given the wave function jψðnÞi of the entire system, which can be described by the corresponding density matrix ρ ¼ jψðnÞihψðnÞj, all local properties within a length scale l can be understood by considering the reduced density matrix of a subsystem of l adjacent spins (Fig. 1); we denote this by ρ l , which is given by where all the other L − l degrees of freedom have been integrated out. Even though ρ is a pure density matrix at any n, the reduced density matrix ρ l is mixed since the pure state jψðnÞi typically gets more entangled as n increases.
The reduced density matrix of l adjacent spins when jψðnÞi has the form given in Eq. (12) is determined in terms of the correlation functions of the c fermions at these l sites [41]. For free fermions, since all higher-point correlation functions are expressible using two-point correlations from Wick's theorem, the reduced density matrix can be determined from the knowledge of two l × l matrices [42], C and F, whose elements are constructed from u k ðnÞ and v k ðnÞ as follows: where i, j refer to sites in the subsystem. Using these expressions, we construct the 2l × 2l matrix C n ðlÞ given by The reduced density matrix ρ l is then completely fixed by the eigenvalues and eigenvectors of C n ðlÞ. For example, the entanglement entropy S ent ðlÞ equals where p k denotes the k-th eigenvalue of C n ðlÞ.
If a nonequilibrium steady state exists as n → ∞, then the convergence of the local properties to their final steadystate values as n increases can be characterized by defining a distance measure D n ðlÞ as We note that 0 ≤ D n ðlÞ ≤ 1, and it is identically zero only if C n ðlÞ ¼ C ∞ ðlÞ. For a periodically driven protocol with dT ¼ 0, C ∞ ðlÞ can be calculated using the appropriate p-GGE; the distance measure decays to zero as D n ðlÞ ∼ n −αðTÞ for large n (and l ≪ L), with αðTÞ ¼ 1=2 or 3=2 depending on the time period T of the drive [43]. We note that D n ðlÞ does not decay to zero as n → ∞ if l=L is finite; this explicitly shows that only local quantities can be described by the p-GGE.

III. COARSE GRAINING IN MOMENTUM SPACE
Even though the entire wave function jψðnÞi needs the specification of (u k ðnÞ; v k ðnÞ) T at k ¼ 2πm=L with m ¼ 1=2; 3=2; …; ðL − 1Þ=2, C n ðlÞ, and consequently, the reduced density matrix for any l ≪ L depends only on certain coarse-grained variables defined as follows [44,45]: These variables are defined using N c ð≫1Þ consecutive k modes that lie within a cell (k cell ), which has an average momentum k c in momentum space and a size δk where 1=L ≪ δk ≪ 1=l. With this condition on δk, the sinusoidal factors in Eq. (18) can be replaced, to a very good approximation, by for all the k modes within the cell. The sum over momentum in Eq. (18) can thus be carried out in two steps-first summing over the consecutive momenta in a single cell, and then summing over the different momentum cells. This gives where N cell ≪ L=2 represents the number of momentum cells after the coarse-graining procedure. Note that even though quantities like ju k ðnÞj 2 and u Ã n ðkÞv n ðkÞ continue to display Rabi oscillations and thus have no well-defined n → ∞ limit, the coarse-grained variables ðju k ðnÞj 2 Þ c and (u Ã k ðnÞv k ðnÞ) c must reach final steady-state values for a nonequilibrium steady state to exist.
If the system is locally described by an infinite-temperature ensemble, then ρ l ¼ I, where I is the identity matrix, for any l ≪ L. This immediately implies that A similar behavior has been found for repeated quenching of the transverse field in the 1D spin-1=2 XY model [46].) We calculate the n dependence of such coarse-grained quantities numerically by taking a sufficiently large system size L (typically L ¼ 131 072, but L ¼ 524 288 in some cases), so that the allowed momentum modes are sufficiently dense in the interval k ∈ ½0; π, and then dividing the Brillouin zone into N cell ¼ 32 cells, each of which has N c ≫ 1 consecutive momenta. For the periodic drive case where dT ¼ 0, we see that although ðju k ðnÞj 2 Þ c indeed converges to a finite value (Fig. 2) and thus has a well-defined nonequilibrium steady state, this constant depends on the value of the average momentum of the coarse-grained cell (k c ) and is different from 1=2, showing that the local description is different from that of an infinite-temperature ensemble.

A. Probability distribution on the Bloch sphere
At each k, the pseudospin state jψ k ðnÞi obtained after n drives can be parametrized in terms of three angles (θ k ðnÞ; ϕ k ðnÞ; β k ðnÞ), The overall phase β k ðnÞ is not important, and we will ignore it henceforth. We concentrate on the angles θ k ðnÞ and ϕ k ðnÞ, which define the motion of a point on the unit Bloch sphere for each momentum k. From Eq. (14), it is clear that U k ðTÞ is an element of the group SUð2Þ (namely, 2 × 2 unitary matrices with determinant equal to 1), and it can be written as where 0 ≤ γ k ≤ π andê k ¼ ðe 1k ; e 2k ; e 3k Þ is a unit vector. For a drive protocol that is perfectly periodic, (θ k ðnÞ; ϕ k ðnÞ) traces out a circle on the Bloch sphere formed by the intersection of this unit sphere, with a plane that contains the point (θ k ð0Þ; ϕ k ð0Þ) and whose normal vector equalsê k . It is then clear that the quantities I k defined as are all independent of n for a periodic drive protocol, and hence they define the extensive number of stroboscopically conserved quantities. The behaviors of the coarse-grained quantities in Eq. (22) depend on the simultaneous positions of (θ k ðnÞ; ϕ k ðnÞ) on the unit sphere of all the momentum modes that lie within a coarse-grained cell. Since their number N c ≫ 1 when L ≫ 1, it is useful to instead look at the probability distribution of these N c points on the unit sphere. In particular, we study the probability distribution of ju k ðnÞj 2 − jv k ðnÞj 2 ¼ cos θ k for each such coarsegrained momentum cell, which we denote as P n ðcos θ k Þ, as a function of n. We can also define another probability distribution from the motion of a single k mode for a large number of drive cycles and denote it byPðcos θ k Þ. Remarkably, we find that for large n, where k c denotes the average momentum of a coarse-grained cell that has N c momenta on the left-hand side of Eq. (28), and the right-hand side is obtained from the motion of a single momentum mode at k c . For an infinite-temperature ensemble, we expect that the N c points will cover the unit sphere uniformly at large n even though all these points started from the south pole at n ¼ 0 because of the choice of jψð0Þi. Thus, P n ðcos θ k c Þ at large n should be independent of the value of the average momentum of the cell, k c , and thus equal 1=2, if the local description is that of an infinite-temperature ensemble. We derivePðcos θ k c Þ for the perfectly periodic case in the next subsection and see that, numerically, P n ðcos θ k c Þ approachesPðcos θ k c Þ for large n. It turns out that Pðcos θ k c Þ depends on k c ; this again shows that the description is not locally that of an infinite-temperature ensemble for periodic driving. B. Probability distribution for periodic driving: dT = 0 We first derive the form ofPðcos θ k Þ when there is no randomness, i.e., when dT ¼ 0. In this case, the time-evolution operator after one time period T for the momentum mode k (where 0 ≤ k ≤ π) is given by Since U k ðTÞ is an element of the group SUð2Þ, it can be written in the form of Eq. (26), where γ k and the unit vector e k ¼ ðe 1k ; e 2k ; e 3k Þ can be derived by using the expression given in Eq. (29). Then, the two-component spinor obtained after n drives is given by Hence, where we have used the identity e 2 1k þ e 2 2k þ e 2 3k ¼ 1. We now assume that γ k =π is an irrational number; this is justified since rational numbers lying in the range [0, 1] form a set of measure zero. Hence, the set of numbers 2nγ k mod π covers the region ½0; π uniformly as n → ∞. We now use the fact that if 2nγ k is uniformly distributed from 0 to π, the probability distribution of the variable v ¼ cosð2nγ k Þ is given by Þ.] Using Eq. (32) in Eq. (31), we find that the probability distribution of cos θ k is given bȳ for −1 ≤ cos θ k ≤ 1À2e 2 3k , whilePðcos θ k Þ ¼ 0 for 1 − 2e 2 3k < cos θ k ≤ 1. We observe that this distribution has square-root divergences at two points given by cos θ k ¼ −1 and 1 − 2e 2 3k , and the range of cos θ k is given by 2ð1 − e 2 3k Þ. Furthermore, the distribution is correctly normalized so that R 1 −1 dðcos θ k ÞPðcos θ k Þ ¼ 1. Given the values of the parameters g i , g f , T, and k, one can calculateê k and therefore e 3k from Eqs. (29) and (26). Equation (33) then gives the probability distribution of cos θ k , which is generated by a large number of drives, as shown in Fig. 3. We also calculate P n ðcos θ k Þ for the 1D transverse-field Ising model with L ¼ 131072, where the positive momenta in the Brillouin zone are divided into N cell ¼ 32 cells, each thus containing N c ¼ 2048 consecutive momenta; we display the result for a large n (¼5 × 10 5 ) in Fig. 3, which validates Eq. (28).
Equation (28) can be analytically understood as follows. For large n, we have argued above that the right-hand side of Eq. (28) is given by Eq. (33) with k ¼ k c . Now, we look at the left-hand side of Eq. (28); here, we have to consider a large but fixed value of n and average over a range of momenta δk around k ¼ k c . We begin with the expressions U k ðTÞ ¼ expð−iγ kêk · ⃗σÞ; For δk ≪ 1, γ k will be close to γ kþδk andê k will be close tô e kþδk . However, (U k ðTÞ) n ¼ expð−inγ kêk · ⃗σÞ; (U kþδk ðTÞ) n ¼ expð−inγ kþδkêkþδk · ⃗σÞ; ð35Þ and it is clear that if n is very large, nγ k and nγ kþδk will not be close to each other. In other words, a small value of γ kþδk − γ k gets amplified as n increases but a small value of e kþδk −ê k does not get amplified. For δk ≪ 1, we can therefore make the approximation of settingê k ¼ê kþδk , but we cannot set γ k ¼ γ kþδk if we are interested in the n → ∞ limit. Hence, imply, using Eq. (31), that Now, we can write γ kþδk − γ k ¼ g k δk, where g k is a number of order 1. Suppose that n is large enough that 2ng k δk ≫ π. This implies that as k goes over all the momenta in a cell lying in the range k c − δk=2 to k c þ δk=2, 2nγ k mod π will cover the range ½0; π uniformly. Hence, in this cell with an average momentum k c , the distribution of the variable v ¼ cosð2nγ k Þ will again be given by Eq. (32). Arguments similar to the ones leading from Eq. (31) to Eq. (33) then show that the left-hand side of Eq. (28) is also given by Eq. (33). We now derive the value to which ðju k ðnÞj 2 Þ c converges as n → ∞ (Fig. 2). It is straightforward to see that 1=2 − ðju k ðnÞj 2 Þ c ¼ −ð1=2Þ( cos θ k ðnÞ) c , where the coarsegrained quantity ( cos θ k ðnÞ) c is given by As n → ∞, using Eq. (28) and Eq. (33), we get This is in agreement with the results shown in Fig. 2 for 1=2 − ðju k ðnÞj 2 Þ c , which, in the limit n → ∞, is equal to ð1=2Þe 2 3k c .

IV. RESULTS FOR PERTURBED FLOQUET SYSTEM WITH RANDOM NOISE
We now study the unitary dynamics when dT ≪ T and the τ i 's in Eq. (15) are taken to be any typical realization of a random sequence of AE1's which are chosen with probability 1=2 each. We again look at the behavior of the coarse-grained quantities to understand the nature of the resulting nonequilibrium steady state. We see that in spite of the integrability of the model, the resulting nonequilibrium steady state is locally described by an infinite-temperature ensemble as n → ∞, unlike the case with dT ¼ 0.
The first quantity we look at is the overlap of a reference k mode jψ k ðnÞi, with other k þ δk modes that lie very close to it, i.e., where δk ≪ 1, as a function of n when the unitary dynamics is of the form given in Eq. (15), with the τ i 's taken from a typical realization of a random sequence of AE1's. Figure 4 shows a plot of jhψ k ðnÞjψ kþδk ðnÞij 2 versus n. At n ¼ 0, all the jψ n ðkÞi start from the same state ð0; 1Þ T . We see that while jhψ k ðnÞjψ kþδk ðnÞij 2 stays localized between 1 and another value (which is derived in the next paragraph) for the perfectly periodic drive with dT ¼ 0 (see Fig. 4), the dT ≠ 0 case is extremely sensitive to initial conditions; we then find that jhψ k ðnÞjψ kþδk ðnÞij 2 covers the entire range [0, 1] at large n when the overlap with different k þ δk modes is considered, where all δk ≪ 1 (see Fig. 4). Note that if δk ≪ 1, then each individual U kþδk ðT þ τ n dTÞ is only slightly different from U k ðT þ τ n dTÞ of the reference k mode. This behavior of the overlaps already suggests that the coarse-grained quantities behave completely differently in the dT ≠ 0 case at large n.
For the perfectly periodic drive with dT ¼ 0, we can understand the behavior of jhψ k ðnÞjψ kþδk ðnÞij 2 versus n as follows. Following the arguments given in Eqs. (34)-(36), we find that jhψ k ðnÞjψ kþδk ðnÞij 2 For n ≫ 1=jγ k − γ kþδk j, sin 2 ½nðγ k − γ kþδk Þ will oscillate between 0 and 1. Equation (40) then implies that ; v k ð0Þ) T ¼ ð0; 1Þ T for all momenta. Here, k ¼ π=2 is taken as the reference mode, and its overlap calculated with the next ten consecutive momenta is shown, where the momentum spacing equals π=65 536. jhψ k ðnÞjψ kþδk ðnÞij 2 will oscillate between 1 and e 2 3k , and the range of values is given by 1 − e 2 3k . We note that this is exactly half the range of values of cos θ k , which is given by 2ð1 − e 2 3k Þ according to Eq. (33). We now study what happens to the probability distribution P n ðcos θ k Þ for the unitary process with dT ≠ 0. To evaluate this quantity, we use a chain with L ¼ 524 288 spins, where each coarse-grained cell in momentum space contains N c ¼ 8192 consecutive momenta (hence, N cell ¼ 32) and the spinor at each k equals ð0; 1Þ T at n ¼ 0. The results for one such cell (with average momentum k c ¼ 31π=64) are shown in Fig. 5. For small n, we see that P n ðcos θ k Þ approaches the distribution for the perfectly periodic drive with the same T but with dT ¼ 0, which has characteristic square-root divergences as shown in Eq. (33). However, as n increases further, the probability distribution starts deviating strongly from this form and approaches a constant distribution eventually at large n (see Fig. 5); this shows that the N c points that constitute the coarse-grained momentum cell are uniformly spread over the unit sphere. The lower panel of Fig. 5 shows the positions of the N c points on the unit sphere, from which P n ðcos θ k Þ is constructed, for n ¼ 100, n ¼ 10 000, and n ¼ 200 000. We now provide an understanding of this spreading based on the idea of random walks.
We define two unitary matrices U k ðT − dTÞ and U k ðT þ dTÞ given by These are the time-evolution operators for the two possible time periods T − dT and T þ dT, respectively. The time-evolution operator U nk for n drives is obtained by multiplying n matrices, each of which is randomly chosen to be either U k ðT − dTÞ or U k ðT þ dTÞ with probability 1=2 each to obtain a particular realization of this random process as shown in Eq. (15).
Next, we parametrize the two matrices as For dT ¼ 0, these matrices reduce to U k ðTÞ given in Eq. (29); in that case, we have seen that the evolution operator U k ðTÞ corresponds to a rotation about a unit vector e k given in Eq. (26). If dT is nonzero but small, the unit vector corresponding to each drive will have a small random component. To quantify this, we note that if dT is small, the unit vectorsâ k andb k given in Eq. (42) differ by small angles Δϕ 1k ¼ arccosðê k ·â k Þ and Δϕ 2k ¼ arccosðê k ·b k Þ from the unit vectorê k . These two small angles are of the same order, and we can denote their average as Δϕ k . For dT ¼ 0 or sin k ¼ 0, it is clear from Eqs. (29) and (41) that the unit vectorsê k ,â k , andb k are identical, and therefore Δϕ k ¼ 0. Hence, we see that when either dT or sin k → 0, and Eq. (43) holds independently of the exact form of the drive protocol. If the drives have time periods randomly given by T þ dT and T − dT, we expect from the theory of random walks that after n drives, the unit vector corresponding to U nk will differ fromê k by an angle proportional to ffiffiffi n p Δϕ k , assuming that Δϕ k is small (this will be true if dT is small). However, since the unit vectors lie on a compact space (given by the points on the surface of the unit sphere), the deviation of the unit vector of U nk fromê k cannot go to ∞ as n → ∞; rather, the unit vectors will cover the unit sphere uniformly giving rise to a uniform probability distribution Pðcos θ k Þ ¼ 1=2. The behavior of the probability distribution P n ðcos θ k Þ shown in Fig. 5 for large n is thus consistent with Eq. (28). Furthermore, from the above argument, the spreading of the overlaps jhψ k ðnÞjψ kþδk ðnÞij 2 as a function of n (as shown in Fig. 4) when many neighboring momentum modes very close to the reference mode k are considered is controlled by a time scale nT, where nðΔϕ k Þ 2 ≃ nðdTÞ 2 ðsin kÞ 2 ∼ Oð1Þ; hence, the time scale nT varies with the momentum k. We explicitly verify this by showing the behavior of the corresponding coarse-grained variables ðju k ðnÞj 2 Þ c [defined in Eq. (22)] in Fig. 6. The data are consistent with an exponential decay to 1=2; however, the rate of decay is clearly sensitive to the values of both dT and k c (the average momentum of the coarse-grained cell) as we see in Fig. 6(a). However, when the data are plotted versus nðdTÞ 2 ðsin k c Þ 2 [see Fig. 6(b)] instead of n, we clearly see that the time scale of the exponential decay to 1=2 is controlled by τ k c ;dT ¼ 1=½ðdTÞ 2 ðsin k c Þ 2 as predicted by the random walk argument. We have further checked that the coarse-grained quantity (u Ã k ðnÞv k ðnÞ) c [Eq. (22)] relaxes to zero for large n for any value of k c .
The behaviors of both P n ðcos θ k Þ in Fig. 5 and ðju k ðnÞj 2 Þ c in Fig. 6 as functions of n show the irreversible approach to an infinite-temperature ensemble when local quantities are probed. Finally, we consider the implication of the dependence of the time scale τ k;dT ¼ 1=(ðdTÞ 2 ðsin k c Þ 2 ) [which controls the relaxation of ðju k ðnÞj 2 Þ c ] on dT and k c for the rate at which the trace distance of the reduced density matrix of a subsystem of l consecutive spins from an infinite-temperature ensemble [using the definition D n ðlÞ in Eq. (21)] approaches zero after a large number of drives n. Since the trace distance receives a contribution from all momenta k, it should approach zero as where C ∞ ðlÞ in Eq. (21) is calculated using an infinitetemperature ensemble. For large n, the integral in Eq. (44) is dominated by the regions in k where τ k;dT is large, namely, k ¼ 0 and π where sin k → 0. Absorbing n in k (or π − k) in the exponential, we see that the trace distance will scale with n as 1=ð ffiffiffi n p dTÞ when dT is small. Thus, D n;ITE ðlÞ ∼ F l ðnðdTÞ 2 Þ, where F l ðxÞ ∼ 1= ffiffi ffi x p for x ≫ 1. Second, F l ðxÞ ∼ Oð1Þ when x ≪ 1 since in this limit the system must relax to the p-GGE that emerges for dT ¼ 0. Thus, the prethermal regime where the system initially approaches the p-GGE exists for time scales that depend on dT as n ∼ 1=ðdTÞ 2 for dT ≪ T. The numerical data for dT ≠ 0 for the behavior of D n;ITE ðlÞ as a function of x ¼ nðdTÞ 2 are consistent with both these expectations, i.e., F l ∼ 1= ffiffi ffi x p for x ≫ 1 and F l ∼ Oð1Þ for x ≪ 1, as shown in Fig. 7 for a fixed T ¼ 0.2 and different values of dT. It is also instructive to look at D n;GGE for nonzero values of dT (Fig. 7, inset) where C ∞ ðlÞ is calculated using the dT ¼ 0 p-GGE as a function of n. For time scales n ∼ 1=ðdTÞ 2 , D n;GGE decreases as a function of n, which shows that the system approaches the p-GGE in the prethermal regime. However, at larger values of n, D n;GGE starts increasing as a function of n, which explicitly shows that the large-n steady state is not described by a p-GGE. The initial relaxation to a p-GGE in the prethermal regime scales either as D n;GGE ∼ n −3=2 (the case shown in Fig. 7, inset) or as D n;GGE ∼ n −1=2 depending on which dynamical phase the periodic drive protocol (dT ¼ 0) belongs to [43].
For completeness, we point out that Ref. [47] (see also Ref. [48]) also considered a noisy 1D transverse field Ising model, though not periodically driven on average, and showed that the asymptotic steady state is an infinitetemperature ensemble. However, they considered the behavior of quantities averaged over different realizations of noise, which resembles an open quantum system with nonunitary dynamics. Here, we show that if one is restricted to understanding only local quantities as functions of n in a perturbed Floquet integrable system, even the unitary dynamics of a single typical realization of a random sequence leads to an infinite-temperature ensemble.

V. RESULTS FOR PERTURBED FLOQUET SYSTEM WITH SCALE-INVARIANT NOISE
We now study the unitary dynamics when dT ≠ 0, where the τ i 's in Eq. (15) are not a random sequence of AE1 but are instead given by the scale-invariant Thue-Morse sequence in Eq. (16), which is neither periodic nor random. The average time period of such a drive protocol continues to equal T. To this end, we first show the coarse-grained quantity ðju k ðnÞj 2 Þ c in Eq. (22) as a function of n, where the τ i 's are chosen to be AE1 according to the Thue-Morse sequence, with T ¼ 0.2 and dT ¼ 0.05 (Fig. 8). We see that these quantities behave entirely differently from either the periodic case (dT ¼ 0) in Fig. 2 or the randomly perturbed case in Fig. 6. From Fig. 8, it is not even clear whether a well-defined nonequilibrium steady state exists as n → ∞. We show below that, in fact, a nonequilibrium steady state does exist in the large-n limit if the local quantities are observed not as a function of n but instead as 2 n (hence, geometrically instead of stroboscopically). New emergent conserved quantities, similar to Eq. (27) when 2 n is used instead of n, appear in this unitary dynamics but only when n becomes sufficiently large, which then allows for the construction of a "geometric" GGE for the resulting steady state.
We begin by noting that the evolution operators for drives corresponding to T − dT and T þ dT at a particular k are, respectively, given by the matrices U k ðT − dTÞ and U k ðT þ dTÞ given in Eqs. (41); for simplicity of notation, we will not write the subscript k in this section henceforth. Furthermore, we denote UðT − dTÞ by A 0 and UðT þ dTÞ by B 0 here. Then, the operators corresponding to the pair of drives ðτ i ; τ iþ1 Þ ¼ −1; þ1 and ðτ i ; τ iþ1 Þ ¼ þ1; −1 are given by A 1 ¼ B 0 A 0 and B 1 ¼ A 0 B 0 , respectively. We now focus on the case where n ¼ 2 m , with m ¼ 0; 1; 2; Á Á Á. It is then easy to see that the unitary dynamics can be entirely expressed in terms of the new matrices A 1 and B 1 . We illustrate this in the second line of Eq. (45) for the case of m ¼ 4 (thus, n ¼ 16). This motivates us to define matrices A m and B m recursively as for all m ≥ 0. We can then show that the evolution operator after exactly 2 m drives is given by A m . [Equation (45) illustrates this for m ¼ 4, i.e., n ¼ 16.] Note that this provides a very efficient method to calculate the dynamics after n ¼ 2 m , which only requires OðmÞ matrix multiplications at each k instead of Oð2 m Þ. We show that, remarkably, A m ¼ B m when m becomes sufficiently large, and thus, there is an emergent periodicity when the system is viewed geometrically (i.e., after every n ¼ 2 m ). We now study how the matrices A m , B m defined in Eqs. (46) behave in the limit m → ∞. To this end, we parametrize these matrices as We also define the angle ϕ m between the unit vectorsâ m andb m , namely, with 0 ≤ ϕ m ≤ π. Using Eqs.
From Eqs. (49), we can prove the following results. First, we find that α m ¼ β m for all m ≥ 1. We therefore use only the variable α m henceforth. Next, we find the recursion relations Since α m ¼ β m for all m ≥ 1, Eqs. (49) imply that This means that the unit vector that points in the direction ofâ m þb m is the same for m and m þ 1, up to a sign; hence, up to a sign, the direction of this unit vector does not change at all with m, all the way from m ¼ 1 to ∞. These arguments will break down if eitherâ m orb m is ill defined; according to Eq. (47), this can only happen if α m or β m becomes exactly equal to 0 or π. However, these cases require very special choices of A 0 and B 0 , which form a set of measure zero; we therefore ignore such special cases.
To summarize, we find that the dynamics given by Eq. (46) conserves two quantities at every iteration: α m ¼ β m , and the direction of the unit vectorâ m þb m is conserved for all m ≥ 1. The existence of these conserved quantities considerably simplifies the analysis as we will see below.
For sufficiently large m, we find numerically that ϕ m → 0, although ϕ m does not approach zero monotonically (see Fig. 9). We can understand this as follows. If ϕ m → 0, Eq. (50) implies that α mþ1 → 2α m mod π. Assuming that α m =π is irrational, we see that α m will cover the region ½0; π uniformly as m → ∞. Next, assuming ϕ m , ϕ mþ1 to be small, we obtain from the lowest-order expansion of the terms appearing in Eq. (51). Since α m covers ½0; π uniformly, j tanðα m Þj varies all the way from zero to ∞; this explains why ϕ m does not monotonically approach zero. In fact, logðj tan α mþj jÞ; ð54Þ and for large N, a uniform distribution of α m implies that P N−1 j¼0 logðj tan α mþj jÞ → ðN=πÞ R π 0 dα logðj tan αjÞ, which is equal to zero. Thus, we would expect that as m → ∞, ϕ m should stay at a constant value if we stopped at the firstorder expression in Eq. (53). We therefore have to go to the next order in the expansion of Eq. (51). We find that The fact that the last factor on the right-hand side is always less than 1 explains why ϕ m eventually goes to zero as m → ∞; this explains the numerical behavior shown in Fig. 9.
As m → ∞, the fact that ϕ m → 0 means that a m −b m → 0. The discussion following Eq. (52) then implies thatâ ∞ ¼b ∞ is equal, up to a sign, to the unit vector that points in the direction ofâ 1 þb 1 . Hence, the value ofâ ∞ ¼b ∞ can be found right in the beginning when we calculateâ 1 andb 1 .
We now examine what happens if we view the system geometrically after n ¼ 2 m drives where m ¼ 0; 1; 2; Á Á Á. In Fig. 9 (inset), we show the components of the unit vector a m versus m for g  (Fig. 9, main panel). We therefore conclude that after an extraordinarily large number of drives given by 2 m , where m ≳ 1300 (we have checked this for other values of k also),â m andb m become identical and independent of m.
Since the unitary dynamics for n ¼ 2 m can be represented by a single unitary matrix A m for the Thue-Morse sequence, this immediately shows that J k defined as equals J k ≡ hψ k ð0Þjâ ∞ · ⃗σjψ k ð0Þi and becomes independent of m for sufficiently large m at each momentum k. We should stress here that, unlike I k defined in Eq. (27) for the dT ¼ 0 case which is conserved as a function of n, the quantities J k are conserved geometrically, i.e., when n ¼ 2 m . Furthermore, one requires m to be larger than a threshold value (m ≳ 1300 for the drive parameters shown in Fig. 9) for these conserved quantities to emerge at all momenta k, and these conserved quantities are not manifest in the unitary dynamics below this threshold value of m. We have also verified numerically thatâ m andb m become identical and independent of m at large m not only for dT ≪ T but for any value of dT=T. The motion of the N c ð≫ 1Þ consecutive momenta contained in any of the coarse-grained momentum cells on the unit sphere as a function of n provides another way to geometrically see the emergent conserved quantities of this unitary dynamics (with the τ i 's chosen according to the Thue-Morse sequence). As shown in Fig. 10, when the system is viewed geometrically, i.e., with n ¼ 2 m , the N c points that belong to a momentum cell with average momentum k c initially perform a complicated dynamics on the unit sphere but eventually settle down to a very simple motion along a circle in the plane defined byâ ∞ (at k c ) for large enough values of m beyond a particular threshold.
Finally, we show in Fig. 11 the components of δn k 1 a k ðT; dTÞ −â k ðT; 0Þ versus k when bothâ m andb m have settled down to the same value for a unitary process with the τ i 's chosen according to the Thue-Morse sequence for a given T and dT. Sinceâ k ðT; 0Þ represents the unit vector that corresponds to the unitary matrix U k ðTÞ for the perfectly periodic case (dT ¼ 0) given in Eqs. (26) and (29), we must have δn k → 0 as dT → 0. In Fig. 11(a), we show the variation of δn k versus k for dT ¼ 0.05; T ¼ 0.2, and the corresponding results for dT ¼ 0.8, T ¼ 0.2 (thus, dT > T) in Fig. 11(b). We see that indeedâ ∞ ð¼b ∞ Þ equals (up to a sign) the unit vector in the direction ofâ 1 þb 1 . The nature of the variation with k is rather different for dT ≪ T [ Fig. 11(a)] and dT ≫ T [ Fig. 11(b)]. Also, we note that δn k equals 0 for k ¼ 0 and k ¼ π as expected since ja 3 j ¼ 1 and a 1 ¼ a 2 ¼ 0 independently of the value of dT at these two momenta.
Knowing the value ofâ ∞ ¼b ∞ as a function of k (Fig. 11) allows us to construct the local description of the final nonequilibrium steady state. The construction is completely analogous to the perfectly periodic driven case (dT ¼ 0), where the relevant integrals of motion are now J k [in Eq. (56)] at each momentum instead of I k [in Eq. (27)]. The density matrix of the g-GGE then equals where the Lagrange multipliers λ k are fixed by the condition and the normalization constant Z ensures that Tr½ρ g-GGE ¼ 1. The final steady-state value of any local operator can then be evaluated by using this ensemble.
We show the results for one such local quantity hψðnÞjs x jψðnÞi in Fig. 12. Figure 12(a) shows the convergence of this quantity to the final steady-state value calculated using ρ g-GGE for two different cases (dT ≪ T and dT ≫ T) as a function of n ¼ 2 m ; the convergence to the final g-GGE is extremely slow, as we have already discussed in this section. Figure 12(b) shows the data for the case dT ≪ T as a function of n from which it is clear that even after a large number of drives n ∼ 10 5 , the system does not seem to relax to the final nonequilibrium steady state in any simple manner if dT ≠ 0 [whereas the relaxation of local quantities at dT ¼ 0 is much faster,  as is clear from the inset of Fig. 12(a)], and the relaxation to the nonequilibrium steady state takes place over a much longer time scale, as we see in Fig. 12(a). We note here that this extremely slow relaxation of local quantities to their final steady state appears to be quite insensitive to the value of dT as long as it is nonzero. Figure 12(b) also shows the perfect agreement of the data calculated at n ¼ 2 m using the unitary matrix A m at each k to the results of the more direct calculation at each n [using Eq. (15)] without using the recursive structure of the Thue-Morse sequence.

VI. CONCLUSIONS AND OUTLOOK
In this work, we have considered a prototypical integrable model, the one-dimensional transverse-field Ising model, which is continually driven with a time-dependent Hamiltonian that results in a unitary dynamics. Since the model is translation invariant, the modes can be labeled by their momentum k. For a subsystem whose size is much smaller than that of the total system (shown schematically in Fig. 1), the reduced density matrix of the subsystem is controlled by certain coarse-grained quantities in momentum space as defined in Eq. (22). Even though the corresponding quantities for a single momentum continually fluctuate as n → ∞ and do not have a well-defined large-n limit, these coarse-grained quantities have a welldefined n → ∞ limit if a nonequilibrium steady state exists. Furthermore, we define a probability distribution on the unit sphere at each average momentum k c by employing a pseudospin representation at each k and using the large number of momenta lying in a small cell in momentum space centered around a momentum k c during the coarsegraining procedure. This probability distribution at each k c changes in an irreversible manner to the final n → ∞ distribution as a function of n. At large n, this distribution generated from the points in a cell (centered at k c ) at that n has the property stated in Eq. (28) that it equals the probability distribution generated from the motion of the pseudospin at a single momentum k c (thus, with no coarse graining in momentum) for a large number of drives.
We first recapitulate the results known for the wellstudied case in which the system is driven perfectly periodically with a time period T. Since the model is completely integrable with an infinite number of conserved quantities I k [defined in Eq. (27)], the system does not locally heat up to infinite temperature; instead, this driving leads to the periodic generalized Gibbs ensemble. We also calculate the behavior of the coarse-grained quantities in momentum space and the probability distributions on the unit sphere as a function of k c , which immediately shows that the nonequilibrium steady state is not an infinitetemperature ensemble.
We then consider two new driving protocols, both of which involve driving with two possible periods T þ dT and T − dT, such that the unitary time-evolution operators corresponding to the two periods do not commute with each other; these driving protocols can be thought of as perturbations of the perfectly periodic driving if dT ≪ T. In the first driving protocol, the successive drive periods are randomly chosen to be T þ dT and T − dT, with probability 1=2 each. In this case, we find that after a large number of drives, the system eventually approaches an infinite-temperature ensemble. Remarkably, this ensemble also has the property that the probability distribution of a coarse-grained quantity averaged over many momenta in a cell is equal to the probability distribution of the same quantity generated by many drives for a single momentum. The second driving protocol that we consider is one in which the time periods T þ dT and T − dT are chosen according to the Thue-Morse sequence. This sequence is neither periodic nor random, but it has the property of being scale invariant if the number of drives is doubled. In this case, we find that the system approaches a nonequilibrium steady state with conserved quantities J k [defined in Eq. (56)], which are smooth functions of k if we view it after n ¼ 2 m drives rather than after n drives. We call this steady state a geometric generalized Gibbs ensemble; to the best of our knowledge, this is the first known example of such an ensemble. We note, however, that this ensemble emerges only when m is large (of the order of 1300 for one particular set of parameters), which implies that the number of drives n required is astronomically large.
An interesting problem for future studies may be to understand how the local properties relax in time to the final geometric generalized Gibbs ensemble that we have discussed here. Such a study would lead to a deeper understanding of the fixed point of a sequence of SUð2Þ matrices, which is generated from two arbitrary initial matrices A 0 and B 0 (which do not commute with each other) by the recursion relations given in Eq. (46). The numerical results that we have presented in Sec. V provide only a glimpse of this general problem, which may possibly involve a very rich mathematical structure.
More generally, it will be interesting to understand the conditions that determine whether a continually driven integrable system locally heats up to an infinite-temperature ensemble or not [49]. For the latter case, it will be instructive to find other examples of quasiperiodic drive protocols [50], apart from the one we have pointed out here, where well-defined nonequilibrium steady states may emerge without the restriction of a time-periodic drive.