Long-lived interacting phases of matter protected by multiple time-translation symmetries in quasiperiodically-driven systems

We show how a large family of interacting nonequilibrium phases of matter can arise from the presence of multiple time-translation symmetries, which occur by quasiperiodically driving an isolated quantum many-body system with two or more incommensurate frequencies. These phases are fundamentally different from those realizable in time-independent or periodically-driven (Floquet) settings. Focusing on high-frequency drives with smooth time-dependence, we rigorously establish general conditions for which these phases are stable in a parametrically long-lived ‘preheating’ regime. We develop a formalism to analyze the effect of the multiple time-translation symmetries on the dynamics of the system, which we use to classify and construct explicit examples of the emergent phases. In particular, we discuss time quasi-crystals which spontaneously break the time-translation symmetries, as well as time-translation symmetry protected topological phases.


I. INTRODUCTION
Many-body quantum systems give rise to a vast array of interesting phases of matter. The last decade has seen a dramatic expansion and refinement in our understanding of the landscape of such phases [1][2][3][4]. Recently, it was understood how time-translational symmetry (TTS) itself can give rise to and protect intrinsically out-ofequilibrium phases of matter in isolated quantum systems [5][6][7][8]. Arguably the simplest example is the discrete time crystal [7,[9][10][11][12][13][14], characterized by the spontaneous breaking of the discrete TTS of a periodic (Floquet) drive. This is manifested in physical observables exhibiting robust, long-lived oscillations at an integermultiple of the base driving period. Experimental signatures of this behavior have been reported in various platforms [15][16][17]. Going beyond discrete time crystals, a large number of other Floquet phases in which TTS plays an essential role, including topological phases beyond the equilibrium classification, have also been discussed.
The richness of Floquet phases naturally raises the questions: are there fundamentally different nonequilibrium interacting phases beyond Floquet, which are not dynamically engineered analogs of static phases? What is the role of TTS in characterizing these phases?
Apart from theoretical interest, these questions are placed upon us by the dramatic experimental advances in controlling and manipulating isolated quantum manybody systems, such as cold atoms [18], trapped ions [19][20][21], nitrogen-vacancy centers [22,23], and superconducting qubits [24,25]. These systems provide a natural platform to realize dynamical protocols, allowing us to systematically study physics in out-of-equilibrium settings, * d else@mit.edu † wenweiho@fas.harvard.edu ‡ pdumitrescu@flatironinstitute.org including thermalization and equilibriation. Classifying nonequilibrium phases tells us exactly what long-time, dynamical collective behaviors are possible and the universal features defining them. Driven interacting systems are, however, generically expected to heat up to a featureless infinite temperature state due to a lack of energy conservation [26][27][28][29][30]. Thus, to meaningfully define phases of matter in such settings, systems must be protected against heating, at least for some long timescale. For Floquet systems, this challenge can be overcome by applying high-frequency drives leading to exponentially long-lived prethermal regimes [31][32][33][34][35] or by applying strong disorder leading to many-body localization (MBL) [30,36,37]. Generalizing these ideas to more generic driving scenarios remains an important open question.
In this paper, we consider interacting quantum manybody systems subject to a quasiperiodic drive that consists of several incommensurate frequencies and is smooth in time. We rigorously show that under such driving scenarios, the system is protected at high driving frequencies from heating for a parameterically long time, giving rise to a so-called 'preheating' regime. The heating time scales as a stretched exponential of the ratio of the drive frequency to local coupling strengths. We demonstrate through a recursive construction that there is an effective static Hamiltonian governing time-evolution in the preheating regime, which generalizes the Floquet analysis of [12,35] to a large class of new dynamical systems.
The presence of a preheating regime in quasiperiodically-driven systems opens up an avenue to realize novel long-lived, nonequilibrium phases of matter. We provide a set of general driving conditions to realize these phases and discuss how they are distinguished by a notion of multiple time-translation symmetries (TTSes) of the drives; thus, they are fundamentally different from those in static or Floquet settings. In particular, we classify two exemplars of such phases: the discrete time quasi-crystal (DTQC), which spontaneously breaks the multiple TTSes, as well as quasiperiodic symmetryprotected and symmetry-enriched topological phases, which are protected by them. Our results showcase the richness of the landscape of quasiperiodically-driven phases, and excitingly opens up new directions in the rapidly developing field of nonequilibrium quantum matter.
Before we proceed, let us note that the study of the dynamics of quantum systems under quasiperiodic driving has a venerable history, encompassing diverse applications from experiments in chemistry and physics, to the basic structure of first-order differential equations [38][39][40][41][42][43][44][45][46][47][48][49][50][51]. Interesting dynamical behavior related to topology in few-body or non-interacting scenarios have also been reported [52][53][54]. Many-body quasiperiodically-driven quantum systems with interactions have received increasing attention comparatively recently [55][56][57]. Our approach is distinguished from previous studies in that we explicitly establish the stability of the nonequilibrium phases we discuss, by rigorously providing a bound on their lifetimes. We additionally demonstrate the robustness of their universal properties against small changes in the driving protocol, which justifies their characterization as 'phases of matter'.
Remark on notation: the term 'time quasi-crystal' (TQC) has been used for systems that show a quasiperiodic response, arising from either a quasiperiodic [55] or a periodic drive [58,59]. In this manuscript, we will restrict use of the term time quasi-crystal to the first sense, where a quasiperiodic drive gives rise to a response with a different quasiperiodic pattern (see Sec. IV).  In this paper, we show the existence of long-lived nonequilibrium phases of matter protected by multiple time-translational symmetries (TTSes). The natural scenario to consider is quasiperiodically-driven systems, which are defined as follows. Consider an at-least piecewise continuous Hamiltonian H(θ) parametrized by the m-dimensional "standard" torus T m θ = (θ 1 , · · · , θ m ), which is 2π-periodic in each angle θ i . Additionally, let us pick a vector of rationally independent frequencies ω = (ω 1 , · · · , ω m ), so that n · ω = 0 for any non-zero integer vector n ∈ Z m . When m = 2, it suffices to choose the ratio ω 2 /ω 1 to be an irrational number, such as the golden ratio ω 2 /ω 1 = ϕ = (1 + √ 5)/2. The dynamics of the system under the time-dependent Hamiltonian is time-quasiperiodic for m ≥ 2, and constitutes a quasiperiodic drive. Here, t is the physical 'single' time and θ 0 ∈ T m are some fixed initial angles (see Fig. 1). Note that this class of drives encompasses Floquet driving as the case m = 1, which will enable us to make a direct comparison of our analysis to previous work. At first glance, it is puzzling how one could obtain phases protected by TTSes in quasiperiodically-driven systems. After all, the incommensurate nature of the driving frequencies imply that the Hamiltonian Eq. (1) has not even a single time-translational symmetry (there is no period T = 0 for which H(t) = H(t + T ) for all t), let alone multiple TTSes. However, H(t) derives from an underlying H(θ) on T m through Eq. (1), which is symmetric under translations θ → θ+τ . Here the translation vectors belong to a lattice τ ∈ L = 2πZ m generated by m independent symmetries. It is thus conceivable that these symmetries have meaningful and nontrivial implications for the single-time system. Since for m = 1 this is nothing but the time-translation symmetry of a Floquet system, we will still refer to these symmetries as "timetranslation symmetries" for quasiperiodically-driven systems (m ≥ 2).
In what follows, we will make the above statements precise by interpreting the consequences of the multiple TTSes in a certain class of physical systems. Specifically, we consider quantum many-body systems defined on lattices in arbitrary spatial dimensions with bounded local Hilbert spaces (i.e. spins or fermions), and which respect a sense of locality -the interaction strength decays sufficiently fast with distance. In particular, we allow for a Hamiltonian having interactions with amplitude that are at least exponentially-decaying with distance (termed 'quasilocal'), see Sec. VII. We will show that: 1. For such strongly interacting many-body systems under some non fine-tuned quasiperiodic driving θ on the trajectory θ(t) = ωt + θ0 (blue arrows). As the frequency vector ω is a set of rationally independent frequencies, the trajectory never returns to itself, but covers the torus uniformly. Shown here is the case m = 2, with tan α = ω2/ω1 irrational.
conditions, the multiple time-translation symmetries of H(θ) give rise to an actual symmetry of the effective time-independent Hamiltonian that describes the dynamics in a long-lived preheating regime. This enables the existence of novel nonequilibrium phases of matter protected by these symmetries.
2. The classification of quasiperiodically-driven manybody phases of matter -both spontaneous symmetry breaking and topological phases -is the same as the classification of equilibrium phases with a symmetry group extended by Z ×m . This is a direct generalization of the m = 1 Floquet results of Ref. [60].

B. Long lifetimes in quasiperiodically-driven systems
Owing to the lack of energy conservation, a generic ergodic interacting driven system is expected to heat to a featureless infinite-temperature state, where symmetries act trivially and a discussion of phases of matter is moot. Therefore, before we can even discuss new phases realizable with multiple TTSes, we must establish that there exist suitable quasiperiodic driving conditions where such deleterious heating is controlled, at least for some parameterically long time.
In Floquet systems (m = 1), the conditions needed to achieve this are relatively mild. Energy absorption or emission between the system and the drive can only take place in integer multiples of the driving frequency, i.e. ∆E = nω for integer n. By suppressing resonances between energy eigenstates connected by such discrete energy levels, heating can be slow. Both Floquet-MBL and Floquet prethermalization involve suppressing heating in a such a way, though through different physical mechanisms. The former uses strong disorder to directly curtail the probability of local resonances [30,36,37]. The latter entails driving at such high frequencies compared to local energy scales J ω that the system can only absorb the large drive quanta nω by performing a multiple-spin rearrangement. This rate is heavily suppressed giving rise to a long heating timescale t * ∼ e const.ω/J [31].
In quasiperiodically-driven systems (m ≥ 2), by contrast, energy absorption or emission occurs in units of ∆E n = ω · n for any integer vector n ∈ Z m . As the frequencies ω = (ω 1 , · · · , ω m ) are rationally independent, the set of all possible such quanta ∆E n is dense on the real line. Superficially, it seems impossible to avoid immediate heating.
A more careful consideration shows however that this is not necessarily an insurmountable problem. Consider, as an example, a static Hamiltonian H 0 weakly driven by a quasiperiodic perturbation V (t) = V (ωt + θ 0 ). Expanding V (θ) in a Fourier series one sees that V n induces transitions between energy levels of H 0 separated by energies ∆E n in linear response theory. The set of all ∆E n with |n| below some cutoff becomes ever more closely spaced as the cutoff increases. Nonetheless, as long as the amplitude V n decays fast enough with |n| as compared to this spacing, resonances will not proliferate at large |n|. This observation suggests a restriction of V (θ) to be 'sufficiently' smooth on the torus, which in turn translates to smooth drives in time.
Resonances arising from small |n| processes also need to be suppressed -but this should be achievable with the same mechanisms as in the Floquet case (strong disorder or high-frequency driving). One of the key aspects of our present work is to make this plausible stability statement concrete. Specifically, we will consider the case corresponding to a quasiperiodically-driven many-body Hamiltonian H(ωt+ θ 0 ) under two conditions. First, that it is smooth in time, in the sense that the amplitude of the Fourier coefficients H n of the driving Hamiltonian decay as H n e −const.|n| at large |n|. Second, that the driving frequencies ω i are large compared to any local energy scales J of the Hamiltonian. Under these conditions, we will show that the heating time t * is rigorously bounded for any > 0, and for all except a set of measure zero choices of frequency vectors ω by Here ω = |ω| is the norm of the driving frequency, C, C are dimensionless numbers depending on the numbertheoretic properties of the irrational ratios ω i /ω j , and m is the number of incommensurate frequencies. While the stretched exponential form of t * can be intuitively understood within linear response theory (Sec. VI), we will prove the bound Eq. (3) using a recursive construction beyond linear response (Sec. VII). Note that this result is analogous to the exponentially-slow heating bound known for Floquet systems; however in those cases there were no restrictions on the smoothness of the drive.

C. Description of preheating dynamics
The existence of a long timescale Eq. (3) implies a long-lived preheating regime, and opens up the possibility to define phases of matter in this time interval. How can we concretely describe the dynamics, and eventually characterize phases, in the preheating regime?
Let us briefly recall the Floquet scenario (m = 1), where generally the dynamics in a preheating regime is approximately governed by an effective, static, quasilocal Hamiltonian D. More precisely, the time-evolution operator where P (t) = P (t+T ) is a unitary change of frame which is periodic in time, and the approximate equality reflects an omission of small local terms in the Hamiltonian that do not affect the dynamics up to the heating time t * [32][33][34][35]. A quantum state's dynamical evolution can therefore be understood as comprised of two parts: time-evolution generated by the static Hamiltonian D, and an additional 'micromotion' governed by P (t). If we only consider the state of the system at stroboscopic times, i.e. integer multiples of the driving period, then the dynamics is generated by the Floquet operator U F := U (T ). This operator can be written [if we choose P (0) = P (T ) = I] as U F ≈ e −iDT , in which case the study of the dynamical evolution of the system up to time t * entails studying the eigenstates of the Hamiltonian D.
We strongly emphasize here the importance of D being quasilocal. In fact, Floquet-Bloch's theorem asserts that the decomposition Eq. (4) exists with an exact equality if D is replaced by the Floquet Hamiltonian H F . However, H F will be highly non-local in a generic ergodic manybody system -it must after all describe the eventual heating to an infinite-temperature state -and therefore is not very insightful to use when studying the preheating regime, as compared to D.
One scenario where a long heating time emerges and the approximate decomposition (4) holds is in the limit of high-frequency driving, in which case t * e const.ω/J . The small ratio of the local energy scale to the driving frequency J/ω naturally enables schemes for an orderby-order expansion of D and P . For example, the commonly-used "Floquet-Magnus" expansion [5,61] with P (0) = P (T ) = I gives at lowest orders and is generally an asymptotic series. Refs. [32,33] showed that if truncated at some optimal order, Eq. (4) is satisfied with an error ∼ 1/t * = O(e −const.ω/J ).
Ref. [34] proved a similar rigorous statement for an effective Hamiltonian D obtained through a different expansion, using techniques that are not exactly a perturbative expansion in J/ω. We now return to quasiperiodically-driven systems (m ≥ 2). For the smooth high-frequency drives that we consider, we will prove that a similar decomposition of the unitary time evolution operator U (t) as Eq. (4) exists, A quantum state's dynamics is again effectively comprised of time-evolution by some static Hamiltonian D and a micromotion given by some unitary timequasiperiodic change of frame P (t) := P (ωt + θ 0 ) with underlying P (θ) smooth on the torus. In Sec. VII, we will demonstrate how to construct the effective Hamiltonian D and unitary P through an iterative renormalization procedure of the driving Hamiltonian H(t) that can be understood as a generalization of the methods of Ref. [34], as well as bound the optimal order to which the procedure should be carried out. This gives rise to an optimal D such that the description Eq. (6) is valid at least for times t t * , with t * satisfying Eq. (3). While it seems natural to assume that the decomposition Eq. (4) carries over from the Floquet case to quasiperiodically-driven systems, this is far from obvious. It is known rigorously that the decomposition Eq. (6) with an exact equality (i.e. Floquet-Bloch's theorem) is not guaranteed in general quasiperiodic systems, there being obstructions to defining a generalized Floquet Hamiltonian [40]. To understand why one does expect (6) to hold in the quasiperiodically-driven case with conditions given in Sec. II B, observe that the highfrequency assumption suggests that an expansion analogous to the Floquet-Magnus expansion Eq. (5) can be written down. Representing H(θ) as a Fourier series H(θ) = n H n e in·θ and assuming the form Eq. (6) with P (t) quasiperiodic, one can perform a formal expansion in powers of 1/ω of the effective Hamiltonian D = D (0) + D (1) + · · · as well as the unitary P (θ) = exp(Ω (1) (θ) + Ω (2) (θ) + · · · ), whose leading order terms read (see Appendix A) Here Ω (1) (θ) = − n =0 H n e in·ω /(n·ω). However, while relatively simple to construct, even the low order terms in Eq. (7) already signal a difficulty not present in the Floquet case: The denominator n · ω can be arbitrarily small, leading to possible divergences and bringing into question the validity of the expansion. This is precisely the manifestation of the denseness of resonances discussed in Sec. II B. As before, we observe that this issue can potentially be circumvented if the size of the Fourier coefficients H n decay sufficiently rapidly with |n|, such as in the case of smooth driving, so that the small denominators are suppressed. Note also that while similar to Eq. (5), the expansion of Eq. (7) does not reduce to it upon setting m = 1, as we explain in Appendix A. This point which will be important in the discussion of emergent symmetries. A central contribution of our paper is to show how imposing the smoothness conditions on the drive indeed leads to a meaningful high-frequency expansion which can be used to construct an effective static Hamiltonian. The expansion we develop, which is different from that of Eq. (7), is given in Sec. VII.

D. Equilibration and steady states in the preheating regime
Let us now describe the kind of dynamics and 'steady states' one can expect in the preheating regime, given an effective, static, quasilocal Hamiltonian D via Eq. (6).
If D is generic and non-integrable, one expects dynamics to lead to thermalization with respect to D, when viewed in the time-dependent frame defined by P (t). Explicitly, the system locally approaches an equilibrium distribution ρ β ∝ e −βD , where the inverse temperature β depends on the energy of the initial state as measured by D. This happens provided local relaxation timescales t r are much less then the heating timescale t r ∼ J t * , which always occurs for high-frequency drives. As the system is expected to eventually thermalize to an infinitetemperature state, ρ ∝ I for t > t * , due to the corrections in Eq. (6) that cannot eventually be neglected, one refers to equilibriation to ρ β in the preheating regime as prethermalization. In particular, we can then talk about 'prethermal quasiperiodically-driven phases of matter' in this steady state. Of course, one must remember to include the effects of the time-quasiperiodic unitary P (t) upon moving back to the laboratory frame. However, as P (t) when constructed in the high-frequency limit is perturbatively close to the identity, this simply endows the steady state ρ β with additional micromotion of amplitude ∼ J/ω; see Fig. 2.
We can also consider the case where D is non-ergodic, such as when it is highly disordered leading to many-body localization (MBL). In this case, there is a complete set of quasilocal integrals of motion τ z i ('l-bits') satisfying [τ z i , D] = 0. The system evolving under D, when viewed in the rotating frame P (t), will not thermalize, but in- stead exhibit MBL phenomenology. This includes logarithmic entanglement growth, initial state memory and localization-protected quantum order [62][63][64][65].
In the laboratory frame, the dynamics is rather more interesting: Owing to the rotating change of frame, one finds that the l-bits are never integrals of motion. Instead, we find that the "reverse Heisenberg" evolutioñ τ z i (t) R = U (t)τ z i U (t) † depends on time quasiperiodically, i.e.τ z i (t) R =τ z i (ωt + θ 0 ), where we have defined the l-bits in the laboratory frame asτ z i := P (0)τ z i P (0) † , which is localized near site i. Indeed, we haveτ z i (t) ≈ P (t)e −itD τ z i e itD P (t) † = P (t)τ z i P (t) † . Since ωt + θ 0 can get arbitrarily close to θ 0 on the torus for certain times t, this implies thatτ z i (t) R can have arbitarily close recurrences such that τ z i (t) R ≈ τ z i . This also implies that the usual forward Heisenberg evolutionτ z i (t) = U † (t)τ z i U (t) has arbitrarily close recurrences at the same times, though the latter is not necessarily quasiperiodic. Thus, we see thatτ z i (t) R andτ z i (t) never spread out too far from site i but also never precisely return tõ τ z i . This behavior is different from Floquet-MBL where [τ z i , U F ] = 0 so l-bits are truly integrals of motion at stroboscopic times. Nonetheless, the dynamics is still far more restricted than in an ergodic system where any local operator spreads ballistically or diffusively. What we have described is thus a new kind of dynamical localization that can be dubbed "quasiperiodically-driven MBL", albeit the fact that we have only shown that it is stable until t * bounded by Eq. (3). Proving whether the quasiperiodically-driven MBL is stable to beyond this time, perhaps even forever, remains an interesting direction for future work.
In this paper, we will consider quasiperiodically-driven phases of matter realizable in one or the other of the scenarios described above: prethermalization or (stretchedexponentially long-lived) MBL.

III. EMERGENT SYMMETRIES PROTECTED BY MULTIPLE TIME TRANSLATION SYMMETRIES
Having motivated quasiperiodically-driven systems and outlined their dynamics in suitable regimes, we now analyze what kind of new phases of matter can arise in these systems. As a first step, let us consider the scenario where a direct high-frequency drive is applied to a system. This procedure is often referred to as 'Floquet engineering', as the drive is used to modify and control interactions of an underlying Hamiltonian. Indeed, the ground states of the effective static Hamiltonian D that is generated in a high-frequency expansion can be different from those of the original undriven Hamiltonian [5,[66][67][68][69][70].
However, from a phases of matter point of view, a direct high-frequency drive will not yield fundamentally new long-time collective behavior that is not already reproducible in some -possibly complicated -static system at equilibrium. This is because a quantum state's evolution is effectively governed entirely by D and never has any significant nontrivial micromotion during its time evolution. Precisely, this stems from the fact that the unitary frame transformation P (t) in the description of the time evolution operator Eq. (6) is perturbatively close to the identity. To uncover novel phases, especially those that are inherently out-of-equilibrium, we need to go beyond.

A. Frame-twisted high-frequency limit in Floquet systems
For Floquet systems (m = 1), Ref. [12] provided general periodic driving conditions which do give rise to fundamentally new non-equilibrium phases. We briefly review these here.
The main idea is to consider periodically-driven Hamiltonians that approach the high-frequency limit, but only when viewed in a certain rotating frame, a so-called 'frame-twisted high-frequency' limit. Consider a manybody driven system with Hamiltonian of the form H(t) = H 0 (t) + V (t). Here H 0 (t) = H 0 (t + T ) is a sum of quasilocal terms, with associated time-evolution opera- Evolution under H 0 is taken to have the special property U 0 (N T ) = U 0 (T ) N = X N = I, for some positive integer N > 1 and an operator X, which is not itself the identity. The term V (t) = V (t + T ) describes interactions assumed to have local energy scale J ω/N . Since X is not perturbatively accessible from the identity, a strong drive is required to realize this evolution -the local energy scale of H 0 is ∼ ω/N , and thus increases with ω. One therefore cannot naively apply the high-frequency expansion Eq. (5).
In the rotating frame defined by the interacting picture of H 0 (t), time-evolution is governed by the interaction Hamiltonian This is a quasilocal Hamiltonian which is still timeperiodic, albeit with period N T . Since J ω/N , a high-frequency expansion can be meaningfully applied to it, and one can see how there is a long-lived preheating regime in this frame of reference.
However, the frame-twisted high frequency limit imposes a stronger condition than just long-lived preheating. Specifically, as shown in Ref. [12], in the laboratory frame the Floquet unitary U F := U (T ) takes on the special structure where D is a quasilocal Hamiltonian which additionally satisfies [D, X] = 0 identically. Here V is a small time-independent quasilocal unitary. The approximate equality reflects an omission of small time-dependent local terms whose effects only become relevant after times t * ∼ O(e const.ω/J ). The physical statement is that a system periodically driven under these conditions always has an emergent Z N symmetry generated by X. One can add small, potentially time-dependent perturbations to H(t) as long as they respect the time-periodic nature of the drive, and the structure of Eq. (9) will be unchanged. The emergent symmetry is therefore robust and underpins the stability of inherently nonequilibrium Floquet phases of matter that can now emerge. From Eq. (9), one sees that when observed at times t that are integer multiples of N T (i.e. at times which are stroboscopic with respect to the longer period), the system, in the time-independent frame described by V, settles into an equilibrium state of D distinguished by the emergent Z N symmetry. However, when viewed after every time interval T (the original period of the driving Hamiltonian), due to the action of the symmetry operator X, the state of the system can transform nontrivially. This happens, for example, if the system spontaneously breaks the symmetry. This additional periodic action is precisely what makes the longtime collective behavior of this system inherently out-ofequilibrium, and the robustness of the phenomenology justifies the terminology of them being called fundamentally nonequilibrium phases of matter. We reiterate that this remarkable result is a consequence of the discrete TTS of the Floquet drive and guaranteed to happen with no additional symmetry requirements.

B. Twisted time-translation symmetries and emergent symmetries
In quasiperiodically-driven (m ≥ 2) systems, it is natural to look for an expression of the form Eq. (9). Since there is no single time-translation symmetry, however, there is no analog of the Floquet operator U F and the construction of Ref. [12] does not carry over.
Our key observation is that one can rederive the results of Ref. [12] for Floquet systems in a considerably simpler way that does accord an extension to quasiperiodicallydriven systems. This relies on realizing that the interaction Hamiltonian H int (t) described in the previous section possesses a symmetry as a consequence of the TTS of the original laboratory frame Hamiltonian. We will refer to this symmetry of H int (t) as a 'twisted time-translation symmetry'.
Precisely, in the Floquet setting, we say that a timeperiodic operator O(t) = O(t + T ) has a twisted TTS, if there is an integer N > 1 and a unitary operator g satisfying g N = I, such that for all t whereT = T /N . In terms of the Fourier modes O n in O(t) = n O n e inωt , the twisted-TTS states that gO n g † = e 2πin/N O n . Note that H int (t) in Eq. (8) has a twisted-TTS with unitary g = X † , provided we rescale time t → t/N so that the periodicity of H int (t), originally N T , becomes T .
To gain some intuition as to why this concept is useful, suppose that we had a Hamiltonian H(t) with a twisted-TTS and we were to construct the effective Hamiltonian D from the high frequency expansion given by Eq. (7) with m = 1. One immediately sees from the action of twisted-TTS in Fourier space that D constructed this way commutes with g to all orders. Additionally, it can be shown that the change of frame P (t) will also inherit the same twisted-TTS as H(t), i.e. P (t +T ) = gP (t)g † . Note that the usual Floquet-Magnus expansion Eq. (5) will not give this result, see Appendix A. The physical conclusion is that dynamics of a driven Hamiltonian with a twisted TTS can always be viewed in some frame as effectively governed by a time-independent Hamiltonian D with an emergent Z N internal symmetry, generated by g.
Applying these considerations to the rotating frame Hamiltonian H int (t) in Eq. (8) to construct the effective Hamiltonian D using Eq. (7), ones recovers the statement Eq. (9) of Ref. [12] in a transparent fashion: Since U F = U 0 (T )U int (T ) = XT e −i T 0 Hint(t )dt , we can write it as with [D, X] = 0, and where we have used XP (T ) = P (0)X. The twisted-TTS concept immediately generalizes to quasiperiodically-driven systems (m ≥ 2). Recall that a time-quasiperiodic operator O(t) = O(ωt + θ 0 ) derives from an operator O(θ) living in a higher-dimensional space, where O(θ + τ ) = O(θ) for any τ ∈ L = 2πZ m . Suppose, there is additionally some finite translation vectorτ and a unitary operator gτ satisfying g Ñ τ = I for some integer N > 1, such that We then say that O(θ) has a gτ -twisted time-translation symmetry. Note that as O(θ + Nτ ) = O(θ), Nτ ∈ L.
In terms of Fourier modes O n , the twisted-TTS acts as gτ O n g † τ = e in·τ O n . Furthermore, since O(θ) is defined on a torus with dimension m ≥ 2, it can have multiple independent twisted-TTSes corresponding to different translation vectorsτ and unitary operators gτ .
As in the Floquet case, for a Hamiltonian with twisted-TTSes, the effective Hamiltonian D constructed through Eq. (7) manifestly commutes with gτ . Also, P (θ) will inherit the same twisted-TTSes P (θ +τ ) = gτ P (θ)g † τ . These properties will hold for other high-frequency expansions as well, like the one we develop in Sec. VII.
C. Frame-twisted high-frequency limit in quasiperiodically-driven systems Although Eq. (11) seems like an obscure condition, analogous to the Floquet case, twisted TTSes can arise naturally in a frame-twisted high-frequency limit of a drive that does not need to satisfy additional symmetry constraints beyond time-quasiperiodicity itself. Indeed, we will see how m twisted-TTSes can manifest from m 'untwisted' time-translations of the original driving Hamiltonian H(θ), when viewed in a suitable frame of reference.
As discussed in Sec. II A, the many-body Hamiltonian of the system H(t) = H(ωt + θ 0 ), derives from a Hamiltonian H(θ) on the standard torus T m . Assume now this has the form with H 0 (θ + τ ) = H 0 (θ) and V (θ + τ ) = V (θ), for any τ ∈ L = 2πZ m . Let Γ 1 , . . . , Γ r be a set of r mutually commuting operators, each of which is a sum of quasilocal terms and has integer eigenvalues. We take Unless otherwise stated, the index summation convention is implied. Here f i (θ) are real-valued functions satisfying where Q ij is a dimensionless r × m matrix with rational entries. The interaction  (12), has periodicity on the standard torus (grey shaded area). The interaction Hamiltonian Hint(θ) generally has different periodicity. Shown is an example where its periodicity is on a sheared torus (red shaded area), given by the translation vectors τ 1 = 2π(2, 1) and τ 2 = 2π(1, 2) which defines a lattice L (red points). (b) There exists an invertible linear transformation R mapping θ → Rθ so that (τ 1 , τ 2 ) → (Rτ 1 , Rτ 2 ) = (τ 1, τ 2). The Hamiltonian H int (θ) := Hint(R −1 θ) has periodicity on the standard torus (red shaded area), and additionally has gτ -twisted TTSes is itself quasiperiodic, although defined on a larger torus than H 0 (θ). Specifically, U 0 (t) = U 0 (ωt + θ 0 ) for some function U 0 (θ), which satisfies Here L is a sublattice of the original lattice L defined by τ ∈ L if and only if τ ∈ L and e iQij τ j = 1. A simple example is where m = r and Q = (1/N )I m for some integer N > 1; in that case L = 2πN Z m = N L, so that the basic original cell is enlarged by N in each direction. The explicit form of the time-evolution operator is . We emphasize that Eq. (14) holds only due to our special form of H 0 (t) from Eq. (13). In general, one does not expect U 0 (t) to have any quasiperiodicity, even if H 0 (t) does.
We are now in a position to see how twisted TTSes emerge. By transforming into the interaction picture of H 0 (t), the interaction Hamiltonian is and has local energy scale J. Furthermore, it derives from a Hamiltonian H int (θ) = U † 0 (θ)V (θ)U 0 (θ), which has periodicity in the larger unit cell H int (θ+τ ) = H int (θ) with τ ∈ L . Choosing Q to have rational entries ensures that the new unit cell is still finite (see Appendix C).
Since the Γ i commute, Eq.
for τ ∈ L. Together with V (θ) = V (θ + τ ), this yields This is almost the twisted-TTS condition Eq. (11), except the periodicity of H int (θ) is not on the standard torus T m . However, there is always an invertible linear transformation R on θ-space ( Fig. 3) so that the Hamiltonian H int (θ) := H int (R −1 θ) is periodic on the standard torus and Eq. (11) then holds exactly for H int (θ). We will generally not need to invoke the transformation R, except when we discuss estimates of heating times and our construction of effective Hamiltonians in the preheating regime, which take as input high-frequency Hamiltonians that are periodic on the standard torus. In those cases, we should remember that the coordinate transformation means the the frequency vector is also rescaled ω → Rω when we consider the single-time evolution, since H int (ωt + θ 0 ) = H int (Rωt + Rθ 0 ). This implies that, the dynamics governed by the effective Hamiltonian D constructed from H int (θ) in a high-frequency expansion assuming J R ij ω j will last for a long time t * ∼ O(e const.×(|Rω|/J) 1/(m+ ) ); see Sec. VII. Then, for times less than t * , the evolution operator in the interaction frame can be written as where P (θ) is periodic with respect to translations τ ∈ L .
According to the discussion of Sec. III B, the effective Hamiltonian D in the preheating regime will have emergent symmetries arising from the twisted TTSes, Moreover, P (θ +τ ) = g τ P (θ)g † τ for τ ∈ L. Analogous to the Floquet case, these symmetry properties of D are robust to small, potentially time-dependent perturbations to the driving protocol, as long as they respect the timequasiperiodicity of the system.
Although Eq. (20) holds for each translation vector τ ∈ L, not every such τ corresponds to a different operator g τ . In fact, g τ = I if and only if τ ∈ L . The unitary operators g τ therefore belong to a finite Abelian group of emergent symmetries, G ∼ = L/L . In the simple case where Q = (1/N )I m , we find that m = r, G = Z ×m N . As a slightly less trivial example, consider r = m = 2 and . This results in a lattice L as seen in Fig. 3, and G = Z 3 . We refer the reader to Appendix C where we show, given some rational matrix Q, how to extract to compute G from the Smith decomposition of Q. Since G is a finite Abelian group, it is always of the form G = Z q1 × Z q2 × ... [71]. Collecting all the ingredients discussed above, we are now in the position to realize novel, inherently out-ofequilibrium phases of matter. Time evolution in the laboratory frame for t < t * , under the driving scenarios outlined in this section, is governed by the evolution operator U (t) = U 0 (t)U int (t), which can be written as Here that has periodicity on the standard torus T m , that is, for translations τ ∈ L. Another way to state Eq. (21) is that in the rotating frame defined by V(t)U 0 (t), time evolution of the system is simply governed effectively by the static Hamiltonian D. However, if one goes back to the laboratory frame, this frame transformation -being not perturbative close to the identity -can endow the state of the system (in particular, the steady state of D) with large, structured time-quasiperiodic micromotion, giving rise to a panoply of different long-time dynamical collective behaviors whose dynamical signatures are robust and universal.
In the next two sections, we will illustrate the physical implications of our results with two examples of such phases: time quasi-crystals and dynamic quasiperiodic topological phases. We will return to the important task of formalizing the preceding discussions on dynamics as well as explicitly constructing D, in the later sections.

IV. PHASE I: DISCRETE TIME QUASI-CRYSTALS
A discrete time quasi-crystal (DTQC) is a phase which spontaneously breaks some or all of the time-translation symmetries of a quasiperiodic drive [55]. It is characterized by a dynamical response of physical observables, which display stable long-time oscillations with a time-quasiperiodicity that is different from the timequasiperiodicity of the original driving Hamiltonian H(t). This can be diagnosed by computing the power spectra of local observables, which will exhibit robust peaks at frequencies which are shifted from the base frequencies by a fractional amount. Since there are several (m ≥ 2) independent time-translation symmetries, there are a multitude of ways that these symmetries can be spontaneously broken, leading to a variety of patterns and associated DTQC phases. The DTQC generalizes the discrete time crystal, a phase which spontaneously breaks the single time-translation symmetry of Floquet systems [9,11]. Note that the concept of a DTQC as well as some aspects of its phenomenology have been proposed in [55], which numerically observed a DTQC-like signal in a quasiperi-odic step-drive with disorder; albeit with a slow logarithmic decay of the envelope in time. Our present work shows how the logarithmic decay can be avoided (even without disorder) through smooth driving, hence rigorously proving the stability of the phase up to the long heating time t * . In addition, we provide drives that generalize to a large class of symmetry breaking patterns.
To understand DTQC phases, consider the time evolution of a quantum state |ψ 0 with a quasiperiodic Hamiltonian of the type discussed in Sec. III C. Then one can take a frame-twisted high-frequency limit, so that for times t < t * , the time-evolution operator can be decomposed as in Eq. (21). Recall that the effective time-independent Hamiltonian D in the preheating regime possesses multiple unitary symmetries g τ where τ ∈ L = 2πZ m , which belong to some finite Abelian group G. Now let us consider times where the system has prethermalized, so that in the rotating frame V(t)U 0 (t) the state is locally described by a thermal state ρ β ; see the discussion in Sec. II D. In the laboratory frame, the state of the system when probed by local observables is We see that ρ(t) is time-quasiperiodic, since ρ(θ) at least satisfies ρ(θ + τ ) = ρ(θ) for τ ∈ L . But does ρ(θ) have the periodicity of the original drive H(θ) of Eq. (12), characterized by the lattice L = 2πZ m ? This turns out to depend on whether or not ρ β is symmetric under the emergent symmetries g τ . To see this explicitly, we can write where we use that V(θ + τ ) = V(θ) and U 0 (θ + τ ) = U 0 (θ)g † τ . Therefore ρ(θ +τ ) = ρ(θ) for τ ∈ L if and only if g τ ρ β g † τ = ρ β . In other words, if ρ β is a state that preserves all the symmetries g τ of the effective Hamiltonian D, then ρ(θ) preserves all m multiple time translation symmtries of the driving Hamiltonian H(θ).
If on the other hand ρ β is not invariant under g τ , g τ ρ β g † τ = ρ β , then the emergent symmetry g τ is said to be spontaneously broken in the thermal state of D. This can, of course, happen for multiple g τ at the same time. From Eq. (23), ρ(θ) will then have a periodicity different from the original Hamiltonian H(θ), and consequently ρ(t) = ρ(ωt + θ 0 ) will have a different timequasiperiodicity than the driving Hamiltoinian H(t) = H(ωt + θ 0 ) -the hallmark of a DTQC phase. The precise connection between the spontaneous breaking of g τ and of TTS reflects the fact that g τ was a manifestation of the TTS in the first place.

A. Observable consequences
The spontaneously broken TTSes in quasiperiodicallydriven systems manifest themselves most clearly through periodicity changes in the m-dimensional θ-space. The interplay of multiple time translation symmetries gives a large variation in the number of different symmetry breaking patterns and associated DTQC phases. However, there will also be measurable signatures of these patterns in terms of the dependence of the system on physical time t, for example, in the Fourier spectrum (or power spectrum) of local observables (see also [55]). These are analogous to probing quasicrystalline structures in space through their diffraction patterns [72,73].
Consider the regime described above, where the state of the system is described by Eq. (23). Then, the expectation of a local observable o(t) := ô(t) can be written as o(t) = o(ωt+θ 0 ), where o(θ) := Tr(ôρ(θ)) has periodicity that depends on the periodicity of ρ(θ). Let L SSB be the largest lattice such that ρ(θ + τ ) = ρ(θ) for all τ ∈ L SSB , which describes the symmetry-breaking pattern in θ space. Then we can expand o(θ) as a Fourier series where the sum is over the reciprocal lattice vectors α ∈ L * SSB , which are the vectors α satisfying e iα·τ = 1 for all τ ∈ L SSB . (Note that for smooth driving, the Fourier coefficient o α will decay exponentially with |α|.) Consequently, the power spectrum of o(t) has peaks at frequencies Ω α = α·ω. Of course, not every peak will necessarily appear for any choice of observablê o, since it is possible that o α = 0 for some α's. Note also that any finite observation time will broaden the peaks of the spectrum; for the interacting many-body phases we consider, the heating time t * sets a fundamental limit on resolution. Now, in a DTQC, since the symmetry lattice L spontaneously breaks to the proper sublattice L SSB , the reciprocal lattice L * SSB is also a proper superlattice of L * = Z m . This implies that some of the frequencies Ω α = α · ω for α ∈ L * SSB are not derivable from integer linear combinations of the base driving frequencies (ω 1 , · · · , ω m ), i.e. they do not correspond to the base harmonics Ω n = n · ω for n ∈ Z m . This, is the dynamical signature of the spontaneous breaking of the time translation symmetries. Of course, the frequencies associated with original drive harmonics Ω n = n · ω are dense on the real line, so they can lie arbitrarily close to those frequencies Ω α reflecting the symmetry-breaking. However, the peaks in the power spectrum at frequencies Ω α will nevertheless be well resolved from those at Ω n . This is because the weights of peaks at frequencies Ω n approaching Ω α become ever more strongly suppressed as a consequence of the smoothness of the drive, as they involve very large |n|. Therefore, the presence of well-defined peaks at frequencies Ω α constitutes a sharp dynamical signature of the DTQC phase. Additionally, an important signature of the DTQC phase is that the location of these peaks Ω α in the power spectrum reflecting the spontaneous symmetry-breaking, is robust against small perturbations to the driving protocol, such as in changing f (θ) → f (θ) + (θ) for small, smooth (θ), or by adding small time-quasiperiodic terms to the Hamiltonian V .
In the following subsection as well as Appendix D, we provide examples of Hamiltonians that exhibit DTQC phases.
B. Example Hamiltonian: Z2 DTQC Consider a system of spin-1/2 degrees of freedom on a lattice, evolving with the m = 2 time-quasiperiodic Hamiltonian Here σ x , σ y , σ z are the standard Pauli matrices, and we will choose the couplings J ij to be ferromagnetic so that the Ising Hamiltonian i,j J ij σ z i σ z j has an ordered phase at finite temperature. This could be achieved using nearest-neighbor couplings in two or greater dimensions, or with power-law decaying interactions in one dimension; see Sec. VIII A for a discussion of the dynamical consequence of power-law interactions.
The driving profile f (θ) is given by with Q = 1/2 1/2 and a function parameterized by an integer N , see Fig. 4(a,c) for an illustration of the driving protocol for N = 20. The Hamiltonian Eq. (26) falls into the class of Hamiltonians described in Sec. III C. It comprises of two terms: first, V describes pairwise Ising interactions between spins with amplitude J ij , as well a longitudinal field in the z-direction with strength h. The couplings are assumed to satisfy J ij , h |ω 1 − ω 2 |/2. Second, H 0 (ωt) describes a quasiperiodic drive on the system in the xdirection, with frequency vector ω = (ω 1 , ω 2 ). In an experimental platform, this can be implemented, for example, by external pulses using lasers or microwaves. For concreteness, we will henceforth take ω 1 = 1 and ω 2 = (1 + √ 5)/2, the golden ratio. Due to the finite number of Fourier components, f (θ) is a smooth function on the torus. The function 1 + F N (θ) is related to a so-called 'Fejer kernel' and can be intuitively thought of as a smoothed-out Dirac delta comb 2π m∈Z [δ(θ + 2πm) − 1], with N governing the degree of smoothness.
Under these premises, let us now take the frametwisted high-frequency limit. Going into the interaction frame of H 0 (θ), we can compute the time-quasiperiodic interaction Hamiltonian H int (θ) = U † 0 (θ)V U 0 (θ). The periodicity of H int (θ) is on the lattice L , generated by the translation vectors τ 1 = 2π(1, 1) and τ 2 = 2π(−1, 1). Here L is a sublattice of the original lattice L = 2πZ m of symmetry vectors of H(θ). One also sees that H int (θ) possesses a single nontrivial twisted-TTS corresponding to a translation byτ = 2π(1, 0) orτ = 2π(0, 1): We next construct the effective time-independent Hamiltonian D from the H int in a high-frequency expansion, using our approach in Sec. VII. However, for our purpose here of understanding its steady states, it suffices to understand the leading order Hamiltonian D (0) in the high-frequency expansion: Since D and D (0) are perturbatively close by construction, the steady states of D and D (0) are in the same universality class. For our expansion in Sec. VII, the leading order term D (0) is the average of the interaction Hamiltonian We note in particular the bounds on integration, corresponding to a unit cell of L . We find that  30) is Ising-symmetric, that is, The steady states of D in the preheating regime are as follows. Provided the initial state has energy density (measured with respect to D (0) ) below some critical energy density, the system will prethermalize to a Gibbs state ρ β that spontaneously breaks the Ising symmetry g = i σ x i . This means that, the expectation value of an operator odd under the symmetry, such as σ z i , is generically not zero, i.e. Tr[σ z i ρ β ] = 0. From the discussion in Sec. IV and Sec. IV A, we see that the state ρ(θ), Eq. (22), has periodicity on the lattice L SSB = L , and the corresponding reciprocal lattice L * SSB is generated by the vectors α 1 = (1/2, 1/2) and α 2 = (1/2, −1/2). Thus, the power spectrum of a local observable will generically have peaks at the frequencies where n 1 , n 2 ∈ Z.  Which peaks dominate, however, depends on the operator measured and its symmetry properties under g. As a concrete example, suppose we were to measure the local observable σ z i (which is odd under g), whose expectation value in time can be written as Because of the periodicity of ρ(θ), the underlying s(θ) = Tr[σ z i ρ(θ)] will have periodicity in L SSB too. To leading order in inverse frequency (i.e. treating V(t) ≈ I), the signal can be analytically derived, and is given by where and derives from an operator G N (θ) parameterized on the torus which is periodic with translations τ in L.
In Fig. 4(b,d), we plot both s(t) and the underlying function s(θ) from which it is derived from, assuming Tr[σ z i ρ β ] = 0.8 and Tr[σ y i ρ β ] = 0. Fig. 4(e) illustrates the corresponding power spectrum P(Ω) of the signal s(t). We see that the spectrum contains peaks at the particular frequencies Ω = (n 1 + 1/2)ω 1 + (n 2 + 1/2)ω 2 where n 1 , n 2 ∈ Z; furthermore, those frequencies given by small values of n 1 , n 2 contribute the most. These frequencies are a subset of the ones described by Eq. (31) (in particular, they do not include the harmonics of the original driving frequencies, i.e. Ω = n 1 ω 1 + n 2 ω 2 ), but at higher orders in the inverse frequency we expect peaks to occur at all values Eq. (31).
While the Z 2 DTQC discussed here is a particularly simple case, we can construct more complicated DTQCs such as ones characterized by the spontaneously broken emergent symmetry groups G = Z 2 × Z 2 (Fig. 5) or Z 3 × Z 2 (Fig. 6) in an analogous fashion. We give the explicit form of the systems and drives to realize these in Appendix D. and Ω = ω2 (A2). Dominant peaks of the observable s(t) correspond to: Ω = ω1/2 (B1), Ω = ω2/2 (B2), Ω = 3ω1/2 (B3). Going beyond leading order in inverse frequency, an observable's power spectrum can exhibit peaks at the frequencies Ω = 1 2 n · ω for n ∈ Z 2 .

C. The MBL-DTQC
We are not restricted to realizing DTQC phases when the Hamiltonian D is thermalizing. Instead, we can consider the case where D exhibits MBL phenomenology and spontaneously breaks the emergent symmetries g τ . Then, the local integrals of motion τ z i are themselves not invariant under g τ , which is the usual sense of spontaneous symmetry breaking in MBL; see [62,[74][75][76][77]. For the model we considered in Eq. (26), we can pick strongly disordered interactions J ij , which leads to an "MBL spin glass" effective Hamiltonian with Ising symmetry [74,77]. In the laboratory frame, this gives rise to an MBL DTQC, whose properties we now describe.
As stated in Section II D, the defining property of MBL in a quasiperiodically-driven system is that there is a complete set of commuting operators τ z i (the "l-bits"), which evolve quasiperiodically under reverse Heisenberg We can define TTS to be spontaneously broken in a quasiperiodically-driven MBL system if the quasiperiodicity of τ z i (t) R is not the one of the drive but a different one. That is, writing . For discrete time crystals in periodically driven systems, a key feature was the spectral pairing between eigenstates of the Floquet operator U F . For example, in the case corresponding to a discrete time crystal with period-doubling, these eigenstates always come in cat state pairs separated in quasienergy by π. However, generalizing this concept to DTQC is subtle because of a , Ω = (ω1 +ω2)/2 (B3). Going beyond leading order in inverse frequency, an observable's power spectrum can exhibit peaks at frequencies lack of a single time evolution operator like U F . Finally, let us mention what behavior is expected for measurable observables for an MBL DTQC. We know that for MBL systems, local observables relax to a steady state, generically as a power law in time [78]. Once this steady state has been achieved, the power spectrum of the time dependence of local observables will display the same behavior as discussed in the prethermal case, by the same arguments. At shorter times before the steady state is achieved, one expects analogous behavior to what is seen for the discrete time crystal in Floquet systems: in addition to "universal" peaks in the power spectra of observable that persist to infinite times, one also sees other nonuniversal peaks that are dependent on the precise disorder realization of the system.

V. PHASE II: QUASIPERIODIC TOPOLOGICAL PHASES
In addition to spontaneous symmetry-breaking phases like the DTQC, one can also consider topological phases protected by multiple TTSes. In discussing these topological phases, we will focus on situations where the driving Hamiltonian has sufficiently strong disorder, so that the preheating Hamiltonian D is MBL (see Sec. II D). Since MBL eigenstates have properties analogous to the ground state of gapped local Hamiltonians at all energy densities [79], the topological classification of such ground states can be applied to each state of D. Static MBL phases can therefore be distinguished by the topology of their eigenstates; this is referred to as eigenstate order [74,80,81].
In driven systems, the concept of eigenstate order is enriched [60,[82][83][84][85][86] as there are additional topological features arising from the time evolution. Consider the Floquet time-evolution operator Eq. (4) in the pre-heating regime, where we temporarily treat the decomposition as exact. For an eigenstate |Ψ of D with energy eigenvalue , define a time evolution We refer to this as the micromotion of the eigenstate |Ψ . Even when the eigenstate |Ψ itself describes a topologically trivial static phase, its corresponding micromotion |Ψ(t) could still be nontrivial. More precisely, let Ω d denote the space of all possible gapped ground states of quasilocal Hamiltonians in d spatial dimensions (we mod out by global phase factors in the definition of points in Ω d ). In Floquet systems, the micromotion is periodic and defines a loop in Ω d -we say the micromotion is nontrivial if this loop is not contractible to a point. We can generalize the notion of topological micromotion Eq. (34) to the case of quasiperiodically-driven systems, using the decomposition Eq. (6). The micromotion |Ψ(t) = P (t)|Ψ is now quasi-periodic, and can be expressed as |Ψ(t) = |Ψ(ωt + θ 0 ) , where |Ψ(θ) is parametrized by the torus T m . This evolution now defines a map on T m → Ω d . If this map cannot be continuously deformed to the constant map, the evolution is nontrivial.
The question of how to classify maps T m → Ω d , or more generally maps X → Ω d for any space X , remains in principle an open problem. There are, however, good reasons [87][88][89][90][91][92][93][94] to conjecture that the answer to this question is already contained within the frameworks used to classify stationary topological phases. See Appendix E for a more technical discussion.
Here, we will focus on a special class of phases, which are natural generalizations of bosonic symmetryprotected topological (SPT) phases and already display a rich set of behaviors. Recall that equilibrium bosonic SPT phases with unitary symmetry G in d spatial dimensions are believed to be partially classified 1 by the group cohomology H d+2 (G, Z) [95][96][97][98]. We can also write group cohomology as the singular cohomology of the so-called classifying space BG of the group G: H d+2 (G, Z) ∼ = H d+2 (BG, Z). The idea is that maps X → Ω d in the presence of symmetry G should be partially classified by replacing BG → X × BG, i.e. the classification is H d+2 (X × BG, Z). In fact, under general conditions for both SPT and symmetry-enriched topological (SET) phases of bosons or fermions, the classification can be derived from the classification of equilibrium SPT and SET phases simply by replacing BG → X × BG. We give some justification for this in Appendix E.
The following powerful statement follows, using the fact that BZ m = T m . We find that the classification of maps T m → Ω d in the presence of symmetry G is in one-to-one correspondence with the classification of stationary symmetry-protected and symmetry-enriched phases with symmetry Z m × G. The interpretation of the additional Z m symmetry is that they correspond to the "multiple time-translation symmetries" referred to in Sec. III B. We call this the quasi-periodic equivalence principle. The periodic case, which we could call the "Floquet equivalence principle", was discussed in [60,86]; compare also the "crystalline equivalence principle" of Ref. [92].
The simplest case of fundamentally nonequilibrium topological phase are ones where eigenstates |Ψ are themselves in the trivial G SPT phase. In this case, the classification of maps T n → Ω d with a G symmetry imposed has the general decomposition m r=1 (k1···kr) where (k1···kr) indicates a sum over all possible choice of non-repeating numbers k 1 , · · · , k r ∈ {1, · · · , m}. Here C s is the classification of equilibrium SPT phases or invertible topological orders with G symmetry in s spatial dimensions for s ≥ 0; for s < 0 we set C s = π −s (Ω 0 ). For the case of the group cohomology classification H d+2 (Z ×m × G, Z), Eq. (35) can be proven using the Künneth formula [99], although the result holds more generally (see Appendix E). The formula Eq. (35) has a simple physical interpretation: the different terms correspond to cases where r time-translations symmetries are 'essentially' involved in the definition of the corresponding phases. One term corresponds to those maps T m → Ω d that depend on all m incommensurate frequencies, which are classified by C d−m . The remaining terms of Eq. (35) depend on only r < m frequencies, corresponding to micromotions which vary on r-dimensional hyperplanes T r ⊆ T m and are deformable to constant evolutions in the other directions.
For periodic drives (m = 1), Eq. (35) reduces to just C d−1 . In this case, we can think of the non-triviality of the micromotion as a pump per Floquet cycle, which nucleates equilibrium (d−1)-dimensional SPT phases and transports them onto the boundary of the system [60,82,83,100]. For m ≥ 2, C d−m terms in Eq. (35) can be interpreted as a higher-order pump ('pump of pumps'). However, making this notion concrete for observables on the boundary of model systems is beyond the current discussion.
Let us now focus on the micromotions which depend on all m frequencies. The ideas of Sec. III allow us to con- struct non-trivial micromotions in a preheating regime through what we refer to as a "bootstrap construction". Indeed, suppose that we choose Q = (1/N )I m for integer N , so that G = Z ×m N . For concreteness we will consider the case m = 2. Then the effective Hamiltonian D has an enhanced symmetry G × G, and accordingly it can host bosonic SPT phases protected by G × G, which are classified by H d+2 (G × G, Z). Invoking the Künneth formula, we see that this contains a factor Suppose we ensure that the effective Hamiltonian D has eigenstates which are in an SPT phase associated with the emergent G × G symmetry corresponding to an element of Eq. (36). The projection map Z → Z N induces a map in group cohomology where we have used the fact that H 1 (Z, A) ∼ = A for any finite Abelian group A. Note that H d (G, Z) obtained in Eq. (37) corresponds to the classification of SPT phases in d − 2 dimensions, i.e. C d−2 in the above notation. If the image of the element of Eq. (36) in Eq. (37) is nontrivial, it means that the micromotions corresponding to the eigenstates of D are nontrivial in a way that relies on both frequencies.

A. Example Hamiltonian
Let us construct a minimal example of a novel quasiperiodically-driven phase of matter realized through the above considerations. Consider a square lattice in two spatial dimensions with spin-1/2 particles on the links (A), on the sites (B), and on the plaquettes (C), as shown in Fig. 7. Assume that the system has a microscopic Z 2 symmetry generated by X A := l∈A σ x l , where the product is over spins on the links. We then define S B = 1 2 v∈B (σ x v + 1), where the sum is over all spins on the vertices, and S C = 1 2 p∈C (σ x p + 1), where the sum is over spins on the plaquettes. Furthermore, We consider the time-quasiperiodic Hamiltonian Here D 0 is a time-independent Hamiltonian that commutes with X A , X B and X C that we will specify explicitly below and V (t) is some generic perturbation which must commute with the microscopic symmetry X A . All the time-dependent quantities f 1 (t), f 2 (t) and V (t) are quasiperiodic with frequency vector ω = (ω 1 , ω 2 ) assumed to be large compared to the local energy scales of D 0 , V (t).
We now choose the driving functions f 1 (t), f 2 (t) to be for i = 1, 2, picking Q = (1/2)I 2 . Here F N (θ) is the smooth approximation to the Dirac Delta comb, as utilized in the DTQC example of Sec. IV B, given in Eq. (28). This choice of driving function strikes a balance between having a smooth drive, necessarily to achieve a long-lived preheating regime, and the property that the interaction Hamiltonian H int (θ) = U 0 (θ) † (D + V (θ))U 0 (θ) (see Sec. III C) has the term U 0 (θ) † D 0 U 0 (θ) satisfying, roughly, U 0 (θ) † D 0 U 0 (θ) ≈ D 0 . Accordingly, upon taking the frame-twisted high-frequency limit, we find that there is an effective Hamiltonian D = D 0 + δD in a long-lived preheating regime at high driving frequencies, where δD commutes with X A , X B , and X C , and where the local strength of δD is on the order of that of V (t).
Suppose D 0 is an MBL Hamiltonian whose eigenstates are in the SPT phase for symmetry Z ×3 2 corresponding to the non-trivial element of . Indeed, such Hamiltonians can be realized via a decorated domain wall construction [101], where the ground state is a superposition of fluctuating domain walls. The intersections of X B domain walls and X C domain walls (which precisely occur at A spins) carry a −1 charge of X A ; see Fig. 7. Concretely, we can write where the coefficients h v ,h l ,h p are chosen from some random distribution, and where p 1 (l), p 2 (l) are the two plaquettes adjacent to the link l; v 1 (l), v 2 (l) are the two vertices connected by the link l; and |{σ v }, {σ p }, {σ l } is a basis state labelled by the eigenvalues of σ z v , σ z p , σ v l for all vertices v, plaquettes p, and links l. Note that if the symmetric term δD is sufficiently weak, the effective Hamiltonian D 's eigenstates will also belong to the same SPT phase as the MBL Hamiltonian D 0 .
Following the general arguments from before, one finds that the eigenstates of D indeed undergo nontrivial micromotion, classified by the non-trivial element of H 1 (Z 2 , U(1)) = Z 2 , which also classifies SPT phases with Z 2 symmetry in 0 spatial dimensions (i.e. Z 2 charges).
We expect that the topology of these micromotions have direct observable signatures, but we will leave exploration of this for future work. Here we will merely make some speculations based on group cohomology formulae and analogies with other phenomena. Specifically, observe that from the Künneth formula [applied in a different order to the way that would give Eq. (35), for instance], the SPT discussed above arises from the term According to the usual physical interpretation of the Künneth formula, we expect that if we spontaneously break the microscopic Z 2 symmetry on the boundary, then fusing two Z 2 domain walls on the boundary leaves behind a nontrivial 0-dimensional SPT classified by H 2 (Z × Z, Z). The latter corresponds to quantized energy pumping between the two incommensurate frequencies (see below), so we can interpret this as Z 2 domain walls on the boundary binding a half-quantized energy pump.

B. Relation to previous works
Some aspects of topology in quasiperiodically driven systems have previously been studied, primarily in the context of non-interacting systems. We now briefly mention how these phenomena are related to our classification and construction.
Ref. [52] discussed "quantized energy pumping" between m = 2 different frequencies in d = 0 spatial dimensions. Let us interpret the topological invariant that was found. In Ref. [52], the limit of very low frequencies was studied, so that the adiabatic theorem ensured that the system is always in the ground state of the instantaneous Hamiltonian. The latter varies quasiperiodically, H(t) = H(ωt + θ 0 ) for some continuous function H(θ) defined on the 2-torus T 2 . Therefore, we can treat the projector onto the ground state of H(θ) as a function of θ, which defines a micromotion that can be classified according to the general framework discussed above. We find that This Z invariant is the Chern number over the torus T 2 of the Berry connection of the ground state. Note that the adiabatic limit considered is very different from the high-frequency one we considered above -in particular, there cannot be a decomposition of the time-evolution operator of the form Eq. (6) here [40]. However, this does not affect the classification. The topological invariants discussed in Ref. [54] can also be understood in a similar fashion.
One might think that the nontrivial invariant of this kind could also be realized in the frame-twisted highfrequency limit through a bootstrap construction as in Sec. V A. However, this is not the case. The reason is that in the frame-twisted high-frequency limit, the emergent symmetry is always a finite group G = Z ×2 /L . One can check that although H 2 (Z × Z, Z) = Z, the image in this group under the map H 2 (G, Z) → H 2 (Z × Z, Z) induced by the projection map Z × Z → G, is always trivial for any finite quotient G.
In Ref. [53], "Majorana multiplexing" was introduced, where a system in one spatial dimension which is quasiperiodically driven with m incommensurate frequencies, may host 2 m different kinds of boundary Majorana zero modes that can occur simultaneously, while being protected from coupling to each other. This is consistent with our general framework, since the Majorana modes are distinguished by the charge (±1) they carry under each of the m generators of the Z ×m "timetranslation symmetry" group. One can readily show that any combination of such boundary Majorana zero modes can be realized in a frame-twisted high-frequency limit in a bootstrap construction of the general form described above. In particular, this demonstrates that the "Majorana multiplexing" of Ref. [53] can be made stable to interactions, at least up to the parameterically long heating time t * , in the presence of strong disorder and driving at high frequencies.

VI. ESTIMATING THE HEATING TIME IN QUASIPERIODICALLY-DRIVEN SYSTEMS
We return to the important question of dynamics and slow heating in quasiperiodically-driven systems which underpin the existence of the nonequilibrium phases of matter discussed in the preceding sections. We first give schematic arguments within linear response that can be used to understand the stretched exponential scaling of the heating time t * of Eq. (3), by extending the discussion of Sec. II B. This analysis follows that of Ref. [31]. The scaling of t * that we obtain through these arguments will be borne out in the rigorous bounds on heating of Sec. VII.

A. Linear response arguments
As discussed in Sec. II B, the heating rate in a quasiperiodically driven system is governed by the competitition between the decay Fourier series coefficients of the driving term V n with |n|, and the fact that |ω ·n|/|ω| could become ever smaller as |n| increases. In order to estimate the heating time, we need to quantify exactly how small |ω · n| can get as a function of |n|. The key mathematical tool is the following Diophantine condition for all integer vectors n = 0, where c is a constant depending on the ratios of the frequencies ω i /ω j , but not on the overall frequency scale |ω|. It can be rigorously shown (see Appendix F) that all choices of frequency vectors except a set of measure zero obey such a condition for any γ > m − 1 (with a constant c depending on γ).
In particular, Eq. (43) holds (again, for any γ > m − 1, and with c depending on γ) when the ratios ω i /ω j are all irrational algebraic numbers -that is, they are each roots of some polynomial equation with integer coefficients, in which case it is known as the subspace theorem [102,103]. In our analysis we always assume that the frequency vectors ω obey Eq. (43). This will allow us to derive lower-bounds on the scaling of the heating time at high frequencies (recall that this high frequency limit corresponds to taking ω := |ω| → ∞ while keeping the ratios ω i /ω j fixed). For the periodically-driven case, we can use Eq. (43) with γ = 0.
In linear response, a term V n can drive transitions between energy levels of the average Hamiltonian separated by energy ∆E n = ω · n. However, such processes require a rearrangement of at least ∼ ∆E n /J sites, where J is the local strength of the average Hamiltonian, and hence the amplitudes for such processes are suppressed by a factor e −κ(∆En/J) ≤ e −κ(ω/J)|n| −γ (44) for some constant κ > 0, using the Diophantine condition Eq. (43). Since only |n| = 0 processes contribute to heating, the smallest value of |n| is 1, in which case this term evaluates to e −κ(ω/J) , an exponentially small factor at high frequencies. For γ = 0 (the periodic case), this is the whole story since there is no n dependence on the right-hand side of Eq. (44), and we recover the usual exponential scaling of heating time with frequency. For quasiperiodically driven systems, on the other hand, γ > 0 and then the right-hand side of Eq. (44) goes to 1 at large |n|. This would suggest that the heating rate is not suppressed even at high frequencies.
To derive stronger results, the key is to use the smoothness of the driving. If we assume the local strength of V n is bounded by some function g(|n|), then the rate of a heating process governed by a given V n is controlled by both Eq. (44) and g(|n|), which gives a suppression factor To obtain the full heating rate, we sum over all such processes. Asymptotically, it is governed by the fastest heating rate: For a given form of g(|n|), there will be a value of |n| that maximizes Eq. (45), and the heating time then scales like the inverse of the maximum value of Eq. (45). If g(|n|) decays exponentially in |n|, which is the case we consider throughout our paper, then we find t * ∼ exp[const.
(ω/J) 1/(γ+1) ], which is a similar scaling to the bound in Eq. (3). One sees that unlike the Floquet case, for quasiperiodic driving (γ > 0), it is crucial to use the smoothness of the driving. Indeed, if g(|n|) does not decay with |n|, we see that the heating time does not grow with frequency at all if γ > 0. Treating other forms of drives, where g(|n|) has different asymptotic forms at large |n|, is an interesting direction to develop; see Sec. VIII C.

B. Stability to varying the frequencies
Our results have assumed keeping the ratios ω i /ω j constant, and only changing the overall frequency scale |ω| to reach the high-frequncy limit. In practice, the physical time evolution of the system should not depend too sensitively on the precise ratios ω i /ω j and associated constant c. We expect there to be a long timescale before the system can resolve the distinction between a drive with driving frequencies ω and one with nearby driving frequencies ω + δω (for example, a close rational approximation, for which our results technically do not apply).
We can estimate this timescale at the level of the linear response arguments. Consider a vector ω which satisfies the conditions for our results to hold, Eq. (43), and consider perturbing the frequency ω → ω = ω + δω. By rescaling we can ensure that |ω | = |ω|. From the triangle inequality, provided that That is, the new frequency vector still satisfies the required approximation condition (with a constant c that is independent of δω), but only when |n| ≤ n max .
In this context, the effect of processes with |n| ≥ n max will never be felt if we truncate the Fourier expansion V n of the driving so that V n = 0 for |n| > n max . Assuming smooth driving, so that V n decays exponentially with |n|, then this truncation will not substantially affect the dynamics until a time t perturb ∼ e Cnmax = exp C c|ω| 2|δω| for some constant C. After t perturb , the evolution of the system will be governed by a different regime. Therefore, the worst that could possibly happen when perturbing the frequency vector is that the heating time t * will grow with frequency until it reaches t perturb and then stop growing because the system immediately heats at time t perturb . However, this heating scenario is not inevitable, because the dynamics after t perturb could themselves have a much longer heating time t * that grows with frequency. For example, if ω is a rationally related frequency vector, then the system is really periodic and one can invoke Floquet slow heating results.

VII. LONG-LIVED, PREHEATING REGIME IN QUASIPERIODICALLY-DRIVEN SYSTEMS: RIGOROUS RESULTS AND SKETCH OF PROOF
In this section we finally formalize the preceding discussions on slow heating and emergent symmetries, into a rigorous theorem on preheating in quasiperiodicallydriven systems. We do this by explicitly constructing the effective Hamiltonian D and providing bounds on its validity. The construction employed also manifestly allows for emergent symmetries to appear in D, should the driving Hamiltonian have twisted-TTSes as described in Sec. III B. The exact formulation of the proof and heating bounds is somewhat involved and we leave the technical details to Appendix G. Here we shall provide an accessible statement of the theorem, as well as an outline of the proof.

A. Conditions and setup
Recall that we are considering quasiperiodically driven systems defined on lattices in d-spatial dimensions with locally bounded Hilbert spaces, i.e. spins or fermions, with a driving frequency vector ω that is rationally independent, ω · n = 0 for any nonzero integer vector n ∈ Z m . We furthermore assume that the timequasiperiodic Hamiltonian is quasilocal, with the drive performed at high-frequencies and is furthermore smooth in time. What we mean by the high-frequency condition is that each driving frequency ω i is large compared to the local energy scales J of the Hamiltonian H(θ), and what we mean by the smoothness of drive condition is the imposition of the condition that the Fourier modes H n decay exponentially fast with |n|, i.e. H n = O(e −κ |n| ) for some κ > 0. Additionally, let us remind the reader that we will assume that the frequency vector ω obeys the Diophantine condition Eq. (43), which holds for all choices of rationally-independent frequency vectors except for a set of measure zero.
More precisely, our construction and theorem makes use of a notion of a local norm O κ parameterized by a constant κ > 0, appropriate for a many-body operator O(θ) parameterized on the torus T m and which acts on an infinite lattice with bounded local Hilbert space dimension. To define this norm, let us first write (nonuniquely) a many-body operator O(θ) in terms of a sum of local 'potentials' O Z (θ) on T m , where Z is a finite subset of the lattice and where O Z (θ) only acts nontrivially on sites x ∈ Z for all θ ∈ T m . We can furthermore decompose the local potential into its Fourier modes O Z (θ) = n O Z,n∈Z m e in·θ so that We then define O κ to be for some constant κ > 0, where the supremum is over sites x on the lattice, and the norm appearing in the sum is the standard operator norm. The norm measures the strength of local terms making up O(θ), specifically taking into account the decay of the strength of interactions in both spatial extent and Fourier space. Terms corresponding to larger spatial support and higher Fourier modes are weighted more, parameterized by κ which can be understood as the decay constant. Note that the norm is only really useful if it is finite. Thus the local terms appearing in the sum in Eq. (49) have to be decaying at least exponentially fast both in real and Fourier space. The quasilocal and time-smooth nature of the driving H(θ) we imposed ensures that there is some κ > 0 for which H κ < ∞.

B. Theorem and statement of results: long-lived, preheating description of dynamics and emergent symmetries
With this setup, we turn to the object of interest: the unitary time-evolution operator generated by a quasiperiodically-driven Hamiltonian H(t) = H(ωt + θ 0 ) for some θ 0 ∈ T m . We assume there is some decay constant κ 0 > 0 such that H(θ) has a local norm J := H κ0 which is small enough compared to the driving frequency: J ≤ Kω. Here K is some numerical constant that does not depend on the Hamiltonian, geometry of the lattice or driving frequencies, which we give explicitly in Appendix G. We then have the following statements. (A) Existence of a long-lived preheating dynamical description. It is possible to find a decomposition of the unitary time-evolution operator as: (52) where P (t) = P (ωt + θ 0 ) is a time-quasiperiodic quasilocal unitary, D is a time-independent, quasilocal Hamiltonian, and V (t) = V (ωt + θ 0 ) is a time-quasiperiodic, quasilocal Hamiltonian. The functions P (θ), V (θ) are smooth on T m .
Define the decay constant κ = κ 0 /4. The Hamiltonian D is close to the time-averaged Hamiltonian where Here C, K are numerical constants which we compute in Appendix G, that importantly do not depend on the Hamiltonian or driving frequencies.
We now spell out the dynamical consequences.
• Slow heating. The time-averaged Hamiltonian H is an almost conserved energy operator, up to perturbative corrections of order J/ω, captured precisely by for all t ≥ 0. In the above, the standard operator norm is used, and we have divided by the volume "Vol." of the system. In other words, the normalized energy density as measured by H /J grows very slowly, apart from a small quantity of order J/ω. HereK, C are finite constants that do not depend on the Hamiltonian, geometry of the lattice, or driving frequencies. An analogous statement holds replacing H → D as they only differ by small terms on the order of J/ω.
From this we can thus derive that the heating time t * , defined as the time beyond which the quantity 1 Vol. U † (t) H U (t) − H starts growing appreciably, obeys a lower bound for appropriately defined numerical constants C , C. Note that this is only a lower bound for t * , because Eq. (56) is only an upper bound on how fast the energy density can change.
• Effective description of dynamics. If we definê then for any local operator O, we have that where K (O) is a numerical constant not depending on the Hamiltonian or driving frequencies. Thus, the difference Eq. (59) is very small for times less than the heating time t * ∼ 2 q * .
(B) Existence of emergent symmetries in effective Hamilton D. If in addition the original driving Hamiltonian H(θ) has twisted-TTSes generated by gτ (for some set ofτ s), that is, The above statements form the rigorous basis on which the results of this paper rests on. They show that (i) the system does not heat until the long time t * , that (ii) in the preheating regime the decomposition Eq. (6) holds up to small corrections that can be ignored until t * , in the sense made precise by Eqs. (56,59), and that (iii) the effective static Hamiltonian D possesses emergent symmetries as a consequence of twisted-TTSes of the drive.

C. Main ideas and sketch of proof
Let us provide here the ideas underlying our technical procedure that allows us to derive the previous assertions, as well as sketch the proof. To obtain the dynamical statement (A), we employ an iterative procedure 'renormalizing' the initial driving Hamiltonian H(θ) through a series of small rotations that sequentially reduce the norm of time-dependent pieces. This is possible under conditions of high frequencies where the parameter J/ω, the ratio of the local energy scales to the driving frequency ω, sets a natural small parameter. Upon stopping at some optimal order, this will eventually give us the effective static Hamiltonian D as well as a remnant small time-dependent piece V (θ). Note that it is expected that the procedure generically cannot be carried out ad infinitum as this would imply that a driven interacting system has a static local Hamiltonian description. This would go against the unbounded heating we expect to occur at long times.
Such a logic is behind the rigorous prethermalization works of Refs. [34,35], and indeed our manipulations largely follow closely that of Ref. [34] but with a number of technical extensions to handle the quasiperiodicallydriven scenario. Our main contribution, however, and the biggest departure from the earlier works, is that we will employ a renormalization procedure specifically tailored for preserving a twisted-TTS at all stages, which then allows us to obtain our statement (B) on emergent symmetries in D.
Setting up the iteration. -The renormalization process is nothing but a sequence of well-chosen rotating frame transformations, effected by the unitaries P (q) (θ) = P (q) (θ + τ ) where τ ∈ L ∈ 2πZ m and q = 0, 1, · · · up to some cut-off q * to be determined. We start the process by defining the original Hamiltonian as the zeroth-level renormalized Hamiltonian H (0) (θ) ≡ H(θ). The rotation P (q) (θ) defines the q-th level renormalized Hamiltonian at the next level, Here H (q) (θ) is the resulting Hamiltonian in the new frame of reference. Let us also define the unitaries in particular, U (0) (t) = U (t), the unitary time-evolution operator of interest. In terms of dynamics, the iterative procedure is As mentioned, the aim is to reduce the time-dependent terms in H (q) (θ) at each level. To that end, let us define the time-averaging operation so that we can write H (q) (θ) = D (q) + V (q) (θ), where Thus we see that we need to eliminate, or at least reduce, the contributions of V (q) (θ) in the next level Hamiltonian H (q+1) (θ). The high-frequency assumption allows us to choose appropriate 'small' rotations P (q) (θ). To gain some intuition, we write P (q) (θ) = e A (q) (θ) for some antihermitian operator A (q) (θ) assumed to be smaller than V (q) by a factor of 1/ω. Then, using the Duhamel formula to expand H (q+1) (θ), we have and we see schematically that all terms beyond the first line are smaller than V (q) by a factor of J/ω or less. Let us therefore demand that in our iterative procedure, the generators A (q) (θ) are chosen to satisfy Note that there is a freedom of choice in the solution of this partial differential equation, as the initial condition has not been specified. So far, the manipulations have been formally identical to that of Ref. [34] upon reduction to the Floquet case. However (and here is the crucial technical difference), we solve the above equation for A (q) (θ) = n A (q) n e in·θ in terms of Fourier modes of V (q) (θ) = n V (q) n e in·θ , with specific choice This fixes the initial condition by setting A (q) = 0. By contrast, the immediate generalization of Ref. [34] would be to set A (q) (θ 0 ) = 0. Our different choice of initial condition is what allows our iterative procedure to preserve twisted-TTSes, as we will see later.
Eq. (68) highlights the fact that the solution only makes sense should the Fourier modes V (q) n decay fast enough so that A (q) n can be written as a convergent Fourier series and A κ is at least well-defined for some values of κ. This is the technical reason for the imposition of the 'smoothness' of drive conditions. As a consistency check, one can see from Eq. (68) that A (q) is smaller than V (q) by a factor of ∼ 1/ω, as we postulated. Note that this choice also preserves the quasilocality: if V (q) (θ) is quasilocal, then so will be A (q) (θ).
Estimating the optimal order q * . -Having set up the iteration, let us now present the logic behind bounding how far this iterative procedure can be carried out to. While satisfying of the relation Eq. (67) makes H (q) (θ) ever less dependent on θ, there is a price to pay: the extra terms generated at each level in the renormalized Hamiltonians are of ever longer range. To account for this and in order to meaningfully estimate their local strength, we allow for some sequence of strictly decreasing decay constants κ 0 > κ 1 > κ 2 · · · > 0 which we have to pick judiciously, and measure the Hamiltonian H (q) (θ) at the q-level through its norm H (q) (θ) κq . However, at some stage, the smallness of κ q will impede our ability to bound the Hamiltonian H (q+1) (θ) at the next level; this is when the renormalization procedure stops. The aim is to choose a suitable set of decay constants κ q , which allows for the inductive process to be carried out to as high an order as possible, rendering the resulting V (q * ) (θ) optimally small and giving the effective Hamiltonian D := D (q * ) at the stopping order q * .
Our proof indeed follows such a procedure. We skip the heavily technical details in setting up various inductive bounds, as well as choosing the appropriate decay constants κ q , but merely state that we end up with the optimal level of truncation q * as given in Eq. (55). This also yields the claimed bounds on D, V (θ), defined as the optimal D (q * ) , V (q * ) (θ), respectively. The result (52) then follows, with P (θ) := q * l=1 P (l) (θ). We refer the reader to the Appendix G for the full details of our manipulations.

D. Emergent symmetries in effective Hamiltonian D
Thus far, the above discussion was purely dynamical, as related to statement (A). Now let us consider the statement (B) on the emergent symmetries in D.
Our choice of solutions Eq. (68) explicitly preserves the twisted-TTSes of the original driving Hamiltonian H(θ) at every level of the renormalization procedure. This in turn allows for them to be manifested as unitary operators that commute with the effective Hamiltonian D. To see this, recall that a twisted-TTS acts in Fourier space on a Hamiltonian according to gτ H n g † τ = e in·τ H n , so in particular H = H 0 = gτ H 0 g † τ is symmetric. Now if the q-th level Hamiltonian H (q) (θ) has a twisted-TTS, then V (q) (θ) has the same symmetry, while D (q) = H (q) commutes with gτ . But from Eq. (68), this immediately implies that that A (q) (θ) and hence P (q) (θ) will also inherit the twisted-TTS of V (q) (θ), and therefore so does H (q+1) (θ). As this is true for every q, we end up with the statement that the effective Hamiltonian D obtained at the optimal order q * has emergent symmetries gτ .

VIII. EXTENSIONS & FUTURE DIRECTIONS
Having demonstrated how one can achieve novel nonequilibrium phases of matter in quasiperiodicallydriven systems, we now consider extensions and future directions arising from our work. Many of these are directly inspired by the recent development in nonequilibrium Floquet systems.

A. Long-range interactions
One of the assumptions which allowed us to bound heating rates and establish the pre-heating regime, was that the many-body systems we considered had quasilocal interactions -that is, the amplitudes of the interaction terms in the Hamiltonian decay at least exponentially fast with space. The study of prethermalization in Floquet systems suggests that this restriction may be lifted to encompass long range interactions. Refs. [104,105] have demonstrated slow heating for periodic driving in the presence of two-body power-law interactions, provided that the interactions decay with distance as ∼ 1/r α , with α > d where d is the spatial dimension, see also Ref. [106]. The existence of an effective Hamiltonian approximately generating the dynamics in the preheating regime, however, can only be proven for α > 2d. We can immediately combine our proof with the approach of [105] to derive similar results of slow heating and emergent symmetries in quasiperiodically driven systems with power-law interactions. These extensions are particularly valuable for realizing quasiperiodic phases in trapped ion systems, which naturally have long range interactions.

B. Time-independent systems: continuous time quasi-crystal
Although we have focused on applications to quasiperiodically-driven systems, our theorem also has an intriguing implication for systems with time-independent Hamiltonians. In particular, let Γ 1 , · · · , Γ m be a set of commuting operators, each of which has integer eigenvalues, and let ω = (ω 1 , · · · , ω m ) be some rationallyindependent numbers. Consider the time-independent Hamiltonian where Observe that the time evolution generated by H 0 is quasiperiodic U 0 (t) = exp(−iH 0 t) = exp(−itω i Γ i ), since U 0 (t) = U 0 (ωt), where U 0 (θ) = exp(−iθ i Γ i ). If we define the Hamiltonian in the rotating frame generated by H 0 , i.e.
then it also has quasi-periodic time dependence H int (t) = H int (ωt), with We invoke our theorem to construct a rotating frame transformation P (θ) that generates a time-independent quasi-local effective Hamiltonian D, up to corrections that can be ignored until a time t * that scales like a stretched exponential in |ω|. This might not seem very useful, since the original Hamiltonian H was already a time-independent Hamiltonian. However, crucially, H int has a twisted TTS in the sense discussed in Sec. III B. Indeed, we have that but now for any vector τ ∈ R m (thus, we have a continuous twisted TTS rather than a discrete twisted TTS as considered previously). Therefore, our theorem ensures (see Eq. (52)) that the time-evolution operator generated by H int (t) can be decomposed (ignoring corrections that only become important at times t t * ) as where D has the emergent symmetries [D, g τ ] = 0 for all vectors τ , and hence (by taking the limit τ → 0) [D, Γ i ] = 0 for all i. Moreover, P (t) = P (ωt), where P (θ) obeys the twisted TTS, i.e. P (θ+τ ) = g τ P (θ)g † τ . Hence, using the form of g τ , we find that P (θ) = U † 0 (θ)P U 0 (θ), where P := P (0). Substituting into Eq. (74), we find On the other hand, we also know that, by definition, U int (t) = U † 0 (t)e −iHt . Therefore, we find that The right-hand side of Eq. (77) commutes with P Γ i P † for i = 1, · · · , m. That is, the Hamiltonian H has an approximate (because it only holds up to time t * ) emergent U(1) ×m symmetry. The m = 1 case was already proven in Ref. [12] (though the connection with twisted TTS was not identified). Let us discuss two applications of this result. Firstly, we can imagine that the emergent U(1) ×m is spontaneously broken. In that case, we obtain a continuous time quasicrystal -a time-independent Hamiltonian which spontaneously develops a quasiperiodic response with frequencies ω 1 , · · · , ω m until the long time t * . This is a generalization of the prethermal continous time crystal discussed in Ref. [12].
Another application is the topological protection of quantum information. Ref. [107] used the m = 1 version of this result to argue that the decoherence time for a qubit encoded in a Majorana zero mode on the boundary of a one-dimensional topological superconductor could be made exponentially long (in a parameter of the Hamiltonian) at arbitrary energy density, not just in the ground state, even without disorder. Ref. [107] also considered a two dimensional planar code, but they were not able to show long lifetime in the limit of large system size, because the m = 1 result does not allow for the seperate conservation of the number of e-type and m-type excitations separately. With our new m > 1 result, we can now ensure that these numbers are separately conserved, which does lead to a (stretched) exponentially long lifetime for the encoded qubit, even as the system size goes to infinity.

C. Non-smooth drives
In this paper, we have focused on drives that are smooth in time. It would be interesting to further develop our discussions for cases where the Fourier modes of the drive H n decay slower than exponentially with |n|. For example, discrete step-drives decay as some powerlaw with |n|. Within the simple linear response estimates of Sec. VI, it seems natural to expect a heating time t * that scales like a power-law with frequency. However, a more careful calculation of the heating processes along the lines of Sec. VII would be valuable.
Ref. [55] considered driving an MBL system with a discrete Fibonacci step-drive and found two slow heating regimes, one governed by an effective Hamiltonian generated from a high-frequency expansion whose description lasts for power law times, and a subsequent one where there is a slow logarithmic decay of observables that lasts up to an exponentially long timescale. Some of the phenomenology considered is very different from the one discussed in Sec. II D; establishing the relationship to the work here remains an open question.
On a technical level, we note that extending our results of Sec. VII to the case of power-law decay of Fourier modes is more complex then applying them to the case of power law decay of interactions with distance. Considering our results and those of [105], we see that the analog of power-law dependence of Fourier modes are interactions that decay as a power-law in the number of sites contained within the support of the interaction. By contrast, Ref. [105] established results only for interactions which decay as a power-law in the diameter of the support, while still decaying exponentially in the number of sites.

D. Quasiperiocally-driven topological phsaes beyond eigenstate micromotions
In Sec. V, we focused on eigenstate micromotions as a diagnostic of nontrivial topological phases in periodically or quasiperiodically-driven MBL systems. However, it is known that there are topological Floquet phases which cannot be diagnosed in this way, many of which have particularly fascinating phenomenology. For example, the phases discussed in Refs. [108][109][110][111] are characterized by a chiral pumping of quantum information along the onedimensional boundary of a two-dimensional bulk. One expects analogous phenomena in quasiperiodically-driven systems, but we leave such developments for future work.

E. Experimental realizations
Lastly let us briefly mention the possibility of experimentally realizing the quasiperiodically-driven phases of matter we have discussed, in particular, the discrete time quasi-crystal (DTQC) phase described in Sec. IV.
The driving protocol and conditions required, described in III C, can easily be realized in setups such as in synthetic quantum systems of trapped ions, or in solid state systems like nitrogen-vacancy (NV) defects in diamond. Both cases give rise to ensembles of coherently in-teracting effective spin degrees of freedom which are well isolated from the environment. Appropriate sequences of laser or microwave pulses can be engineered to realize particular drives in these systems. Indeed, these were utilized to realize sequences [11,13] that led to signatures of discrete time crystal phases observed in experiments of a 1d chain of trapped ions [15] and in a dense, disordered 3d ensemble of NV centers in diamond [16].
Extending these protocols to realize time-quasiperiodic phases is relatively straightforward, although the timedependence of the drives needs to be smooth in time to avoid fast heating. The different platforms offer different comparative advantages. In trapped ions, the ability to tune the range of long-range interactions between spins implies that it should be possible to realize a DTQC in a prethermal regime even in 1d, provided that the power law exponent satisfies 1 < α < 2. On the other hand, one can utilize the fact that NV centers in diamond are grouped by their orientation with respect to the crystallographic axes of diamond. Using interactions within and between multiple groups, one can naturally realize more complicated TQCs, such as the example Hamiltonian in Appendix. D 2 giving rise to a Z 2 × Z 2 DTQC, which utilizes two collection of spins. Note that if in 3d, the dipolar interacting nature of the spins precludes localization, but thermalization is nonetheless slow in heavily disordered samples (precisely, critically slow, see [14,112,113]). Thus, a long-lived version of a DTQC might still be realizable, protected by critically slow thermalization dynamics.
The dynamical signatures that one would look out for, would be the appearance of peaks at fractional harmonics of the input drive frequencies in the power spectrum of a local observable. Additionally, these signatures should be robust to small perturbations to the drive protocol which still preserve its time-quasiperiodicity. This would signal the spontaneous breaking of multiple time-translation symmetries of the original drive and is the defining characteristic of the DTQC phase (see Sec. IV). The exact pattern that is manifested in the Fourier harmonics, however, depend on the exact symmetry breaking pattern, realized using different driving protocols. As there are multiple time-translation symmetries in quasiperiodic drives, there are myriad complex symmetry breaking patterns, which could be observed even for a given experimental platform.

IX. CONCLUSION
In this paper, we have shown how interacting quantum many-body systems that are quasiperiodically-driven can realize a panoply of long-lived nonequilibrium phases of matter under a general set of driving conditions. These phases are protected by multiple time-translation symmetries, which arise as the driving Hamiltonian derives from a function H(θ) defined on a higher dimensional torus. They are fundamentally different from phases re-alizable in static or Floquet systems and add to the richness of the landscape of possible phases of matter realizable in nonequilibrium settings. As an exemplar, we described the phase of matter obtained by the spontaneous breaking of some or all of the multiple time-translational symmetries -a time quasi-crystal. We also gave a classification of the symmetry-protected topological phases of matter achieved in the quasiperiodic setting.
Key to our results was our ability to identify a dynamical regime, a so-called preheating regime, in which the deleterious effects of heating in driven systems was controlled for a parameterically long time, as well as our analysis of how multiple-TTSes play out in this regime. We provided a class of general driving conditions that realize prethermal or quasiperiodic many-body localized phases in this regime. We emphasize that their dynamical signatures are universal and robust against small perturbations to the driving protocol, as long as they respect the time-quasiperiodicity of the drive: they are thus genuine nonequilibrium phases of matter in interacting many-body quantum systems.
Our results open up exciting new lines of experimental and theoretical research. Indeed, the nonequilibrium phases of matter we have discussed, in particular, time quasi-crystals, are immediately directly accessible in experimental platforms of today. Theoretically, our work establishes new universal structure present in quantum many-body systems in highly out-of-equilibrium settings. What other possible fundamentally nonequilibrium phases of matter remain to be discovered?
We note two aspects of the following discussion: 1) This expansion is not the naive direct generalization of the Floquet-Magnus expansion to quasiperiodicallydriven systems, a point we shall highlight in the derivation. 2) We do not prove any bounds on the validity of the expansion. More precisely, what this means is as follows. As with the usual Magnus expansion in Floquet systems, we do not expect the series that will write down to be convergent for a many-body system -instead, it is an asymptotic series. For the Floquet case, it is possible to analyze the optimal order to which the series should be truncated to [32,33,35], which gives an optimal effective Hamiltonian D and an estimate on the lifetime that the system can be viewed as having dynamics under D. However, we do not give such an estimate, but merely generate the formal series. That said, the analysis of this series does suggest an interesting alternative and complementary derivation of our results, which we reserve for future work.
The logic behind the generalized Floquet-Magnus expansion, or for that matter, the usual Floquet-Magnus expansion in periodically-driven systems, simply involves moving into a suitable rotating frame of reference and deriving the rotating frame Hamiltonian. At high frequencies ω = |ω|, it is possible to organize this transformation as a perturbative series in 1/ω, and use the freedom endowed by the frame transformation to cancel time-dependent pieces in the rotating frame Hamiltonian, order by order. This is also the same logic underlying our technical procedure that yields a different effective Hamiltonian D; thus, the following expansion should be viewed as a different derivation of an effective description of the system.
To start, we consider a quasiperiodically-driven Hamiltonian H(t) = H(ωt + θ 0 ), where H(θ) is a periodic Hamiltonian on the torus. The unitary time-evolution operator obeys the following equation We move into a new frame of reference effected by the unitary P (t) = P (ωt+θ 0 ), for some periodic P (θ) on the torus. In other words, we decompose the time evolution operator as U (t) = P (t)Ũ (t)P (0) † , whereŨ (t) obeys We see that time-evolution is now generated by the Hamiltonian D(t), which is time-quasiperiodic as D(t) = D(ωt + θ 0 ) for a periodic D(θ) which is given explicitly as We will show how P (θ) can be chosen to make D(t) static, i.e. D(t) = D.
To that end we write P (θ) as the exponential of a sum of antihermitian operators, , where we have organized Ω as a series with terms labeled by q. q will turn out to track the order of the high frequency expansion. Now we can write, using Duhamel's formula, Expanding, we get The question now is how best to reduce the θdependent pieces of D(θ). One natural way, is to assume that Ω q (θ) has norm ∼ 1/ω q , i.e. treat it as a high frequency expansion. Therefore we need to pick the generators Ω q+1 (θ) in such a way that this assumption holds at every step. Indeed, if we write G (q) (θ) in terms of a Fourier series then we see that by imposing the condition n =0 that the norm of Ω q+1 is smaller than Ω q by 1/ω, since Ω q+1 ∼ G (q) /|ω| and G (q) ∼ 1/ω q , schematically. The imposition of Eq. (A8) thus defines a family of expansions, differing by a choice of initial condition. Indeed, the immediate, direct generalization of the Floquet-Magnus expansion to quasiperiodically-driven systems corresponds to one particular choice: One solves Eq. (A8) with the condition that Ω q (ωt + θ 0 )| t=0 = 0, so that P (t)| t=0 = P (ωt + θ 0 )| t=0 = I. This gives the solution To see why this is the 'standard' expansion, for a Floquet system we get from the above P (0) = P (T ) = I and so at stroboscopic times t = N T for N ∈ Z we have the familiar relation U (N T ) = e −iDN T , ignoring questions of convergence. Instead, we will choose to solve Eq. (A8) as This uniquely defines our 'generalized Floquet-Magnus' expansion.
Lastly, let us explain how our choice of solution Eq. (A10) giving our generalized Floquet-Magnus expansion manifestly allows for twisted time-translation symmetries to be preserved: it is such that if the original driving Hamiltonian H(θ) had a twisted time translational symmetry, that is, if there exists a unitary operator g and integer N satisfying g N = I, as well as some vector τ on the torus such that H(θ +τ ) = gH(θ)g † , then D(θ) also possesses the same twisted time-translation symmetry. To see this, recall that on the Fourier modes a twisted TTS acts as gH n g † = e in·τ H n . Now, if we assume G (q) (θ) has a twisted-TTS for some q, then our choice of Ω q+1 (θ), Eq. (A10), will also inherit the same twisted-TTS (one just sees that the Fourier modes of Ω q+1 (θ) obey the required transformation). Consequently, G (q+1) (θ), given by sums of conjugations of H(θ) with various Ω k (θ) with k < q, etc. (see Eq. (A6)) will also have the same twisted-TTS. Note the standard Floquet-Magnus expansion corresponding to a choice Eq. (A9), does not have this property.
Let Q be some r × m matrix with rational entries. Let L be the lattice comprising all vectors of the form 2πn for some integer vector n. We are interested in determining the sublattice L ≤ L comprising those vectors τ ∈ L such that e iQij τ j = 1.
To solve this problem, we invoke the well-known fact that any integer matrix Z has a Smith normal form decomposition Z = V DW where V and W are unimodular integer matrices (that is, they are invertible integer matrices whose inverses are also integer matrices) and D is a diagonal integer matrix. Since any rational matrix can be converted to an integer matrix simply by multiplying by some integer, we conclude that there is also a Smith normal form for rational matrices: Q = V DW where V and W are still unimodular integer matrices, and D is a diagonal rational matrix.
The Smith decomposition allows us to reduce the problem to the case where Q is diagonal. Suppose that if Q = D = diag(p 1 /q 1 , · · · , p min(r,m) /q min(r,m) ), where each p i ,q i are coprime integers (with q i positive and p i non-negative). We can choose to set q i = 1 if p i = 0, and we also define q i = 1 for r < i ≤ m. Then, we see that L comprises integer linear combination of the m vectors spanned by τ (i) = q i e i (no summation), for (e (i) ) j = δ ij . More generally, if Q = V DW then L comprises integer linear combinations of the vectors W −1 τ (i) . We also find that L/L ∼ = Z q1 × · · · × Z qm .
Note that a particularly simple case is where r = m and Q = Z −1 for some (not necessarily unimodular) integer matrix Z. In that case, one finds that all the p i 's are 1, the q i 's are the diagonal entries of the Smith decomposition of Z, and L simply comprises integer linear combinations of the columns of Z.

Appendix D: More examples of DTQC
We provide here the quasiperiodically-driven Hamiltonians realizing the Z 2 × Z 2 and Z 3 × Z 2 DTQCs shown in Figs. 5, 6 of Sec. IV.
The model in consideration is a direct extension of the one considered in Sec. IV B, but now involves two groups of spins. Concretely, suppose we have a system comprised of two subsystems A, B. Each subsystem consists of spin-1/2 degrees of freedom placed along either a 1d chain or a 2d square lattice, and we assume that the two subsystems A, B are stacked on top of each other (so that site i of subsystem A is directly above site i of subsystem B).
Consider the following m = 2 quasiperiodic Hamiltonian where now Similarly to the Hamiltonian in Eq. (26), Eq. (D2) describes pairwise Ising interactions of spins between and within subsystems, and a longitudinal field in the zdirection with strength h ω i /2. The interactions have strengths J ij , J ij , and we consider the 2d case where J ij = −Jδ i,j where i, j represent nearest-neighbor pairs of sites on one subsystem, and J ij = −Jδ i,j . We take 0 < J < ω i /2.
The driving Hamiltonian is taken to be and we choose Q = 1 One can understand this driving protocol as two Floquet drives, performed at frequencies ω 1 , ω 2 , acting on different halves of the system. Indeed, if the two subsystems were disjoint, this protocol would simply result in two independent prethermal discrete time crystals. In the present case, the system has interactions between the different subsystems -neverthless, we will show that a stable DTQC phase will be realized.
Repeating the same analysis as in Sec. IV B, the leading order effective Hamiltonian is We observe the following salient features: First, D (0) is manifestly Z 2 × Z 2 symmetric. In particular, the σ z i,A σ z i,B interactions between the two subsystems, which are odd under both Z 2 groups, have been eliminated. In fact, any interactions that remain between the subsystems, even at higher orders, must be Z 2 × Z 2 symmetric. Second, the Ising interactions in the z-direction dominate, in both subsystems.
Thus an initial state that has sufficiently low energy density with respect to D (0) will equilibriate to a thermal state ρ β which spontaneously breaks both Z 2 × Z 2 Ising symmetries. ρ(θ) will then have periodicity on the lattice L SSB = L , whose reciprocal lattice L * SSB is generated by the vectors α 1 = (1/2, 0) and α 2 = (0, 1/2). The power spectrum of the expectation value in time of local operators measured in this state, will then exhibit peaks at frequencies which is the manifestation of the spontaneous breaking of the TTSes of the driving Hamiltonian.
In the high frequency limit, we can treat V(t) = I and derive analytic expressions for the measurement of an observable in time s(θ) = Tr[ŝρ(θ)]. For example, for the observableŝ = 1 Vol. Vol.
i σ z i,B , this works out to be where G N (θ) is the function on the torus for which G N (ωt) = t 0 dt ω1 2 F N (ω 1 t ) and G N (θ) is the function on the torus for which G N (ωt) = t 0 dt ω2 2 F N (ω 2 t ) . We plot this observable in Fig. 5 using Tr[σ z i,α ρ β ] = 0.9 and Tr[σ y i,α ρ β ] = 0.

Z3 × Z2 DTQC
Next we present the Hamiltonian for the Z 3 × Z 2 DTQC. Consider, like in the previous example, a system comprised of two subsystems A, B of square lattices in d = 2, stacked so that sites lie on top of each other. Assume now however the local degrees of freedom of subsystem A (which reside on the sites) to be comprised of three levels (i.e. a 'spin-1') and that of subsystem B to be comprised of two levels (i.e. a 'spin-1/2'). Take the m = 2 quasiperiodic Hamiltonian to be with where the operators µ, η which act locally on A are given explicitly by These are the so-called 'clock operators' that enter into quantum clock models: η increments the state of a 'clock' (the levels which the matrix µ is diagonal in) one step in a cyclic fashion. Thus, V describes ferromagnetic Isinglike interactions between degrees of freedom of A, and also ferromagnetic Ising interactions between degrees of freedom of B. We take the driving Hamiltonian to be and Q = 1 3 0 and assume the driving frequencies are sufficiently large compared to all local couplings. We now take the frame-twisted high-frequency limit. The interaction Hamiltonian H int (θ) has periodicity on the lattice L generated by the translation vectors τ 1 = 2π(3, 1) and τ 2 = 2π(0, 2). Therefore the symmetries that the effective Hamiltonian D has, belong to the group G = Z 6 ∼ = Z 3 ×Z 2 . Viewing G as Z 3 ×Z 2 , it can be checked that group is explicitly generated by the unitary symmetries gτ = i η i forτ = 2π(1, 1) satisfying ( i η i ) 3 = I and gτ = i σ x i forτ = 2π(0, 1) satisfying ( i σ x i ) 2 = I. Repeating the same analysis as before, the leading order effective Hamiltonian D (0) can be computed to be where δ is a small term (compared to the dominant interactions above), given explicitly by It can be explicitly checked that D (0) is Z 3 × Z 2 symmetric, and that the dominant terms are the interactions µ † i µ j on subsystem A and Ising terms σ z i σ z j on subsystem B. Therefore, we expect that the steady state ρ β of this system, in the preheating regime, should exhibit spontaneous breaking of both the Z 3 ×Z 2 symmetries. ρ(θ) then has periodicity on the lattice L SSB = L , whose reciprocal lattice L * SSB is generated by the vectors α 1 = (1/3, 0) and α 2 = (−1/6, 1/2). The power spectrum of the measurement in time of local operators at long and late times, will then exhibit peaks at frequencies with n 1 , n 2 ∈ Z.

Appendix E: Topological phases
As we mentioned in Sec. V, the idea behind classifying quasiperiodic topological phases is to look at the nontrivial micromotion of the eigenstates. First we make the following observation [87][88][89][90][91][92][93][94]: all the existing classifications of topological phases in d spatial dimensions with unitary symmetry G (for example, Refs. [98,[114][115][116][117]) can be expressed in terms of classifying the homotopy classes of maps BG → Θ d for some space Θ d . Here BG is the classifying space for the group G, which is defined as BG = EG/G, where EG is any contractible space on which G acts freely (BG is unique up to homotopy equivalence). In the case that Θ • forms a so-called "Ωspectrum", this amounts to the assumption that the classification is given by some generalized cohomology theory evaluated on BG [91,93], which is the case for all known classifications of invertible topological phases. However, our observation is more general than that, and applies also to classifications of non-invertible phases, for example the G-crossed modular tensor category classification of symmetry-enriched topological phases in two spatial dimensions [88,116,118,119]. Now in all these classifications, Θ d is some abstract space that can be introduced in a formal mathematical way. However, what we want to posit is that this space actually has physical meaning. Specfically, we want to say that Θ d is actually homotopy-equivalent to Ω d , the space of all gapped ground states, introduced in Sec. V. Indeed, Kitaev (Appendix F of [89]) has shown at the microscopic level that any G-symmetric gapped ground state in d spatial dimensions gives rise to a map BG → Ω d , which would explain the origins of the maps BG → Θ d described above.
If these conjectures hold, then it is clear how to classify maps X → Ω d in the presence of symmetry G. First of all, we use the approach of Kitaev to turn this into a map X × BG → Ω d , which is equivalent to a map X × BG → Θ d . In other words, we simply replace BG → X × BG in the formulation of the classification.
We are now in a position to show how homotopy classes of maps T m → Ω d are related to SPT/SET phases with symmetry Z ×m . Indeed, this follows from the mathematical fact that as can be seen by setting E(Z m ) = R m , with Z m acting freely by discrete translations. More generally, in the case of systems with symmetry (other than time-translation symmetry) G, then we want to classify maps T m ×BG → Ω d , and then these are in one-to-one correspondence with topological phases with symmetry G×Z m , as can be seen from the equation Now let us consider cases where the states under consideration are invertible states; that is, they do not admit fractionalized excitations such as anyons. In that case, as we have already mentioned, it is believed that the spaces Θ d for different dimensions are related in a particular way. Specifically, the idea is that, if we define for any space X , h q (X ) to be the homotopy classes of maps X → Θ q , then h • (−) should define a generalized cohomology theory, which obeys some set of axioms [120]. For any generalized cohomology theory, one can show, using the suspension and wedge axioms, and the fact that [120] ΣT d = Σ(S 1 ∨ T m−1 ∨ ΣT m−1 ) (here Σ denotes suspension and ∨ denotes wedge sum), By applying this formula recursively we derive Eq. (35).
Proof. This is a higher-dimensional version of the famous Dirichlet approximation theorem [103]. Like that theorem, it can be proven using Minkowski's theorem [103], which states that if C is some bounded convex subset of R m , whose volume satisfies µ(C) > 2 m , then C ∩ Z m contains at least one point other than 0.
Indeed, using the rescaling invariance of σ to set |ω| = 1, let us define C R,c = {x ∈ R m : |x · ω| ≤ c, |x| ≤ R}, then µ(C R,c ) ≥ cv m R m−1 for some constant v m which only depends on m. Hence we can satisfy the condition to apply Minkowski's theorem to C R,c provided that cv m R m−1 > 2 m . This means that for any R > 0, there exists n ∈ Z d nonzero with |n| ≤ R such that |ω · n| ≤ 2 m+1 v −1 m R −(m−1) . This indeed implies that we must have σ(ω) ≥ m − 1.
In fact, this lower bound is generically saturated, in the following precise sense.
Proof. We extend the proof given for the case m = 2 in Ref. [121]. By the rescaling invariance of σ(ω), it is sufficient to prove that µ(S) = 0, where S is the intersection of W with B, the unit ball in R m .
Finally, let us give a concrete example of a frequency vector ω with σ(ω) = m − 1. Indeed, we have Proposition 3. Let ω be a vector of frequencies that are not rationally related, and suppose that all ratios ω i /ω j are algebraic numbers, i.e. each one is a root of some polynomial equation with integer coefficients. Then σ(ω) = m − 1.
Proof. This is known as the subspace theorem [102,103]. In this section, we provide the full, rigorous formulation of our statements on a long-lived, preheating dynamical description of quasiperiodically-driven systems. Our starting point is the iterative procedure described in Sec. VII, which we will carry out up to some optimal iteration order q * , to be computed.We will develop here the technical machinery needed to prove rigorous bounds on this iteration, and thereby prove our theorem.

Some definitions
We adopt the notion of a potential from Refs. [], which is a formal sum where the sum is over subsets Z of the lattice. In order to analyze quasi-periodic driving, we want to consider potentials that are parameterized by the torus T m , that is, Φ(θ) where θ ∈ T m . It will be prove convenient to analyze the θ dependence in Fourier space. Accordingly, we will define a "colored potential" to be a formal sum Φ(θ) = Z,n Φ Z,n e iθ·n . (G2) where n ∈ Z m represents a Fourier mode. We can define the formal commutator of colored potentials according to where we take the commutator on the right-hand side to be supported on Z ∪ Z whenever Z and Z are nondisjoint (otherwise the commutator is zero). In other words, (G4) We can simplify our notation by defining a colored set to be a pair Z = (Z, n). We define a "∪" operator on colored sets according to (Z 1 , n 1 ) ∪ (Z 2 , n 2 ) = (Z 1 ∪ Z 2 , n 1 + n 2 ), and we declare that two colored sets (Z 1 , n 1 ) and (Z 2 , n 2 ) are "disjoint" if and only if Z 1 and Z 2 are disjoint. Then, Eq. (G4) simply becomes [Φ, Θ] Z = Z1,Z2:Z1∪Z2=Z Z1 and Z2 non-disjoint which is formally identical to the uncolored case. As in the uncolored case, we define the exponentiated action of one potential Θ on another Φ according to where ad Θ Φ := [Θ, Φ]. We also can similarly write the norm Eq. (50) in this succinct notation as where the supremum is over sites x on the lattice, and where we made the following (purely formal) definitions: |(Z, n)| = |Z| + |n|, and x ∈ (Z, n) if and only if x ∈ Z.

Statement of the theorem
Our starting point is the iterative procedure described in Sec. VII C. We can perform this iterative procedure formally at the level of colored potentials. The goal is to prove bounds on the iteration, which are encapsulated in the following theorem. We start the iteration from some colored potential H (0) (θ) = D (0) + V (0) (θ), where D (0) is constant (that is, D (0) Z,n = 0 for n = 0) and V (0) (θ) has zero constant component (that is, V Z,0 = 0). We perform the iteration described in Sec. VII C, giving at the q-th iteration a colored potential H (q) (θ) = D (q) + V (q) (θ), where D (q) is constant and V (q) (θ) has zero average on T m . Then the theorem is as follows. Theorem 1. Suppose that the driving frequency vector ω obeys a Diophantine condition |n · ω| ≤ c|n| −γ |ω| for all integer vectors n where γ > m−1, and c is a constant depending on the ratios ω i /ω j but not on the overall scale of ω = |ω|.
We assume also that the norm of the driving frequency ω is large enough compared to some local energy scale ω 0 , namely there exists a decay rate κ 0 > 0, such that where is the maximum iteration order. Here where the definition of the norm · κ0 is given in Eq. (G7) and K, K numerical constants that do not depend on D (0) , V (0) , the geometry of the lattice, or the driving frequencies, and are given by with a = c −1 γ γ e −γ × max{9 × 2 6+γ , 3 × 2 1+γ κ 0 }. At the q * -th iteration define D := D (q * ) and V := V (q * ) . Then we have for some numerical constant C.

The proof
Our key technical tool is the following: Lemma 1. Let Θ, Φ be colored potentials and assume that 3 Θ κ ≤ κ − κ , with 0 < κ < κ. Then and Proof. This is the colored potential version of Lemma 4.1 from [34], and the proof follows line by line identically to the proof of that Lemma.
Armed with this Lemma, we can now proceed to analyze the iteration in a similar manner to what Ref. [34] did for the periodic case. The key new aspect compared to the periodic case is that the formula for A (q) in terms of V (q) , Eq. (68), involves a denominator which could potentially become very small. However, the Diophantine condition assumption in the theorem allows us to bound this denominator. This gives the following lemma. Lemma 2. For 0 < κ < κ, and ω obeying the Diophantine condition |n · ω| ≤ c|n| −γ |ω| with γ > m − 1 for all integer vectors n ∈ Z m and a constant c, we have where c = c −1 γ γ e −γ .
Let us suppose there is a sequence of strictly decreasing decay constants κ 0 > κ 1 > κ 2 > · · · > 0 (we will define this later), we then have where w q = 2 W (q) κq+1 . Using Lemma 1 and Lemma 2, we see that, provided that then we obtain where we have chosen 0 < κ q+1 < κ q < κ q . In particular, we choose to let κ q = (κ q + κ q+1 )/2 and so we obtain (G32) while the requirement Eq. (G30) becomes Now we turn to the question of how we should choose the decay constants κ q . A good way to do this is as follows (generalizing the approach of Ref. [34], which was an improvement on the original analysis of Ref. [34]). For notational brevity in what follows, we define λ = ω 0 /ω, the small parameter.
Armed with this intuition, we define the decay rates κ q we shall employ, as the values obtained from evaluating the following generalized function κ(q) at integer values q ≥ 1 κ q := κ(q) for q = 1, 2, 3, · · · where κ(q) µ+1 := κ µ+1 1 − (µ + 1)b µ λ µ (q − 1) for q ∈ R. (G40) We remark that the decay rate is only physically meaningful and mathematically useful when κ q > 0, and we will find conditions later on that restrict us to this range. Note also the different constant b from a of Eq. (G39), which will allow us to correct for the fact that we actually want to satisfy the finite difference equation Eq. (G34) rather than the differential equation Eq. (G37). We then find that κ(q) satisfies the differential equation with a replaced with b, i.e.
To ensure that the induction hypothesis on D (q) κq remains satisfied, let us note that Eq. (G34) also ensures that and summing over q ensures D (q) κq ≤ 2ω 0 .
This concludes the proof of the Theorem.